Matheuristics and Column Generation for a Basic Technician Routing Problem

: This paper considers a variant of the Vehicle Routing Problem with Time Windows, with site dependencies, multiple depots and outsourcing costs. This problem is the basis for many technician routing problems. Having both site-dependency and time window constraints lresults in difﬁculties in ﬁnding feasible solutions and induces highly constrained instances. Matheuristics based on Mixed Integer Linear Programming compact formulations are ﬁrstly designed. Column Generation matheuristics are then described by using previous matheuristics and machine learning techniques to stabilize and speed up the convergence of the Column Generation algorithm. The computational experiments are analyzed on public instances with graduated difﬁculties in order to analyze the accuracy of algorithms for ensuring feasibility and the quality of solutions for weakly to highly constrained instances. The results emphasize the interest of the multiple types of hybridization between mathematical programming, machine learning and heuristics inside the Column Generation framework. This work offers perspectives for many extensions of technician routing problems.


Introduction
If mathematical programming and especially Mixed Integer Linear Programming (MILP) are powerful frameworks for modeling a vast diversity of N P-hard combinatorial optimization problems, including complex real-world optimization problems, the resolution with exact methods such as the Branch and Bound (B&B) algorithm is limited in practice for large instances of real-world applications [1]. Metaheuristics find good solutions of large optimization problems without optimality or quality proofs. Highly constrained instances of optimization problems are a major bottleneck of metaheuristics. Metaheuristics and mathematical programming approaches have complementary advantages and drawbacks, and it results in their hybridization, called matheuristics [2,3]. Some matheuristics use exact methods to solve large neighborhoods of an exponential size, which tends to have fewer local minimums of a better quality than by using small neighborhoods [2,4]. Solving the exact neighborhoods allows studying the quality of implied local minimums; such intermediate results are of interest for selecting which types of neighborhoods have to be carefully implemented in a local search heuristic [4]. Matheuristics allow designing many variants of heuristic operators; it takes profit of parallel implementation techniques [5]. Having also other complementary advantages and drawbacks, machine learning (ML) techniques are also recently considered in such hybridization schemes [6,7]. Note that such hybridization can also provide dual bounds for large instances of optimization problems, using aggregation or decomposition techniques to compute dual bounds for heuristically

Problem Statement and Related Work
This section defines our technician routing problem, situates the closest problems from the literature and presents state of the art methods for these related problems, focusing particularly on mathematical programming and meta-heuristic approaches and their hybridization.

Problem Statement and Notation
The variant of VRPTW considered in this paper, namely TRPTW for Technician Routing Problem with Time Windows, is described by using the terminology from [18]. The routing optimization concerns a set of I vehicles with J requests and denoted jobs or interventions, in different locations of customers. Interventions require the specific skills of technicians and have time window constraints on the arrival time of the technicians. In the following developments, a technician denotes equivalently a vehicle, as if each vehicle is driven by a single technician, and interventions require only one technician with the required qualifications. One may consider equivalently a technician team, which is defined a priori in the optimization, sharing the same vehicle and working together in the interventions, and there is no splitting of a technician team for realizing different interventions at the same time. The resulting skills of the technician team are mathematically equivalent to a skill set and a subset of interventions that can be realized, as if we have a fictive single technician with these resulting skills. For the sake of simplifying the presentation, the following developments consider one technician by vehicle, with skills associated to this technician. Furthermore, technicians are starting and finishing their tour in defined depots, with time-window constraints for leaving and going back to the depot.
It is possible to plan only a subset of interventions; a penalty in the objective function is associated for each request that is not served. The objective function to minimize is the total length of vehicles' routes, adding penalties for the interventions that are not realized. Such objective function may be observed as a weighted sum optimization, considering the bi-objective minimization of routing costs and unserved requests. Here, a weighted sum optimization makes sense with relevant values of penalties. Some interventions may be outsourced and realized by other companies, and penalties correspond to outsourcing costs. Some interventions may be urgent, and they are preferred be realized within the time horizon with high penalty, whereas other interventions can be postponed for another day and have low penalties. Having high and equal penalties for jobs is equivalent to a hierarchical optimization that firstly maximizes the number of jobs that are planned and then minimizing the route lengths. Mathematically, it is feasible for the optimization problem to have empty routes and only penalization. However, finding routes that serve all the requests is a NP-hard decision problem.
The notation for sets, indexes and input parameters is gathered in Table 1. The jobs are indexed with j ∈ J = [ [1, J]]. The intervention of job j has to begin in time window [T min j , T max j ], its duration is L j and its outsourcing cost is P j . For the sake of unifying notation, we indexed technicians with i ∈ I = [[J + 1, J + I]] and defined nodes n ∈ N = I ∪ J = [[1, J + I]] to index the different locations of jobs and depots for departure and arrival of technicians' routes. D n,n and T n,n denote, respectively, the kilometric distance and transportation time from the location of node n to its of node n . Note that distances are calculated with Euclidean rules so that we have triangle inequalities D n,n + D n ,n D n,n and T n,n + T n ,n T n,n . Technician i ∈ I starts from his depot after T start i and must be back at his depot before T end i . J i denotes the subset of jobs that technician i can complete. Table 1. Notation, sets, indexes and input parameters. J i Subset of jobs that technician i ∈ I can complete with skill constraints.
N i = J i ∪ {i + J} Subset of nodes that technician i ∈ I can visit (depot and skill constraints).

D n,n
Distance from the location of node n ∈ N to the location of node n ∈ N .
T n,n Transportation time from the location of node n ∈ N to the location of node n ∈ N .
[T start i ,  [20]. The simplest variant of VRPTW considering that some jobs can be realized by only a subset of vehicles is the Site-dependent VRPTW (SDVRPTW) [21]. The extension of the VRPTW considering several depots is known as Multi-depot VRPTW (MDVRPTW) [22]. There is one major difference between TRPTW and a site dependent MDVRPTW: contrary to VRPTW extensions, it is feasible to serve any subset of the requests for TRPTW by considering outsourcing costs, i.e., penalties for interventions that are not realized. MDVRPTW and SDVRPTW are sub-problems of TRPTW considering the same transportation graph, high and equal outsourcing costs for jobs and upper bounds of any route visiting each customer (which can be calculated easily): If the optimal cost of such TRPTW is greater than the upper bound, the original MDVRPTW or SDVRPTW is not feasible; otherwise, we have an optimal solution visiting each customer. We note that such outsourcing costs or equivalently fictive vehicles with the penalization cost are used by CG approaches for VRPTW problems to ensure feasibility for CG iterations [23].
Skill VRPTW is another variant of the VRPTW inspired by technician routing applications [24,25]. Skills model the ability of vehicles (or technicians) to proceed to specific locations (or intervention with specific qualification). Skills are equivalent to site-dependencies, and skills VRPTW and SDVRPTW have the same set of feasible solutions. Skill VRPTW is an extension of SDVRPTW, and costs of vehicle usage are increasing with skills, whereas costs depend only on traveled distance for SDVRPTW. Note that Heterogeneous Fleet VRPTW (HFVRPTW), where vehicles are differentiated in the travel time and costs, is a common extension for skills VRPTW and SDVRPTW [26]. Indeed, site-dependent constraints can be considered with infinite costs (or a trivial upper bound of any feasible solution) in a HFVRPTW. For the sake of having realistic solutions with real-world application, optimal routes in skill VRPTW can be unbalanced; some works consider load balancing in the objective function [27]. This results in a multi-objective extension and minimizing routing costs conjointly with the maximization of route balance for better equity and robustness relative to uncertainty [28]. Recently, an extension of skill VRPTW was proposed by taking into account that intervention time depends on the skill levels [29].

Related Technician Routing Problems
TRPTW can be observed as a simple case of many technician routing problems. A large diversity of technician routing problems was studied, with many variants of additional constraints [30]. The technician routing and scheduling problem (TRSP) and workforce scheduling and routing problem (WSRP) denote such extensions of TRPTW [31,32].
A rich technician routing problem was proposed for the ROADEF Challenge 2007, organized by the French Operations Research society (ROADEF) [18]. This technician routing problem optimizes the daily assignment of technicians into teams, of teams to tasks and of teams to daily routes in a multi-period problem; site-dependency is implied by the teams composition. There are also precedence constraints among interventions; for instance, some interventions must be preceded by the delivery of some equipment. The TRSP considered by [31,32] considers a single period version of the ROADEF Challenge 2007 by analyzing the impact of team building optimization. Recently, the VeRoLog Solver Challenge (VSC) 2018-2019 proposed a rich technician routing and scheduling problem with truck routing and team building optimization [19]. VSC differs from the ROADEF Challenge 2007 with additional scheduling constraints for technicians. For instance, a single technician may perform at most one tour per day and must not perform tours on more than five consecutive days.
In the TRSP considered by [33], vehicles are equivalent to technicians with no team building optimization, which defines site-dependencies. This TRSP includes additional request requirements for spare parts and tools that can be seen as non-renewable resources. Technicians collect tools and spare parts from a central depot at any time during the execution of their routes. Recently, an extension was provided considering a mix of conventional and electric vehicles, which requires considering recharging constraints [34].
Multi-period and dynamic technician routing problems consider many variants of scheduling constraints. The technician routing problem with experience-based service times addresses the possibility that technicians learn and improve their skills with an impact on performance [35]. Dynamic versions correspond to realistic situations where new requests, including emergency interventions, arise during the intervention tour and requires re-optimization [36,37]. Urgent interventions may induce considering the minimization of service disruption over a time horizon, which is different from standard VRPTW objective function [38]. Vehicle routing problems such as SDVRPTW and MDVRPTW are known to be efficiently solved by optimality using CG techniques and Branch and Price (B&P) or Branch, Cut and Price (B&C&P) [11]. Recent trends solve a large class of vehicle routing problems in a generic way [12]. A crucial point is to solve CG sub-problems as efficiently as possible, with significant works solving the shortest path problems with resource constraints with labelling algorithms [39,40]. However, very few specific works concerned the application of the CG algorithm for the SDVRPTW variant. SDVRPTW was recently studied in unified B&C&P applying for HFVRPTW [26]. Specific works in mathematical programming concerned HFVRP [41]. A unified and efficient approach was provided with the dual ascent method with Lagrangian bounds for HFVRP, SDVRP and MDVRP [42]. A B&C&P algorithm for the multi-depot HFVRPTW was provided by [43].
Metaheuristics are widely used to solve vehicle routing problems [44]. As with exact methods, recent trends develop generic heuristics for a unified formulation of several extensions of VRPTW. A unified formulation of MDVRPTW and SDVRPTW is given by [20]. Such unified extension was solved efficiently firstly with tabu search [20], secondly with Adaptive Large Neighborhood Search (ALNS) [45] and lastly with a hybrid Genetic Algorithm (GA) with adaptive diversity management [46]. Ant Colony Optimization (ACO) is another population metaheuristic that was proven efficient for some vehicle routing problems related to technician routing [47].
Additional constraints in rich VRPTW induce difficulties for metaheuristics, and matheuristics are widely used in such contexts, especially CG matheuristics [14]. The authors of [15] provide a taxonomy of matheuristics for vehicle routing problems: Decomposition approaches are constructive heuristics solving iteratively sub-problems as in [48]; improvement heuristics apply local search heuristics with large MILP neighborhoods as in [49]; and CG-based heuristics derive primal heuristics from Lagrangian relaxations or CG and B&P algorithms as in [13,50]. The hybridization of ML techniques with CG algorithm was proven efficient for some vehicle routing problems [16,17].

Solving Technician Routing Problems
Efficient algorithms for solving technician routing problems are similar to the ones for VRPTW extensions. Since recently, technician routing problems were mostly solved by heuristics. An ALNS algorithm was proposed to solve a TRSP as a SDVRPTW extension [31]. A simple Iterated Local Search (ILS) has been proposed for the same problem and showed excellent results in short solving time [32]. Population meta-heuristics showed also excellent results for technician routing problems. Multi-attribute VRPTW is an extension that is also a basis for technician routing problems; hybrid GAs were shown to be efficient [51]. A Particle Swarm Optimization (PSO) algorithm was proven efficient for solving a TRSP [52]. An ACO heuristic was designed to solve efficiently a multi-period workforce scheduling and routing problem with dependent tasks [53].
Recent progresses in MILP solving allow considering exact methods for technician routing problems with CG algorithms, B&P and B&C&P approaches [54][55][56]. In the previously mentioned TRP with site-dependencies, two approaches considered a matheuristic using an MILP set covering formulations [33,57]. In such cases, an MILP computation selects routes among routes that are calculated by heuristics, similarly with the master problem of the CG algorithm. In such cases, the columns are generated by using only heuristics and without using reduced costs, a parallel ALNS for [33] and an hybridization of ALNS and Variable Neighborhood Descent (VND) for [57].
Due to the large diversity of specific technician routing problems that were studied [30], we focus in the following on algorithms solving VSC 2018-2019 [19] and the ROADEF Challenge 2007 [18]. For such problems, several studies deal with exactly the same optimization problem and public instances; it allows relevant comparisons. The winner of VSC 2018-2019 used a trajectory local search heuristic, which is a hybridization of ALNS and Variable Neighborhood Search (VNS) [58]. A similar quality of results was obtained with a set partitioning matheuristic "Columnwise neighborhood search" [59]. Competitive results were also obtained with a population hyper-heuristic [60] and by using an iterated local search using several neighborhood search operators and destroy/repair heuristics [61].
The winner of the ROADEF Challenge 2007 considered a three phase matheuristic [48]. A preprocessing phase estimates the number of technicians required in each team through the first MILP. Then, a second MILP computes a lower bound allowing job preemption. Multiple solutions are constructed by iteratively solving a third MILP by matching formulation with different parameter values but without introducing randomness. This approach was improved after by [62]. Another two-step decomposition matheuristic provided good results and is ranked fifth [63]. Trajectory heuristics provided also excellent results, with ranks 2, 3 and 4, local search with multiple neighborhoods and operators for [64], ALNS for [65] and Greedy randomized adaptive search procedure (GRASP) for [66].

MILP Compact Formulations
In this section, two alternative MILP formulations are provided for TRPTW, with cuts to improve the quality of the Linear Programming (LP) relaxation for a better efficiency of B&B tree search.

First MILP Formulations
We define binary variables u i,n,n ∈ {0, 1} for each technician i ∈ I and each pair of nodes (n, n ) ∈ N 2 , where u i,n,n = 1 if and only if technician i proceeds to node n immediately after node n. For each technician i, such definition unifies routes between jobs locations u i,j,j for each pair of jobs (j, j ) ∈ J 2 , a route between depot and job location, respectively, u i,i,j = 1 and u i,j,i = 1, is defined to perform intervention j ∈ J as the first or last job of the day and u i,i,i = 1 is the case where technician i has no job in his working day. With such a definition, we need only to consider variables u i,n,n for each technician i ∈ I and each pair of nodes (n, n ) ∈ (J i ∪ {i}) 2 = N 2 i , and the other variables are zeros with skill constraints and definition of departure and arrival depots. Furthermore, we introduce continuous variables t j by representing the start time for job j, bounded with the time-windows constraints.
It results in the following MILP compact formulation below. The objective (1) for minimizing is composed of two linear expressions of variables u i,n,n , the total length of the routes and the penalties for postponed or outsourced jobs j having zero values for variables u i,n,j for all i ∈ I and n ∈ N . Constraints (3) are implied by the definition of variables, u i,j,j = 0 for each technician i ∈ I and each job j ∈ J . Constraints (2) are elementary constraints: Each intervention is realized at most one; it can be postponed or outsourced if variables u i,n,j are zeros for i ∈ I and n ∈ N . With triangle inequality, it is sub-optimal to visit a customer twice or more times. Indeed, having triangle inequalities D n,n + D n ,n D n,n and T n,n + T n ,n T n,n , removing jobs planned twice, is still a feasible solution, and it decreases the objective function; there is an optimal solution allocating exactly one technician per job. Elementarity constraints are necessary here in the MILP model as multiple visits induce a negative penalization and, thus, a high bonus in the objective function. Constraints (4) are flow conservation constraints for modelling routes for technicians: the flow leaving and the flow arriving to a node for technician i is the same flow that is leaving from the node. Constraints (5) initialize flow constraints; there is exactly one arc leaving from a depot, and it can be u i,i,i = 1 in the case of an empty route for technician i so that, in any case, the leaving flow is exactly one. Note that constraints (4) and (5) imply that ∑ n∈N i u i,n,i = 1 for all technician i, which is a one value flow for going back to the depot; such constraints are not necessary. Constraints (2) can be seen as flow capacity constraints. Constraints (6)-(8) express time windows constraints, respectively, between two consecutive jobs and for the starting and finishing time of technicians. With constraints (6)-(8), the formulation does not require any sub-tour constraints. Coefficients M 1 j,j , M 2 i,j and M 3 i,j can be chosen as follows: ∀i ∈ I, ∀i ∈ I, ∀j ∈ J i , ∀i ∈ I, ∀j ∈ J i , ∀i ∈ I, ∀n, n ∈ N i u i,n,n ∈ {0, 1}, Infeasibilities due to time windows can written directly on variables u i,j,j for j, j ∈ J i by following rules allowing the removal of variables that could be non zeros in the LP relaxation. If T min j + T j,j + L j > T max j , then u i,j,j = 0; it is not possible to reach intervention j on time after having realized intervention j. If T start i + T i,j > T max j , it is not possible to reach location j on time by starting from depot i, then u i,i,j = 0. If T min j + T i,j + L j > T end i , then u i,j,i = 0; it is not possible to reach depot i on time if j is realized by technician i. For both previous cases, it is actually unfeasible with triangle inequalities that technician i realizes intervention j; all the variables u i,j,n and u i,n,j can be removed.
If constraints (6)-(8) define integer solutions, such constraints are known to induce LP relaxations of a poor quality [67]. In the following, cuts are provided to improve significantly dual bounds without increasing too much the computation time for the dual bounds. Sub-tour cuts, such as in the Traveling Salesman Problem (TSP), can cut continuous solutions induced by the "big M" constraints. One can forbid firstly the sub-tours between two jobs j and j with constraints u i,j,j + u i,j ,j 1. Such cuts can be lifted, observing that cases u i,j,j = u i,j ,j = 0 are feasible here when the vehicle is not leaving the depot. It induces tighter cuts.
Note that it is possible to have u i,j,i = u i,i,j = 1 when a technician does a single intervention in the route, and previous cuts cannot be extended for j, j ∈ N i .

Four-Index MILP Formulation with Ordinality
Many alternative MILP formulations can be designed with the formulation variants of sub-tours in TSP as analyzed in [68]. Similarly to the formulation of Fox-Gavish-Graves (FGG) [69], variables can be modeled with 4-index variables z i,j,n,k ∈ {0, 1}, such that z i,j,n,k = 1 if and only if the k-th job (k ∈ [[1, K i ]]) that is operated by technician i is job j ∈ J i and then the job (of depot) n ∈ N i is reached after intervention j. Without time window constraints and without K i upper bound on the number of intervention per route, this would induce a number of variables in O(I × J 3 ) instead of O(I × J 2 ) for the last MILP formulation. A first, the issue will be to analyze if the better quality LP relaxation compensates the increasing size of the model. The maximum number of interventions in a day K i allows having O(K × I × J 2 ) variables, with K = max i K i . With the variables of the previous formulation, there is the following relation.
The following MILP formulation models TRPTW with variables z i,j,n,k is as follows.
∀j, j ∈ J , t j + L j + T j,j t j + 1 − ∀i ∈ I, ∀j ∈ J i , ∀i ∈ I, ∀j ∈ J i , Note that the objective function is in three parts, and D j,n z i,j,n,k terms do not count the distance from the depot to the first jobs of technicians. Constraints (13) are straightforwardly equivalent to the previous elementarity constraints (2). Constraints (14) are implied by the definition of variables, z i,j,j,k = 0 for each technician i ∈ I, j ∈ J i and k ∈ [[1, K i ]]. Constraints (15) initialize flow constraints, as in (5); the difference is that flow initialization is written after the first job, and this flow is not anymore equal to one as there are no variables indicating straightforwardly if an empty route is chosen. In such case, all the variables related to the technician are zeros. Constraints (16) are the equivalent of previous flow constraints (4) for modeling routes for technician; the balance of flow is here considered regarding the ordinality k. Constraints (17)- (19) express time windows constraints, respectively, between two consecutive jobs and for starting and finishing time of technicians, similarly with constraints (6)- (8). Constraint (17) is the same constraint set as (6) by using relations (11). Constraint (18) is equivalent with (7), observing that j is the first intervention of technician i if and only if there exist n ∈ N i such that z i,j,n,1 = 1, which is equivalent to ∑ K i k=1 z i,j,i,k = 1. Constraint (19) is equivalent with (8), observing that j is the last intervention of technician i if and only if there exist k such that z i,j,i,k = 1, which is equivalent to ∑ K i k=1 z i,j,i,k = 1. As in the previous MILP formulation, time window constraints may induce impossible transitions. Relations (11) allows considering the same preprocessing rules than in the previous MILP formulation. A question is to analyze whether such formulation avoids sub-tours. Actually, feasible continuous sub-tours may occur, for example, z i,1,2,1 = z i,2,1,2 = z i,1,2,3 = z i,2,i,4 = 0, 5. This configuration avoids paying a distance from 1 to 2, with a sub-tour of size 1. The following cuts can also be added to avoid such situations, similarly to previous subtour cuts u i,j,j + u i,j ,j 1 − u i,i,i .
Note that ∑ 1} is the flow arriving at the arrival depot. It is zero if and only if technician i has an empty route; it is the equivalent of u i,i,i in the first MILP formulation.

Matheuristics Based on Compact MILP Formulations
In this section, the previous work on MILP compact formulations is applied in order to design constructive and local search matheuristics.

MILP Neighborhoods for Local Optimizations
Neighborhoods and heuristic operators can be designed generically from an MILP formulation defining a sub MILP of reduced size (that can be parametric) and imposing values relative to some variables. Such local optimization solving sub-MILP are used in "variable fixing" heuristics (as in [5]) or in Fix-and-Optimize (FO) (as in [70]). In special cases, heuristic operators and local optimization can be implemented without using MILP computations; the computational experiments aim to provide results with respect to designing operators for the TRPTW problem and to select the promising operators that require careful implementation.
Commonly to the two last MILP formulations, different variable fixing strategies are processed on assignments x i,j ∈ {0, 1} for i ∈ I, with x i,j = 1 if and only if job j is realized by technician i. x i,j is binary and a linear expression of previous variables u i,j,n and z i,j,n ,k .
A constraint x i,j = 0 for a given i ∈ I and j ∈ J imposes that variables u i,j,n , u i,n,j or z i,j,n,k are zeros for all k and n ∈ N . A constraint x i,j = 1 for a given i ∈ I and j ∈ J imposes that variables u i ,j,n , u i ,n,j or z i ,j,n,k are zeros for all j , k and i = i with elementarity constraints (2) and (13). Such choices of fixing variables allow having balanced sub-problems as in divide-and-conquer strategies instead of the straightforward fixing on original variables u i,j,n and z i,j,n ,k , which results in unbalanced sub-problems. This is denoted with the fixing operator F assign,= .
A variant is the fixing operator F assign, with inequality constraints. A constraint x i,j 0 is equivalent to x i,j = 0 such as previously described, and a constraint x i,j 1 allows another technician to realize job j.
Previous strategies limit the assignments between jobs and technicians. One may wish to impose some orders in the interventions processed by a given technician. Additional constraints t j t j are little nterest as they are almost always non active in the B&B search because of the weakness of big M constraints. Hence, a fixing strategy is based on the implications on the u i,n,n or z i,n,n ,k variables.
With many constraints F order i,j,j , it breaks possibilities of sub-tours, which has positive impacts for the LP relaxation quality and the B&B solving capacities of such reduced MILP.

Greedy Heuristics Iterating along Technicians
Greedy algorithms for TRPTW can be designed by iterating technician by technician. For a given technician, a local optimization can minimize the distances and the waiting times, avoiding holes that can induce scheduling a few number of jobs. Algorithm 1 implements such strategy iterating along technicians with a generic FO strategy where solutions are encoded as the assignment matrix (x i,j ).
Solve the MILP and update assignment matrix in x end for return assignment matrix in x and the solution cost Variable fixing strategy in the local optimization is a crucial input parameter in Algorithm 1. Applying F assign,= on the current solution for the previous technicians, the assignments are kept, and the optimization for the current technician i focus on the jobs that were not planned previously. In order to not seek to reoptimize a solution that was previously defined, one may completely remove this optimization and consider only technician i and the remaining jobs.
By applying F assign, on the current solution for the previous technicians, the current technician can take previously assigned interventions; such situations occur if it improves the distance cost and it is not possible to assign another job to the current technician that is still unplanned. To avoid reoptimization of the previous technicians, one may consider the F order constraints to keep the previous sequence of interventions of the previous technicians and only reinserting some jobs in the planning of the current technician.
In Algorithm 1, one also defines the order to iterate among the technicians. One may iterate the following a randomized order. To take into account the multi-skilled technicians, we iterate following the increasing number of jobs among technicians; the technician can process with skill constraints so that specialized technicians can process the few jobs they have the qualification for, and multi-skilled technicians at the end of the heuristic focus on the jobs that were not previously assigned. This sort of iteration aims to maximize the chances to realize the maximal number of jobs.

Greedy Heuristics Iterating among Jobs
Another idea to process a greedy algorithm for the TRPTW is to process jobs by jobs. Algorithm 2 implements such FO strategy, where parameters include the order to iterate among jobs and an insertion operator with a maximum number of jobs for the insertion.

Algorithm 2
Greedy FO heuristic iterating along jobs Input: -O an order to iterate among jobs, n a maximal number of jobs for the insertion, -Insert, an insertion operator defined with variable fixing strategies F Initialisation: for all the jobs j in the current planning; • K-insertion without order reoptimization: By fixing current orders, k chosen jobs can be inserted in the current planning MILP with the previous variable fixing procedure and F order i,j,j applying for the jobs and technicians of the current planning; • 1-insertion without order reoptimization: The previous strategy with special case k = 1 can be iterated easily without MILP computations. Indeed, considering the best insertion of a given job in the current planning makes at most I + J possible insertions.
Checking the feasibility of a sequence of job is an earliest date scheduling that applies only for improving costs in theses I + J possibilities. Another parameter concerns the sorting rules to consider jobs. The following strategies can be considered:

•
Sorting the jobs by increasing the number of technicians that can process job j with skill constraints; • Sorting the jobs by increasing value of the summed time windows where the technicians can process job j; • Random sorting of jobs can be used to derive GRASP heuristics. We note that we can hybridize deterministic and random strategies with weighted sums to derive GRASP strategies from the deterministic order.

Local Branching Insertions of Jobs
In the previous FO strategy, the insertion order is fixed a priori. Another alternative is to consider the best insertion among the remaining jobs. Algorithm 3 implements such a strategy.
With n = 1, the insertion operator can try to insert each job in the current planning without modifying the order among the previously planned jobs. Enumerating all the possibilities remains polynomial. MILP optimization can be used at each iteration for any value of n, fixing the assignments out of jobToPlan and considering the following constraints, analogously to Local Branching [71].
Additional fixations for variables with MILP optimization can be the order constraints among the already planned jobs.

Algorithm 3 Best local insertion constructive heuristic
Input: n a maximal number of jobs for the insertion, -F an optional and additional fixing strategy for iter from 1 to J n : insert optimally n jobs of jobToPlan in x remove the inserted jobs from jobToPlan end for return assignment matrix in x and the solution cost

VND Local Search Matheuristic
Once a feasible solution is built with previous constructive matheuristics; this section improves the current solution in a local search procedure.

General Algorithm
We present here the general local search algorithm described in Algorithm 4, similarly to [4]. The neighborhoods are defined here using MILP definitions. MILP neighborhoods allow exploring large neighborhoods, using recent progresses of MILP solvers, having fewer and better local optimums than small neighborhoods approaches. Many MILP neighborhoods can be defined; Algorithm 4 changes systematically the choice of the neighborhood within the local search, similarly with multi neighborhood search approaches. It induces that a local optimum for the entire local search is a local optimum for all the neighborhoods considered in Algorithm 4. It ensures having less and better local optimums than classical local search approaches. The stopping criterion could be a maximal time limit or a maximal number of iterations or being in a local extremum for all neighborhoods.

Algorithm 4 VND with MILP neighborhoods
Input: an initial solutions, a set and order of neighborhoods to explore Initialization: currentSol = initSolution, N = initial neighborhood. while the stopping criterion is not met define the MILP with incumbent currentSol and the neighborhood N ) define currentSol as warmstart Usually, MILP neighborhoods are defined for small sub-problems where exact B&B converges quickly to optimality. MILP neighborhoods are defined with three characteristics, for an efficient MILP neighborhood search includes small time limits without the optimality guarantee. Parameterization is empirical for a good trade-off between solution quality and time spent in MILP solving. The current solution is the primal solution given by the last B&B resolution, and it is also defined as warmstart for the next B&B resolution to improve the efficiency of B&B primal heuristics, enabling RINS or Local Branching heuristics from the beginning. This ensures that the solution given by the MILP resolution is at least as good as the current solution at each iteration. This algorithm is, thus, a steepest descent algorithm.

MILP Neighborhoods
Multiple types of large neighborhoods can be defined with the variable fixing strategies previously defined. We consider three main types of neighborhoods: • N tec i : For a given technician i, the neighborhood reoptimizes the routing of technician i considering any job j, even those that were planned in other routes. In other words,

Sequence of MILP Neighborhoods
A key point in Algorithm 4 is the sequence of neighborhoods. Applying Algorithm 4 with a single type of neighborhoods allows analyzing the impact of neighborhoods in the quality and number of local extrema. In this case, different neighborhoods are compared starting with the same initial solutions and analyzing the different qualities of the local minimum where the steepest descent local search converges. A traditional idea with VND is to increase the size of the neighborhoods when a local minimum is reached for a better compromise between computation time and the local minimum reached. Having many types of neighborhoods, the stopping criterion is firstly to reach a local minimum for all the types of neighborhoods. The choice of neighborhoods is nested, alternating deterministically the order of neighborhoods and stopping the resolution when no neighborhood can improve the current best solution.
Using neighborhoods N tec i or N pair i,i , the computations with the single technicians and the pairs are computed successively in a loop. After such a loop without improvement, all the classical insertions and 2-opt neighborhoods were tried so that the local minimum provided is a local minimum for all these classical neighborhoods. Using neighborhoods N jobs J , many partitions of jobs can be designed, similarly with POPMUSIC (Partial optimization meta-heuristic under special intensification conditions, see [72]). A partition induces successive computations of N jobs J along the partition, and several partitions are computed successively.

Parallel Heuristics
Parallelization is a key issue for the efficiency of these heuristics. The MILP-VND local search can be naturally parallelized, exploring several (many) neighborhoods in parallel. This is particularly interesting in advanced phases of the VND when very few neighborhoods give improvements and in the last iteration in a local minimum for all the neighborhoods. The MILP-VND is a pure intensification strategy of local search, without any diversification. One of the advantage of matheuristic construction of solution is that many strategies can be designed or parametrized, as in [5]. Combining the matheuristic of this section is natural in a multi-start heuristic, starting from one solution given by one of the many constructive matheuristics and improving the solutions with MILP-VND local search, similarly to [73]. In this case, parallelization is more useful for parallelizing the multi-start heuristic in a high level team-work parallelization.

Dantzig-Wolfe Reformulation and CG Matheuristics
We investigate now the extended reformulation of Danzig-Wolfe derived from the previous compact formulation and how to derive heuristics.

Extended Formulation and Column Generation
Similarly to [11], we have an extended MILP formulation by enumerating all possible (and non empty) routes P i for each technician i. The cost of a route p ∈ P i is denoted L i,p . The variables of the MILP formulation are z i,p ∈ {0, 1} such that z i,p = 1 if route p ∈ P i is chosen, and y j ∈ {0, 1} with y j = 1 if job j is not planned. For all j ∈ J and p ∈ P i , we used the characteristics function 1 1 j∈p ∈ {0, 1} defined with 1 1 j∈p = 1 if and only if job j is realized in route p. TRPTW can be written as the following extended MILP.
Contraints (28) express that a job j is either postponed in the case y j = 1 or at least in one of the selected routes. Note that, here, elementarity is not a necessary condition; it is implied by the optimality with triangle inequality and not necessary for the validity of the MILP formulation contrary to previous compact MILP formulations. Constraints (29) enforce that at most one route is chosen for each technician, and it is possible to have only zeros for technicians that are not leaving their depots.
The difficulty is that P i cannot be enumerated. The LP relaxation of this extended MILP formulation is solved by the CG algorithm, adding iteratively new routes. By having a subset of routes C i ⊂ P i , the Restricted Master Problem (RMP) denotes the previous LP relaxation applied to a subset of variables.
π and σ denote, respectively, the dual variables associated with the constraints (28) and (29). The RMP is always feasible as it contains a solution with empty routes for technicians and the maximal number of penalization. Using strong duality, the dual problem of continuous RMP (31) is feasible and has the same value as the primal problem.
∀j ∈ J , π j P j (y j ) In the dual areas, cuts have to be added until there is no violation of constraints (33). The reduced cost associated to a route (i, p) appears; RC i,p = L i,p + σ i − ∑ j∈J 1 1 j∈p π j . Variables must be added in the RMP if their reduced cost is strictly negative. If the remaining variables have a positive reduced cost, the RMP gives the LP relaxation of (27)- (29). This negativity question among a non-enumerable set is formulated as an optimization problem, minimizing the reduced costs of the columns that are not in ∪ i C i . These negativity questions are decomposed for all technicians, RC * = min i RC * i , where RC * i is the minimization problem over the possible routes for technician i, and RC * 0 ensures that the RMP gives the LP relaxation of (27)- (29). The computations of RC * i are independent for the technicians and can be formulated using MILP formulations of Section 4.
LP relaxation of (27)-(29) is calculated by Algorithm 5, with an initial set of columns as input. The initial set of columns can be empty, and RMP is always feasible thanks to the penalty costs. Matheuristics can initialize the CG algorithm adding the columns corresponding to primal solutions.

Input:
C set of initial columns. do: solve RMP (31) with columns defined in C store dual variables σ and π and optimal cost from (31) for each technician i ∈ I : solve (35) to optimality with last (σ, π) values if CR * i < 0 then add the optimal column to C end for while: columns are added in C return the last cost of the RMP (31) Remark 1. Constraints (28) are written with inequalities, and the MILP model is still valid with equalities instead. One could also have equalities in (29), defining explicitly empty routes. Using inequalities induces a better stability of the column generation algorithm, with signed dual variables π, σ [11].

POPMUSIC-CG Decomposition
Dual variables associated with jobs can be interpreted as marginal costs to realize such jobs. When continuous RMP does not use any variable visiting a job j, we have π j = P j . When the dual variables are not stabilized, it induces sub-problems generating columns with jobs with the highest values π j in (35). Hence, independent computations of (35) are likely to visit the same jobs, with the highest values of π j . These generated columns are likely to be redundant in the recombination induced by the next RMP computation. It suggests incorporating diversification in the CG algorithm. We investigate combined subproblems for a subset of technicians I 0 ⊂ I (typically two or three technicians), imposing a total diversification in these sub-problems with the following MILP formulation.
Sup-problem (36) can be seen as a concatenation of the technician sub-problems (35) in I 0 ⊂ I. The first constraint enforces that each job can be affected with respect to at least one technician in I 0 for a total diversity among the new columns. One may consider a multi-objective extension and generate columns for I 0 ⊂ I and consider the partial diversity in another objective.
The process gives rise to Algorithm 6, where the partitions of technicians vary each iterations, similarly to POPMUSIC heuristics (Partial optimization metaheuristic under special intensification conditions [72]). Note than POPMUSIC matheuristics are efficient as primal heuristics without computations of dual bounds for vehicle routing problems [74].

Input:
C set of initial columns. do: solve RMP (31) with columns defined in C store dual variables σ and π and optimal cost from (31) compute P I a partition of I in small subsets for each subset I 0 ∈ P I : solve (36) with a matheuristic with last (σ, π) values for each column c with a negative reduced cost add the column to C end for end for while: columns are added in C return the last cost of the RMP (31) and updated set of columns C

Solving CG Sub-Problems
To solve sub-problems (35) with respect to optimality, note that the CG sub-problem related to i is exactly the same as the VRPTW sub-problems with the jobs J i . Constraint (35) can be formulated as an elementary shortest path problem with resource constraints (ESPPRC). If ESPPRC is a NP-complete problem, it is solved efficiently for reasonable instance sizes for VRP problems with labeling algorithms [23,40]. With several technicians in sub-problems (36), the nature of sub-problem changes and exact labeling algorithms will be even more difficult. In this study, generic MILP solvers are used to solve exactly such sub-problems and to validate the contribution of such a decomposition scheme.
However, the computations to optimality are required only to prove the optimality of the RMP at the last iteration. For the first CG iterations, we only need to generate negative reduced cost solutions. The exact computations are time consuming in Algorithm 5; it can be replaced by heuristics to generate quickly negative reduced cost solutions. The choice of algorithm at each iteration seeks for the good trade-off between computation time and quality (negativity) of the generated columns. This trade-off is dealt using different heuristics with different computation times. For instance, the heuristics of Section 4 apply. In Algorithm 2, specific sorts of jobs can be processed for CG sub-problems: One may consider firstly the jobs with the highest dual values π j , especially with the non-assigned jobs in the RMP with a maximal dual value π j = P j , rather than the jobs with a low dual value. Hence, one may iterate in Algorithm 2 the decreasing values of π j . One may also consider π j − D i,j in order to take into account the remoteness of job j.
Sub-problems (36) can be solved to search columns with a negative reduced cost using matheuristics; the number of technician is restricted, and matheuristics restrict the size of MILP computations, removing jobs similarly to Section 4. The POPMUSIC CG scheme applies only for generating columns with a negative reduced cost. To prove the optimality of the RMP, Algorithm 5 applies for the remaining iterations.
Algorithm 7, decomposes the search of columns with a negative reduced cost by using only computations with one technician. These decomposed sub-problems can also be solved as ESPPRC with [23]. This scheme allows generating the best individual columns as in Algorithm 5 and also good complementary columns for these best individual columns. The solutions of (35) with Algorithm 7 are uneven in the reduced costs following the order of priority of the technicians, whereas a straightforward matheuristic search applied to (35) would have a smoother repartition of the reduced costs. A stake of the computational experiments is to determine whether these two generation schemes are complementary or if one dominates the other. To generate both types of columns, Algorithm 7 can be applied to generate initial columns before a smoothing local search phase in order to optimize the objective function of (36) with a MILP-VND similar to Section 4.

How to Select Partitions of Sub-Problems?
A key point for the developments of Section 5.2 and Algorithm 6 is the choice of the partition of sub-problems. Having a subset of technicians with disjoint competences, the joint sub-problem (36) is decomposable, an optimal solution is obtained to concatenate optimal solutions of single sub-problems (35) and there is no benefit to consider grouped sub-problem (36) instead of the standard ones. This highlights the importance of selecting technicians with the most similar skills. Note that technicians differ with the locations of the depots; this point does not allow having exactly symmetrical sub-problems that can be considered as one sub-problem type with a multiplicity of columns that can be used in the RMP, as in [12].

Formulation as Clustering Problems
A first approach is to define a distance among technicians and a null distance for technicians having exactly the same number of skills so that clustering technicians induces minimizing intra-cluster distances. A Hamming distance can be used in that goal. A first approach would be to consider the Hamming distance of the skill subsets. In order to consider the most similar sub-problems (35) in clusters, it is actually better to use the Hamming distance in the vector of J Booleans, indicating the possibility of a technician to realize a given job. With our notation, this distance d among technicians can be written as follows.
A second point before using the clustering algorithm is to define the number of clusters. Here, the crucial point in the application is to deal with a restricted number of technicians in each sub-problems, which is relatively small because of the difficulty of solving sub-problems (36), typically N = 2 or 3 in our numerical experiments. Thus, we consider P = I N clusters, and we consider a clustering problem such as P-means or P-medoids with a cardinality constraint that each cluster is of size at most N.
Having cluster of sizes N = 2 or 3 brings about the fact that cardinality constrained adaptations of standard local search such as Lloyd's heuristic (also called k-means algorithm) are not adapted frameworks. One can use simple greedy heuristics to create such clusters. In Algorithm 6, it is also useful to diversify the partitions. Several good partitions can be constructed by using a randomized greedy algorithm, for instance, using one of the three best possible insertions in a randomized manner instead of choosing the best local insertion.

MILP Formulation with Two Technician by Cluster
In the case where technicians are grouped by at most two, the clustering optimization can be formulated with an MILP formulation similar to assignment problems. By defining binary variables x i,i ∈ {0, 1} for all i < i x i,i = 1 if and only if technicians i, i are grouped in the same cluster, the clustering optimization can be formulated as following MILP.
∑ i >i Constraints (39) enforce that each technician is associated to at most one cluster so that clusters are of size at most two. Constraint (40) enforces that at most one technician is not related to another one. Note that such situation is unavoidable and happens when the number of technician is odd. In the case of even numbers of technician, one can remove constraint (40) and consider an equality in (39). In this case, the problem is very similar with assignment problems, which ensures that the MILP formulation is efficient for large values of I. Even in the general case, the small values render solving ILP easy. One can also consider an extension with at most three technicians by cluster; the bottleneck is the resolution of Constraint (36) here and not optimal computations of such partitions.
As previously mentioned, it is useful in Algorithm 6 to have diversified partitions. By solving the previous MILP formulation, it is easy to have one optimal solution, and many solutions with similar costs can be generated using local moves with 2-opt exchanges. Indeed, considering only competences, there exist many symmetries with 2-opt exchanges with technicians having the same skill sets.

Dealing with Heterogeneity of Sub-Problems
Sub-problems can be very asymmetric in terms of sizes, dealing with different subset of jobs J i with the skill constraints. At the end of the CG algorithms, the sub-problems dealing with less jobs are more likely to have finished the CG process with fewer columns to generate. Some sub-problems can be skipped a priori for some iterations. More generally, CG sub-problem solving should be more frequent for sub-problems with the higher cardinality of J i . Deterministically, the frequency of solving a sub-problem i can be set to J i /J . Another criterion is to focus on the sub-problems that produced the most negative columns in the last iterations. One can also use a Reinforcement Learning algorithm like ε-greedy approach to update the probability of choosing a sub-problem. For such approach, rewards model that columns with a negative reduced cost are found, and one may also count separately as another reward that columns with a significatively negative reduced cost are found.

Related ML Techniques for Decomposition in Mathematical Programming
We mention here the related ideas using ML techniques for decomposition methods in mathematical programming. Learning techniques can also be useful for a CG algorithm to solve sub-problems using the optimal solutions of previous sub-problems, as in [17]. In this work, it allows speeding up the resolution for sub-problems, which is the most time consuming part in CG algorithms. Such ideas are not implemented in our approach; it is complementary with the ML-guided stabilization proposed in this paper. We note also similarities with [9] in the case of a Bender's decomposition structure rather than a DW decomposition. In [9], standard sub-problems are also grouped, and an appropriate clustering problem is defined and solved with ILP formulations an matheuristics.

Tabu Search Intensification
A general difficulty to implement CG algorithms is that dual variables converge erratically. This is closely related to properties of linear optimization as the optimums are extreme points. A consequence for CG algorithm is that many iterations are necessary to have good dual variables to generate the most appropriate columns, with an erratic convergence of dual variables. Stabilization techniques enforce smoothing the convergence of dual variables in order to reduce the number of iterations that converge for the CG algorithm and, thus, reduce the computation time needed for CG convergence. Mathematical properties can not only be exploited for stabilization but also heuristics [75]. This section is motivated by the following fact: Many good columns can be obtained with slight modifications from a good column. At one iteration, high values of π indicate the jobs that are likely used in a new column. A kernel of jobs with high π j values can be combined in several interesting columns regarding the reduced costs. Using such a property, it is interesting to aggressively generate columns with local search moves from interesting columns. This section introduces a Tabu Search (TS) matheuristic to generate quickly a pool of good solutions. Let x i,j = ∑ n ∈N i u i,j,n , the binaries indicating if technician i ∈ I 0 realizes job j. Having M feasible solutions as previously calculated, we denote withx m i,j the value of these variables. To forbid already generated columns, we add the following "no-good-cuts" constraints similarly to [76], which defines a variable fixing strategy.
In order to search around the last solutionx M i,j , allowing K modifications from the M-th solution, it can also be written as a linear constraint that defines a variable fixing strategy.
These constraints allow generating a TS procedure in Algorithm 8 in order to aggressively generate columns. In the case of the classic CG algorithm, the input for I 0 is a singleton. Adding such constraints in the CG sub-problems can be solved with Branch and Bound for small MILP computations. For larger MILP computations, these constraints respect the structures of the MILP matheuristics in order to quickly find good solutions. We note that the aggressive generation of negatively reduced cost columns with a tabu search was also introduced in [77] for a significative acceleration of CG in practice.

Algorithm 8 Tabu search intensification Input:
-I 0 a subset of technicians -the current value in the RMP of dual variables (σ, π) -a set of initial columns c ∈ ∏ i∈I 0 P i with a negative reduced cost in (35) -an integer M ∈ N, a maximal number of TS iterations -an integer K ∈ N, a maximal number of assignment modifications TSINTENSIFICATION(I 0 , K, M) //Initialization: MILP, a MILP formulation for (35) related to technician i p = c initial columns Taboo list of columns l = {c} an integer n = 0 to denote iterations //Loop to generate columns with negative reduced costs do: Add constraint (43) in MILP with columns of p Add constraint (42) in MILP with columns of p solve MILP remove constraint (43) in MILP with columns of p update p, the optimal columns in the last MIP update l adding the columns of p with a negative reduced cost in l while n < M and ReducedCost(p) < 0 return l // the list of columns to add in the next RMP

Dual Bounds Derived from CG
As mentioned before, when Algorithms 5 or 7 are terminated with the non existence of columns with negative reduced costs, the value of the continuous RMP is a lower bound of the problem, which is the LP relaxation of the problem (27)- (29). There is a theoretical advantage in computing the LP relaxation of the extended formulation (27)- (29) instead of one of the LP relaxation of the compact formulations [78]. However, it is possible to have equality among these LP relaxations in the case of a tight formulation existing for sub-problems, as in [79]. It is a numerical issue to test if an improvement of the LP relaxation is obtained considering the DW extended formulation and if it compensates the higher computation times to process the CG algorithm.
Computing the LP relaxation of the extended formulation requires at least solving each sub-problems to optimality with non strictly negative optimal values. Actually, the conditions to derive the lower bounds are weaker. At each iteration of the CG algorithm, one can compute dual bounds by using an equivalence between DW decomposition and Lagrangian relaxation.
It can be reformulated as a min-max problem.
For optimality, it is necessary that λ j P j . Thus, we have the following.
∑ p∈P i z i,p 1 ∀i ∈ I, λ j P j ∀j ∈ J , z i,p , y j ∈ {0, 1}, λ j 0 Swapping minimization and maximization provides a lower bound.
v max The independent optimization in y induces that y j = 0. Thus, we have the following.
∀i ∈ I, ∀p ∈ P i In other words, we have following families of dual bounds: v l(λ), ∀λ ∈ ∏ j∈J [0, P j ] where l(λ) is a Lagrangian bound defined by the following.
∀i ∈ I, ∀p ∈ P i The reduced cost at each iteration of CG algorithm with dual values σ i , π j appears with λ j = π j . The reduced cost minimization RC * i = σ i + min z i,p ,p∈P i (L i,p − ∑ j∈k π j )z i,p appears, and we have the following.
Relation (45) makes a relation between the Lagrangian bound l(λ) related to the dual signal (π, σ) and the optimal values of the sub-problems RC * i . Using only heuristic resolution of sub-problems to reach stabilized value of the dual signal (π, σ), optimal computations of each sub-problem RC * i allow having a lower bound, which is a Lagrangian bound. We note that having only dual bounds for all RC * i , with some LP relaxations or dual bounds given by a truncated B&B search allows computing the lower bounds of the Lagrangian bound (45), which is possible even for large size instances where optimal computations of sub-problems RC * i are not possible. Note that if the CG algorithm terminates, the dual bounds is ∑ j∈J π j − ∑ i∈I σ i ; the Lagrangian bounds the value of the continuous RMP, and this result was also obtained by duality in (32).

Primal CG Heuristics
Once the LP relaxation (31) is computed by CG algorithm, the question is how to repair LP relaxation into integer solutions. The following strategies can be processed: • RMP_INT_HEUR(): Solving the integer RMP with the generated columns is a first quick heuristic. Such a heuristic was used in [33,46]; • LAG_HEUR(): The continuous assignmentsx i,j = ∑ p 1 1 j∈pzi,p ∈ [0, 1] implied by continuous relaxation of RMPz i,p can be repaired to obtain primal solutions, similarly with Lagrangian heuristics (the equivalence between CG iterations and Lagrangian relaxation is shown in previous section). Variable fixing strategies, fixing assignments with values 0 and/or 1, can be a first strategy. Once variables are fixed, straightforward B&B search or specific heuristics may apply in order to quickly obtain solutions; • LAG_RINS(): RMP_INT_HEUR() and LAG_HEUR() can be combined in LAG_RINS(), a RINS lagrangian heuristic fixing common values in the continuous assignments v i,j and in the assignments of the primal solution given by RMP_INT_HEUR(). With such an approach, the initial solution from RMP_INT_HEUR() is still feasible for the RINS optimization for a better warmstarting of the Branch& Bound search or heuristics setting this initial solution.
We note firstly that these repairing heuristics can be processed at each iteration of the CG scheme, which would induce using such repair heuristics in a multi-start environment provided by the new columns and LP relaxations at each iteration. We note secondly that these heuristics can also be used in a B&P tree.

Management of the CG Pool
In the standard CG Algorithm, one column is generated by each sub-problem resolution. When many iterations are processed, it may be useful to pay attention to the number of columns in the RMP. Having a very large number of columns may cause memory errors or prohibitive computation times to solve RMP problems. In practice, the column pool may be cleaned regularly, for example, by defining a threshold in the number of columns to remove after the next RMP computation. Nonbasic columns, i.e., with a null value in the RMP computation, may be removed. It is possible that the removed columns are re-generated later; the benefit of the removal is in the reduction in column pool and its implications. Several criteria can determine which columns to remove, for example, storing the last used column in a RMP computation and removing the columns that were used the less recently or computing for all the columns the reduced cost with the updated value of dual variables and removing the columns with the highest reduced cost [78]. It is advised that the individual should avoid using this procedure too often and to keep a certain number of nonbasic columns in the master problem after column removal [78].
With the developments of this paper, especially with Algorithm 8, the columns may be generated very aggressively. This induces calling the column removing procedure more often. A specific strategy for the column pool management is proposed here, taking into account that some columns are useful also in the integer RMP computations to define an integer solution, even if the column is not interesting for the continuous RMP computations and, thus, the convergence of the CG algorithm. In our management strategy of the column pool, the following rules are provided. The procedure is enabled only when the number of columns reach a given threshold. The procedure is called every C ∈ {3, 5} iterations of the CG scheme. The columns used in the C computations of the continuous RMP are marked and kept for the following. A computation of the integer RMP is processed after these C RMP iterations in order to mark and keep the best primal solution (and possibly removing jobs that are used at least twice as in the post processing procedure of Section 5.6). The remaining columns are deleted following the updated value of the reduced cost with the last computation of dual variables in order to have a constant and defined number of columns to reiterate CG iterations.

Using Parallel Computing and a Many-Core Environment
Although the CG algorithm follows a sequential structure adding iteratively new columns, several points may be parallelized to speed up the computations in practice. Firstly, the computations of sub-problems are independent for each technician. The I independent computations of sub-problems shall be parallelized. The parallelization to speed up the phases of the sub-problems resolution is the main issue for the speed up, as the sequential part with the LP computations of RMP (31) is less time consuming. In order to optimize the load balancing of the parallel implementation, a LPT (Lowest Processing Time) strategy is natural, using increasing preprocessing times which are predicted with the sizes of the sub-problems with skill constraints. Furthermore, the taboo search intensification proposed in Algorithm 8 allows a more fine-grain parallelization. Lastly, parallelization can be also used when dealing concurrently with several heuristic strategies to generate columns with a negative reduced cost and to perform several of the many variants of heuristics. Having many possibilities of hybridization allows taking profit from a manycore environment.

Computational Results
Numerical experiments were computed with a quadcore computer with Intel(R) Core(TM) i7-4790S CPU 3.20 GHz with 8 threads and running Linux Ubuntu, with 16 Gb of RAM memory.  18]]distinguishes instances generated with the same previous characteristics. For each instance, jobs have the same penalization cost. Penalties are larger than transportation costs so that the optimal solutions serve a maximal number of jobs. For only instance gotic_15_20_40_ex7.txt, it was impossible to realize all the jobs; the other optimal solutions visit every customer. In the instances, localisation of depots and customers is given with 2D coordinates, and the distance from localisation coordinates (x, y) and (x , y ) shall be computed using the following rule (different from [19]).
The (symmetric) commuting time between locations (x, y) and (x , y ) is given by the following: where speed is a single speed parameter defined in the instances.

Analysis Methodology
Instances were generated with relations among instances to analyze the trade-off between optimality and feasibility for the different optimization algorithms. Instances with d ∈ [[11, 18]] are generated from instances with d − 10 ∈ [ [1,8]] with a higher speed parameter. Hence, each feasible solution for an instance with d ∈ [ [1,8]] is feasible for the corresponding instance with parameter 10 + d, For each group of instances with fixed values of parameters a and b defined, a maximal number of competences are defined with c max . The other instances with c < c max are defined from the one with c max : skills in the instance s are grouped modulo c; if skill q is required in instance c max , the corresponding skill in instance c is q mod c. Such reduction keeps the feasible solutions from the previous model and adds potentially new solutions. Hence, instances with c < c max are less constrained than the original ones with c = c max and have a better optimal value. In the case c = 1, there are no skill constraints, and the model is similar with a MDVRPTW. For computational analyses, we investigate the difficulty of instances regarding speed and competence parameters. This defines four datasets:  18]]) and no skill constraint (c = 1).
Hence, datasets s1 and S1 are MDVRPTW instances, with a correspondence and more constrained instances due to speed with dataset s1. Dataset S1 is less constrained than the three other datasets. Original dataset sc max is the more constrained dataset, comparable to the three other datasets. Comparing the difficulty induced for instances sc max and S1 allows analyzing the difficulty induced by speed and rare competence separately.
In order to compare optimization algorithms in a given dataset, the following indicators are provided:

Implementation and Characteristics of MILP Solving
The numerical experiments reported here used Cplex version 12.6 to solve MILP and LP. The analysis of B&B characteristics was provided by using the OPL modeling language to call Cplex. The matheuristics of Section 4 were implemented with OPL script to iterate successive MILP computations. A first CG scheme was implemented with OPL script also by using the populate mode to observe the impact of keeping several columns with negative reduced cost per sub-problem resolution. The final CG scheme and heuristics for subproblems were coded in Python by using the PuLP library to call for MILP and LP. Using PuLP allows using Cplex or open source tools (such COIN-OR) the same programming interface. However, using PuLP resulted in important computation time to load MILP models, similarly to https://python-mip.readthedocs.io/en/latest/bench.html (accessed on 25 October 2021 ). PuLP is not a relevant choice to solve iteratively many MILP models in short computation times, and the loading time of the MILP model may be much longer. Hence, the computation times were not representative with such implementations. The results of the study conclude the types of heuristics that can be used, which include hybridized and parallelized heuristics for future and better implementations.

LocalSolver
As a benchmark local search solver for our study, we used LocalSolver version 7.5. We note that the results are better now by using the last versions of LocalSolver; we provide here the results with version 7.5, and our instances were given to the LocalSolver Company for their internal benchmark. For the latest versions of LocalSolver (at least with version 10), BKSs are found for most of the instances within 600 s solving time. Note that our instances may be interesting for benchmarking general local search solvers, especially to test the accuracy of local search facing highly constrained instances where some multi-start strategies may be useful. Benchmark codes with LocalSolver are also given in https://github.com/ndupin/SiteDepTRPTW (accessed on 25 October 2021). For the efficiency of the LocalSolver model, we used the definition of variables using LocalSolver lists, and TW constraints are modeled as an objective function accumulating TW infeasibilities, with a hierarchical optimization considering firstly feasibility and then the costs of the routes. Such a choice is similar with LocalSolver's VRPTW example and was validated by LocalSolver experts. Two variants were implemented. The general case considers only the disjoint constraint among Localsolver's lists: the routes of each technician visit disjoints sets of jobs. For the special case where a solution exits without outsourcing costs, one can model it with a partition by adding a cardinality constraint with LocalSolver. The reported results show the results with the general variant; failures are cases where a non-maximum number of jobs are assigned to technicians.

Solving MILP Compact Formulations
Analyzing the characteristics of the B&B straightforward resolution for the different variants of MILP formulations of Section 3 is a key issue for determining how to choose the MILP formulation for the next developments. Table 2 compares the quality of the dual bounds for the MILP formulations (1)-(8), (12)- (19) and (27)-(29) on instances solvable to optimality with a straighforward B&B resolution. For the compact formulations (1)- (8) and (12)- (19), the bounds are compared with and without the cuts (10) or (21) for the LP relaxation and the dual bound at the root node of the B&B tree before branching. For the extended formulation, only the LP relaxation is computed with the CG algorithm. This table highlights that the quality of the linear relaxation of the compact formulation is very weak, with a large gap relative to the Lagrangian bounds. The cuts of Cplex provide a significant improvement of the LP relaxation in all the cases. The cuts (10) or (21) provide also significant improvements that are not redundant with the cuts for Cplex. The FGGbased formulation with the cuts (21) gives significantly the best dual bounds among the compact formulations. Table 2. Comparison of the average gaps to dual bounds given by some variants of MILP formulations to the optimal solutions on a set of 55 instances solvable to optimality with Cplex. The results are gathered on instances with the same number of technicians, skills and interventions. LP denotes the linear relaxation, whereas nod1 denotes the bounds obtained after the cutting plane passes of Cplex at the root node before branching.

MILP Formulations Results
By choosing the MILP formulations and the cuts to add for the integer resolution, these results on the dual bounds shall be compared with the computation time. Indeed, adding cuts (10) or (21) for choosing the FGG-based formulation increases the size of the MIP, which induces more computation time for computing LP relaxations. Let us analyze the best compromise in the improvement of LP relaxation relative to the time spent to compute the dual bound. Table 3 compares the computation time for solving small problems to optimality with the previous variants of compact MILP formulations. The FGG-based formulation has significantly worse results than formulation (1)-(8) with the cuts (10). Adding the cuts (10) or (21) is justified also in the computations to optimality by Table 3.
An explanation of the huge gap between the lightest and the FGG-based MILP formulation comes with the efficiency of the generic primal heuristics. For the different compact formulations, the generic primal heuristics provided bad results to find first solutions, and convergence of primal bound is very long. It is indeed difficult to find solutions with the maximal number of jobs allocated in the first stages of the resolution even for medium size instances. Giving a good initial solution as warmstart to solvers allowed improving significantly the computation time for solving instances to optimality, cutting off earlier a significant number of nodes in the B&B tree. These facts are related to the poor quality of the LP relaxation. The empirical results showed that the efficiencies of primal heuristics were significantly worse with the FGG-based MILP formulations, which explains the difference observed in Table 3.

Solving CG Sub-Problems
If CG sub-problems with one technician having a similar polyhedral structure with the entire compact MILP formulation, the efficiency of straightforward MILP solving is very different. For some mid size instances, it can be less time consuming to solve the entire MILP with a compact formulation rather than solving a single technician CG sub-problem. It is paradoxical, and DW reformulation decomposes the original MILP with iterations of smaller problems. The reason for such paradoxical situations is that sub-problems are of a different nature. Seeing CG sub-problems as shortest path problems, many distances are negative in providing a negatively reduced cost. In the original MILP problems, all the distances are positive. For sub-problems, cycles with negatively reduced cost may appear, and this tends to produce many sub-tours with negatively reduced costs in continuous relaxations. When solving sub-problems with many sub-tours of negative costs with an MILP solver, the lower bounds are of a poor quality and many sub-tours are generated in many nodes of the B&B tree. This induces that the convergence to optimality is difficult in this case. Table 3. Comparison of computation time (in seconds) of both compact MILP formulations of Section 3, with or without the 2-sub-tour cuts.
(1)- (8) (1)-(8) Note that this critical situation appears more often in the first iterations of the CG algorithm when the dual values are not stabilized and have values of π j ; in the worst case π j = P j occurs when penalization is used in the RMP in the case of unbalanced dual variables. In the last iterations of the CG algorithm, having more stable dual variables induces less complications, and straightforward MILP solving is more efficient at these iterations compared to the first iterations. A good point is that primal heuristics of MILP solvers (relying on branching, in particular, diving heuristic) remains efficient for finding columns with a negative reduced cost for sub-problems in the first iterations, which is enough to continue the CG algorithm without solving optimality the first sub-problems. Note that for a given CG iteration, the computation efficiency of sub-problems i significantly varies following the number of jobs that the technician considered in a sub-problem can realize. At the end of the CG algorithm, it tends to solve quickly the sub-problems with the less jobs, proving that no column of negative reduced cost exist, and the sub-problems with remaining columns to add are the sub-problems with the maximal number of jobs with the competence constraint.
With the different nature of MILP optimization searching for columns with negative reduced costs, the differences between the FGG formulation and the most compact MILP formulation are emphasized, with larger gaps in the quality of LP relaxations. In order to provide dual bounds of sub-problems, as in Section 5.7, it is better to use the FGG formulation. For using the primal heuristics and solving in truncated time in matheuristics, the most compact formulation (1)-(8) is more efficient. Note that the efficiency of the sub tour cuts is lessened (10). Table 2 shows that the quality of the lower bounds of the extended DW reformulation outperforms compact formulations at the LP relaxation and also at the root noot of the B&B tree considering the generic cuts implemented in Cplex. Note that this was not the case in [79], and these results are coherent with the known results for the class of vehicle routing problems [11].

Results with Extended DW Reformulation
For many instances, the LP relaxation of the extended formulation has no gap with the integer solution, and the integer RMP already provides an optimal integer solution. For such instances, no branching is needed for a B&P enumeration scheme. Note that such results occurs more often for the highly constrained instances. For lowly constrained instances, continuous combinations of columns often exist for decreasing the RMP objective function. Skill constraints induce some decoupling of technician routes and are helpful for CG approaches. Combined with TW constraints and a limited number of jobs per route, it is in favor of CG approaches with the generation of few and diversified columns. Table 4 compares the efficiency of the constructive matheuristics over all instances. These results can be compared to the other approaches over the same set of instances in Table 5; it illustrates that these heuristics are not efficient in terms of accuracy in finding feasible solutions and in terms of the quality of the feasible solutions found. Generally, reoptimizing previous decisions involves F assign, instead of F assign,= , such as standard greedy approaches. Contrary to [5], considering in parallel these different constructive matheuristics and taking the best solution found induces significant over costs for most of the instances, and it is not enough to find a feasible solution for all the instances.

Constructive Matheuristics Iterating Technician by Technician
Constructive heuristics of Section 4.2 iterating technician by technician are particularly inefficient, and many infeasibilities and mostly solutions of a poor quality result. The local optimization tends to produce dense plannings to assign a maximal number of jobs in the first iterations and last iterations with few jobs to assign. This shape of solution explains the poor results of these heuristics. Firstly, the optimal solution may be very different from this shape of solutions with unbalanced solutions, especially for the easiest instances. Even the reoptimization with F assign, , it is not efficient to balance the routes of technicians; a reason for this is that the first routes have also been optimized with distances. Thus, removing jobs induces few impacts. Such variable fixing does not allow placing another job in the previously computed routes. Secondly, the local optimizations do not take into account that some jobs have to be assigned in the planning for rarity reasons; it explains the poor accuracy of these heuristics. These results conclude that the construction of solution "technician by technician" shall be avoided in a specific implementation.

Constructive Matheuristics Iterating by Buckets of Jobs
The constructive heuristics of Section 4.3 iterating with buckets of jobs are more efficient. The results of Table 4 illustrate that slight improvements are obtained using F assign, rather than F assign,= and that considering the TW amplitude to possibly realize jobs is better than considering only the number of technicians able to realize the job. The reason is that rarity of jobs is a key element in the order of local optimizations. Increasing the number of jobs in the local optimizations tends to increase the quality of solutions, and some counter-examples are found where the 10 jobs optimization results in worse results than using the five jobs optimization. Indeed, even if the 10 jobs optimization contains twice the amount as five jobs optimization and would be seen as better local optimization, the optimality of local optimization in greedy algorithms is not a guarantee of the quality of the final solution. Local Branching insertions of Section 4.4 provides aggressive heuristics that are able to find much better solutions but also frequently results in solutions that do not allocate the maximal number of jobs, especially for the most constrained instances. Table 4. Comparison of the statistical indicators presented in Section 6.1.2 using several constructive matheuristics for all the instances. njob denotes Algorithm 2 grouping in subsets of n jobs, and 1tic denotes Algorithm 1. For job insertion, the initial sorting of jobs is based on the number of technicians that can realize the jobs (nbTec) or with the cumulated time considering TW constraints (timeJob). Variable fixing in MILP local optimization includes F assign, or F assign, .

CG Heuristics
As previously mentioned, the LP relaxation of the extended formulation has no gap with the integer solution for many instances, especially highly constrained instances. The integer RMP already gives an optimal integer solution; thus, the post-processing heuristics furnish an optimal solution. Table 5 shows the average results for all the larger instances, and the partial results on the subsets for the comparative analyses are provided in Table 6. In the cases where the integer RMP is not optimal, the solution quality is still excellent and much better than the constructive matheuristics based on the compact formulations. Significant improvements are provided, removing jobs that are realized twice or more in the integer RMP solution, such situation appeared multiple times. Significant and quick improvements are then provided with RINS diving using compact formulation, providing high quality solutions with a majority of BKS.
The difficulty with CG heuristics is the time required to reach the convergence of the CG algorithm. Compact matheuristics are much faster. Speeding up the convergence of the CG exact algorithm is a crucial issue, and analyses are reported in a following section. Using only heuristics and matheuristics to solve sub-problems allows having quick heuristics. However, it induces a significant decrease in solution quality for continuous and integer RMP. Using only heuristics without matheuristic operators for solving subproblems almost keep some jobs not realized even in the continuous RMP. After such quick phase with heuristic solving of sub-problems, using matheuristics (including the matheuristic stabilization operators) allows accurately finding feasible solutions for all of the instances. These two phases are worthwhile for reaching the quickest computation time such as the stages of the CG convergence compared to the use of exact resolutions of subproblems. To derive the CG heuristics, some few iterations of the CG algorithm are required with exact solutions sub-problems with the stabilization operators. Solving sub-problems with LocalSolver in a short time limit and applying TS matheuristic intensification with Algorithm 8 was also very efficient in that context.

Results with MILP-VND Local Search
Tables 5 and 6 provide the comparative results among local search approaches. For the MILP-VND matheuristics, it allows comparing the quality of the local minimum obtained by using different initialization heuristics. The comparative analyses on subsets of instances in Table 6 illustrate the impact of the difficulty of constraints for the efficiency of local search approaches.

Stabilization of CG Algorithm
If the extended formulation and the CG algorithm allows providing lower and upper bounds of an excellent quality, the critical issue is to make the CG algorithm converge quickly, which is the bottleneck of this method. As previously outlined, using quick heuristics in the first phase and then matheuristics in the second phase with the TS intensification as matheuristic operator allows quicker iteration than using an exact method, and it is the quickest way to reach advances phases of the CG algorithm, but it is not enough to reach the convergence of the CG algorithm. Note that this is particularly interesting for avoiding the exact methods at the first iterations with respect to solving the most difficult sub-problems, as described in Section 6.2.1.
This section provides the results applied to decreasing the number of iterations of the CG algorithm in terms of converging faster in order to validate the different CG schemes. Two experiences with exact methods are provided in this section, and the validation of CG schemes has also an impact for designing the first iterations of the CG algorithm with heuristic and matheuristic resolution of CG sub-problems. Table 6. Comparison of the solution quality of local search approaches for the 17 mid-size instances of dataset sc max , Sc max and S1. LocalSolver results are reported for defined time limits, CG heuristics and VND matheuristics at termination. For VND matheuristics, different initialization are provided, and the multi-start VND considers the reported initialization for a parallel or multi-start implementation. and reoptimization of buckets of 10 jobs N jobs in different partitions allowed considering larger neighborhoods than the traditional ones, contains the other main neighborhood and allowed quick resolution times. As in [4], a very short number of iterations of the MILP-VND using iteratively all the considered neighborhoods were enough for converging to a local minimum; for most instances, 2/3 of such iterations were needed, and the worst cases used five such iterations. MILP-VND allows having quick convergence and scaled better for all the larger instances than exact methods. Contrary to [4], such large neighborhoods are not sufficient to converge to local minimums of excellent quality. The quality of the local minimum varies significantly following the different initialization. Considering the four matheuristic initialization and the null initialization, followed by the MILP-VND provided a multi-start matheuristic that is accurate for finding good solutions, the most difficult dataset sc max provided the worst solutions. Dataset sc max is generally the more difficult one for local search approaches, with less feasible neighbor moves because of the hard constraints. This was not the case of the CG-based heuristics. Using MILP-VND after the LAG_RINS heuristic of an excellent quality allowed polishing the results and increases the number of BKS found, illustrating the advantages, drawbacks and complementarity of local search approaches and CG heuristics depending on the structures of the instances.

Impact of CG Initialization and Multi-Column Generation per Sub-Problem
In this section, we analyze the results considering the standard Algorithm 5 with an exact resolution of sub-problems and generating one or several columns per sub-problem. The second variant is implemented by using the populate mode of Cplex with OPL; we note that it is possible to use all the columns with a negative reduced cost found using callbacks, which is not available with OPL but is now available in some MILP modeling interfaces and libraries such as PythonMIP or Julia JUMP. The generation of several columns is used at most for the five first iterations, and the remaining iterations are processed with the standard Algorithm 5 with at most one column generated per sub-problem. The maximal size of the population of Cplex is set to ten: at most, 10 new columns are added in the RMP. Figure 1 illustrates the convergence of integer and continuous RMP in both variants; the results for some small instances are given with the number of iterations to converge in Table 7. Figure 1 illustrates that CG is naturally stable for TRPTW, with a few degenerated iterations. Degenerescence occurs more frequently for VRPTW [11]. A reason is that the I sub-problems are highly heterogeneous with difference competence sets and location of depots; a degenerated iteration would induce the generation of I degenerate columns, whereas in VRPTW, only one sub-problem is generated when taking into account the symmetry among vehicles, and the probability of degenerate iteration is larger [11]. Figure 1 shows that integer and continuous RMP are very close; it happens for highly constrained instances and also with generation of few columns. It explains the high quality of the RMP_INT_HEUR heuristic. Using several solution found by the MILP solver induces many benefits for the CG convergence: The difference in the quality of continuous and integer RMP is significant, and the dual variables have a less erratic convergence. This is a virtuous circle: Having better RMP values induce better dual values that generate more relevant columns. This explains how the experiments with only five iterations of aggressive CG produce such difference. The first iterations are the most unstable. Table 7 considers, furthermore, an initialization with the optimal integer solution. In many cases, the convergence value of the continuous RMP is obtained here, and there are only degenerated iterations with the optimal iterations. In such cases, the dual initial values are very good, but it is not sufficient to guarantee the fast convergence of the CG Algorithm 5. If the optimal initialization makes the number of iterations decrease on average, there exists some instances where there are more degenerate iterations with the optimal initialization than iterations starting from no column. Table 7. Number of iterations for the convergence of CG algorithm to solve the LP relaxation of the extended formulation compared when only one column is generated per sub-problem, improvements with multi-column generation and providing an optimum of the MILP as the initial solution.
the reduced costs) with a VND local search matheuristic, generating all the columns with a negatively reduced cost obtained from these two phases.
• CG5: This corresponds to the strategy CG4 with TS intensification activated with K = 3 and M = 5.
The dominant computational results are illustrated in Figure 2. The known hierarchies among CG schemes are significant: CG1 is significantly dominated by the other schemes, POPMUSIC and TS intensification induce a very significant acceleration in terms of CG convergence. Comparing CG3 and CG4 allows concluding that the optimization of the total reduced cost with diversification in sub-problems (36) induces a better CG convergence than the CG with the hierarchical decomposition of Algorithm 7. Combining POPMUSIC and TS intensification dominates the single stabilization strategies: POPMUSIC column generation is mostly useful for the first iterations in terms of recombining columns in the next RMP, where it can be difficult to find a solution allocating all the jobs (and, thus, paying no prohibitive penalty). TS intensification is crucial when dual variables have a relative stability, which is helpful for the latest phases of the CG convergence. When the dual variables are changing slightly, it is likely that there exist many columns of very good quality among neighbors of the optimal solution of the previous sub-problem. Aggressively using TS intensification prevents generating neighbors from the previous iterations in further iterations, which explains the gain when using TS intensification. To compute the extended LP relaxation to optimality, TS strategies are prominent for reducing the number of CG iterations. TS is particularly efficient for the final iterations of CG algorithm, but it is also useful in the first iterations with respect to having a deterministic scheme to aggressively generate solutions in the first iterations, as in the first experience of stabilization.

Conclusions and Perspectives
This work provided several results and offers different perspectives. The CG algorithm is efficient for computing dual and primal bounds for the TRPTW coherently with

Conclusions and Perspectives
This work provided several results and offers different perspectives. The CG algorithm is efficient for computing dual and primal bounds for the TRPTW coherently with the results obtained for VRPTW problems. The LP relaxation of the DW extended formulation efficiently guides us toward primal solutions of an excellent quality.Matheuristics based on the compact ILP formulation may be inefficient even with combinations of large neighborhoods in a VND local search, contrary to [4], or with parallelization using a portfolio of matheuristics, contrary to [5]. We note also a similarity with [25]: improving the exact methods and MILP formulations is a first step before applying MILP models to matheuristics, and the first step is crucial for the efficiency of the resulting matheuristic.
If CG heuristics are particularly efficient in terms of quality, the key point is to speed up and stabilize the convergence of the CG algorithm. Two matheuristic components are proposed and proven efficient for TRPTW. On the one hand, the heterogeneity of subproblems allows a POPMUSIC decomposition to generate columns by avoiding redundancy in the served requests for a better re-combination in the RMP computation. ML techniques are used to generate appropriate partitions, illustrating the impact of hybridizing ML and ILP techniques. On the other hand, a Tabu matheuristic intensification for sub-problems avoids generating similar columns in neighbor iterations of the CG algorithm, decreasing the number of iterations to converge. POPMUSIC matheuristic is efficient for the first CG iterations where it is difficult to serve all the requests, whereas TS intensification remains efficient for the last phases of the CG algorithm where dual values are quite stable. For both schemes, heuristic solving schemes are designed to speed up the search of negative reduced cost columns, and light heuristics are firstly processed and then heavier matheuristics are used speed up the time of first CG iterations. A perspective would be to improve the exact solutions of such sub-problems for the validated CG scheme by using labelling algorithms. Note that TS intensification and many matheuristics are compatible with labelling algorithms; a perspective is to implement such matheuristics without using an MILP solver. Possessing many strategies and hybridizations for generating columns is a good characteristic with parallelization in a multi-core environment.
We note that our final approach investigates and combines the three types of matheuristics for VRPTW problems in the taxonomy of [15]: CG-based heuristics as a general framework, with computation of primal and dual bounds; improvement heuristics with local search procedures to improve the integer RMP heuristics and to solve CG sub-problems); and decomposition approaches to initialize CG and also in the sub-problem resolution of CG. In such hybridization, mathematical programming is hybridized with the trajectory of local search metaheuristics. A perspective is also to investigate the collaboration between CG heuristics and population metaheuristics. The aggressive CG scheme and benefits or diversification highlighted in this paper offer perspectives to solve CG sub-problems and generate several columns with GA, PSO or ACO heuristics. Recent advances for these population meta-heuristics are promising for such applications [80][81][82][83].
Our instances are now public, and the graduations of difficulty and highly constrained instances with fewer possible neighbor moves are useful for challenging local search implementation. The first results with LocalSolver were already good, with some difficulties on highly constrained instances; we note that our instances are now in their benchmark instances to help the parametrization of their generic solver. This is a perspective for other generic codes with local search applied for vehicle routing problems.
Other perspectives are to extend these results with extensions of VRP and technician routing problems with more constraints (as in [18,19,34]) or considering multi-objective or robust optimization extensions, as in [84][85][86]. The matheuristic stabilization of CG algorithm may also be extended for other problems where CG algorithm is efficient, even without having heterogeneous sub-problems. It raises the question whether the matheuristic stabilization techniques proposed in this paper would also be efficient in the case of symmetrical sub-problems.