{'fraction_non_alphanumeric': 0.026831080978598838, 'fraction_numerical': 0.008690904127041908, 'mean_word_length': 4.749716578006454, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Exact and hyper‐heuristic solutions for the distribution‐installation problem from the VeRoLog 2019 challenge

This work tackles a rich vehicle routing problem (VRP) problem integrating a capacitated vehicle routing problem with time windows (CVRPTW), and a service technician routing and scheduling problem (STRSP) for delivering various equipment based on customers' requests, and the subsequent installation by a number of technicians. The main objective is to reduce the overall costs of hired resources, and the total transportation costs of trucks/technicians. The problem was the topic of the fourth edition of the VeRoLog Solver Challenge in cooperation with the ORTEC company. Our contribution to research is the development of a mathematical model for this problem and a novel hyper‐heuristic algorithm to solve the problem based on a population of solutions. Experimental results on two datasets of small and real‐world size revealed the success of the hyper‐heuristic approach in finding optimal solutions in a shorter computational time, when compared to our exact model. The results of the large size dataset were also compared to the results of the eight finalists in the competition and were found to be competitive, proving the potential of our developed hyper‐heuristic framework.


INTRODUCTION
VeRoLog, the Euro Working Group on Vehicle Routing and Logistics optimization, has been organizing challenges for the routing community, where each challenge aims to promote the design and development of an applicable and effective algorithm for a particular vehicle routing problem (VRP). The organization of the challenges is a collaborative effort with companies, such as PTV, the leading German company in developing software solutions and consultations in traffic, transportation, and logistics, who organized the first and second editions in 2014 and 2015. The third edition in 2017 was organized by ORTEC, one of the largest providers of advanced planning and optimization solutions and services. They provided a real-world VRP problem involving the pickup and delivery of tools to measure milk quality at a number of farms, for a cattle improvement company.
This paper addresses the problem introduced in the VeRoLog Solver Challenge 2019, 1 which was organized by VeRoLog in cooperation with ORTEC. The challenge presented a new and exciting variant of a vehicle routing problem (VRP), based on a real-world problem of one of ORTEC's clients. This problem is about the delivery of equipment, such as vending machines, to satisfy customers' requests, and the scheduling of a number of technicians each with a certain set of skills, who are required to install the equipment at least a day after the delivery. More formally, the overall problem is an integrated version of the capacitated vehicle routing problem with time windows (CVRPTW) [41], and the service technician routing and scheduling problem (STRSP) [20]. The problem is highly complex, combining two interacting and co-dependent NP-hard routing problems into a single model, each problem having its own set of constraints, making it a unique and challenging topic for VRP researchers and practitioners. The key issue in this challenge is to find efficient and low cost routes for both the delivery and installation, while scheduling the appropriate technicians according to their skills and working days for a given instance. The main objective is related to the total costs, aiming to reduce the total number of dispatched trucks and hired technicians on a single day, and within the planning horizon, in addition to minimizing the penalties incurred due to equipment awaiting installation. Many companies operating in the delivery and equipment installation businesses would benefit significantly from an effective and efficient solution to the proposed problem.
For such real-world complex problems, exact approaches, such as, mathematical programming often fail to provide solutions as the size of the instances get larger, and so heuristic-based search methods are preferred. As we have investigated different solution methods for the integrated problem from the VeRoLog Solver Challenge, we had a similar observation. Hence, in this study, we provide a mathematical model for the problem and report the results on a number of small instances. Additionally, we also present a population-based hyper-heuristic approach for the large scale problem instances. The proposed solution methods indeed provide competitive results with respect to the highest ranked scores for each challenge instance.
The remainder of the paper is organized as follows: Section 2 reviews the previous literature on the VRP, focusing on the CVRPTW and the technician routing problem. Section 3 provides a description of the tackled problem and Section 4 formulates the mathematical model, defining the objectives and the problem constraints. Section 5 describes the applied hyper-heuristic framework. Finally Sections 6 and 7 present the results and conclusion, respectively.

RELATED WORK
The problem we deal with from the VeRoLog Solver Challenge 2019 is a unique vehicle routing problem (VRP), that combines a pickup and delivery problem with time windows with the scheduling of technicians to install the delivered equipment. However, it has much in common with some well-known variants of the VRP, thus in the following subsection we provide a brief overview of relevant previous studies. Following this, we cover selection hyper-heuristics in relation to this study.

An overview of history and variants of VRP
The VRP (vehicle routing problem) is a generic term for a class of combinatorial problems that are concerned with the design of efficient routes for a number of vehicles serving a set of customers' requests, originating and ending at a depot location. A solution to a VRP instance is a tour for every vehicle, such that all customers' locations are visited, and each vehicle finishes its tour at a depot. An optimal solution is the one in which the total distance of all tours is minimized along with the associated costs. The term VRP was first introduced in [23] as a truck dispatching problem, where they modeled the problem of how a homogeneous fleet of vehicles can serve the oil demand of a number of gas stations from a central hub, with a minimum traveling distance. This method was then generalized in [19] to a linear optimization problem: how to serve a number of customers located around a central depot, using a fleet of vehicles with varying capacities. This has been known since then as the VRP problem which is one of the most researched topics in the field of operations research and logistics. Over the past decades, the VRP has grown more popular in the literature, and other variations have been developed to model real-life scenarios.
The capacitated VRP (CVRP), which imposes a capacity constraint on the size of vehicles, is considered one of the elementary variants of VRP from which other variants originated. Among exact mathematical approaches to the CVRP, branch-and-cut [7,45] and algorithms based on the set partitioning formulation [9,30] are the most popular approaches. Many studies also applied heuristic methods [29] and genetic algorithms (GAs) have also been widely used, usually by combining them with local search techniques [46,49,50,55].
A generalization of the CVRP is the CVRP with Time Windows, which imposes a time interval ("time window") on the delivery of each customer's request. Exact methods for the CVRPTW have been successful for cases with up to 100 customers [39], and as a result heuristic and metaheuristic methods have been preferred for solving instances of larger scale. Examples of heuristic methods applied to the CVRPTW problem can be found in [54,57,59,60], and other studies that applied metaheuristic methods such as genetic algorithms, ant colony optimization, tabu search, and simulated annealing are in [10,18,24,62].
In our proposed problem, customers' deliveries are scheduled within a planning period, with each customer requiring one or more visits during this period, similar to the multiperiod VRP problem (MPVRP). In contrast to the MPVRP, where service days are known and the frequency of customers visits is predetermined, the visits in our case are scheduled within predefined time windows. The MPVRP is a well-studied variant in the literature of VRP. In a paper by Rahimi-Vahed et al. [56], a path relinking algorithm is applied to the multidepot periodic vehicle routing problem. This is done by generating a reference set of elite solutions, and combining characteristics from those solutions to find even better solutions. The computational results show that this method produces good results in both run-time and solution quality. Archetti et al. [6] present three ways to formulate the multiperiod vehicle routing problem with time windows, then they solve the problem using a branch-and-cut algorithm. The algorithm was able to find good solutions for small problems with 10 customer orders, but was unable to find many good solutions for larger problems. Alonso et al. [2] present a tabu search algorithm for the periodic vehicle routing problem with multiple trips and accessibility restrictions such that not every vehicle can visit every customer. When tested on randomly generated test cases, it performed reasonably well with regards to solution quality. Furthermore the computation time was manageable for instances up to 1000 orders. We refer the reader to [15] for more literature on MPVRP.
The multicompartment VRP, and the multicommodity VRP concentrate on delivering different types of commodities to the customer either using a single vehicle, or by splitting them to several vehicles, thus requiring multiple visits to the same customer. Mirzaei and Wøhlk [47] conducted research on two variants of the MCVRP, one concentrates on split deliveries for different commodities, and the second focuses on delivering all commodities by a single vehicle. They proposed a branch-and-price method and compared the optimal costs of the two variants. The computational results were presented for instances with up to 100 customers, and the algorithm optimally solved instances with up to 50 customers and four commodities. Heuristic examples include [16] who proposed an iterated local search algorithm for solving the multicommodity multitrip VRP with the objective of minimizing the number of used vehicles. In [32] the authors addressed the commodity constrained split delivery VRP, where multiple commodities can be mixed in a single vehicle while satisfying the capacity constraint, and similar to our problem, each customer can be visited more than once, and each visit should deliver only one commodity type. They proposed a heuristic based on adaptive large neighborhood search (ALNS) and tested their approach on benchmark instances. Among the metaheuristic methods applied to the MCVRP, genetic algorithms are the most common so far. One example can be found in [69], which describes a VRP encountered in frozen food delivery. Similar to our model, they associated a penalty cost for late delivery based on the types of products. A GA was proposed for solving the model on instances with real data.
Our problem involves loading-unloading operations from the depot to a customer location, similar to the classical VRP. The difference in our case is that a single customer demand can be split into several requests, if the customer requires more than one machine type, thus several visits by different vehicles may be required. This draws a similarity to the class of VRP problems with split deliveries (SDVRP) [26]. In the SDVRP, the constraint that each customer is visited by only one vehicle is relaxed, and thus customers' demand can be split between several vehicles for delivery. The first heuristic approaches to solve the SDVRP were introduced by Dror and Trudeau [26,27]. After these studies, most of the subsequent work focused on metaheuristic or hybrid schemes. One example is the work of Archetti and Speranza [3] who applied a tabu search algorithm, and Boudia et al. [13] used a genetic algorithm combined with a local search procedure. Hybrid algorithms have since grown in popularity. Examples are found in [4,18]. Many exact models have also been proposed for this problem and one example is the study in [5]. For further literature on SDVRP we refer to [3].
Recently, attention has been paid to a class of VRPs that model multiple vehicle trips under the name "multiple trips vehicle routing problem" (MTVRP). A clear improvement can be obtained by allowing a single vehicle to perform multiple trips, especially when the vehicle capacity is limited and the replenishment of stock is required. The problem assumes that trucks can visit the depot more than once in the time horizon of the problem. Example of studies that modeled multiple depot visits are [61,63]. A survey paper on the literature of MTVRP for further details can be found in [16].
Finally, we mention the workforce routing and scheduling problem (WSRP), which is the focus of the second part of our model that involves the routing of the technicians for the installation of machines at the delivery points. There are some scenarios where personnel are required to perform tasks in certain locations, and hence require some sort of transportation to these locations. These scenarios are known as the workforce routing and scheduling problem (WRSP), which has become a widespread term used by many service providers.
A similar problem to the WRSP, is the service technician routing and scheduling problem (STRSP), which involves designing the least cost routes for vehicles carrying a number of service technicians. Each task demands the allocation of a technician with the required skills set, and this may also be associated with a time window. One of the pioneering publications in this field is the work of Cordeau et al. [20] who solved a real life technician scheduling problem for a large telecommunication company set as a competition by the French Operational Research Society in 2007. In this paper an adaptive large neighborhood search algorithm is implemented. In [68], the authors concentrated on a field technician scheduling problem in the telecommunications industry, and their purpose was to maximize the number of served requests as well as considering the request's priority and the technician's skill level. A local search algorithm, a Greedy randomized adaptive search procedure (GRASP) and a greedy heuristic algorithm were proposed to solve the problem. Kovacs et al. [40] studied the service technician routing and scheduling problem with the objective of minimizing the total routing and outsourcing costs. The authors used an adaptive large neighborhood search algorithm for solving the problem on artificial and real-world instances. Pillac et al. [52] proposed a parallel matheuristic approach for solving a variant of the TRSP in which a number of technicians with a set of accompanying skills, tools and spare parts need to be scheduled and routed within given time windows. The study dealt with the availability of tools and spare parts for the technicians and routing them to the depot for the replenishment of tools. Xie et al. [67] used an iterated local search algorithm to solve the TRSP. They studied a variant where it was given which technicians can serve which orders. The algorithm was benchmarked on instances ranging from 25 to 100 orders and compared to an ALNS algorithm, where it was found that it performs significantly better on large instances with fast computational times.
In addition to technician routing, similar problems can also be found in other fields where scheduling is important, such as the home health care [11], and the scheduling of security personnel [48].
The problem studied in this paper is a unique version considering its set of constraints and the integration with the staff rostering and routing problem. To the best of our knowledge there is no similar version investigated in the literature, and the best matching study to our problem description is by Bae and Moon [8] where they extended the multidepot vehicle routing problem with time windows to a problem of service vehicles used for delivery and installation of electronics. They developed a mixed integer programming model, a heuristic and a genetic algorithm, and compared their performances. There are differences between this study and our problem, for example we consider a longer planning horizon, and deal with multiple types of machines. We also allow multiple visits to the customer by different vehicles, while their model only allows a single visit for both delivery and installation. They also assign a maximum period of time between delivery and installation that must not be exceeded, while in our model we restrict this time by imposing a penalty.

Selection hyper-heuristics
Hyper-heuristics are general purpose search methodologies for solving difficult combinatorial problems. They operate at an abstract and higher level than heuristics, that is, over the low level heuristics space [25]. One of the earliest hyper-heuristic frameworks proposed requires that hyper-heuristics should not use any specific knowledge from the solution domain [21], a feature which makes them applicable to problem instances with different characteristics or even different problems, without a need for further algorithmic or parametric adjustments. This feature forms a principle concept of hyper-heuristics in past and modern research. Hyper-heuristics are defined as "heuristics to choose heuristics" in [21], although the first attempt to design hyper-heuristics dates to as early as 1963 [28]. Burke et al. [14] identifies two categories based on the nature of the heuristic search space: selection and generation hyper-heuristics. In the former class, a heuristic is selected from an existing repository of heuristics to try to discover the behavior of these heuristics in order to enable/disable some of them during the search process; while in the latter, new heuristics are built by discovering the characteristics of the input heuristics. The approach in this study uses the former type of hyper-heuristics, which are based on a single-point based search framework with two consecutive operations that work iteratively to improve a single initial solution through heuristic selection and move acceptance. With the existence of a defined number of low level heuristics, the selection method chooses one of these heuristics and applies it to a solution in hand, generating a new solution. The move acceptance decides on the acceptance of the new solution based on the fitness/objective evaluation. Heuristic selection can be carried out using simple methods such as random selection or by selecting from a predefined ordering of the low level heuristics, or it can incorporate learning by defining some probability measures for the performance of low level heuristics. For the most recent advances and classification of selection hyper-heuristics we refer to [25].
There are different criteria to classify selection hyper-heuristics, and one of them is the solution nature, where selection hyper-heuristics are classified based on this measure into single point or multiple point. Multiple point (population) selection hyper-heuristics utilize multiple current solutions during the search, while single point (single solution) hyper-heuristics are based on a single current solution that is iteratively improved during the search. The majority of the previous studies on selection hyper-heuristics present approaches based on single-point-based search, and only a few used a population of solutions or a mixed approach alternating between using single and multiple solutions for the search. Moreover, those previously proposed population based approaches are mostly a hybrid between a selection hyper-heuristic and an evolutionary algorithm framework.
Cowling et al. [22] investigated a genetic algorithm based on hyper-heuristics for the personnel scheduling problem. A GA is implemented and applied as a high level selector, and a set of low level heuristics are used at each generation to locally improve the quality of each individual, where the low level heuristics are applied in any sequence. Sabar and Kendall [58] proposed a Monte Carlo tree search hyper-heuristic framework that tries to identify good sequences of heuristics using the Monte Carlo search tree. A memory mechanism containing a population of solutions is utilized, and at each iteration a solution from the population is selected, and the population is subsequently updated using several updating rules. Lei et al. [43] proposed a memetic algorithm based on hyper-heuristics to solve an examination timetabling problem. Their approach constructs several heuristic lists based on graph coloring heuristics and applies evolutionary operators to generate new lists. A local search method is used to further optimize the solutions. Hsiao et al. [33] implemented a hyper-heuristic based on variable neighborhood search (VNS) iterating in two stages, first using a population of solutions, and the second stage uses only a single solution. Their approach consists of two main steps, shaking and local search. The shaking phase improves the exploration of the search space, and the local search step looks for the local optima. A population of solutions is used in the shaking stage, where the authors argued that the diversity of solutions is important in the first stages of the search to explore the right search path, and after a period of time the best solution is picked from the population. Tournament selection is used to filter unfit solutions from the population. Lehrbaum and Musliu [42] introduced a hyper-heuristic that alternates between working on a single solution and a population of solutions. Their algorithm starts by scoring the available local search heuristics, and a serial phase working with single solutions starts by applying the heuristics sequentially according to their quality scores. A parallel phase uses a population of solutions, and a heuristic is applied to each individual in the population. The algorithm switches back to the serial phase whenever a global improvement is found (i.e., better than the best found solution so far) .
We also discuss here the successful application of hyper-heuristics in different variants of VRP problems, since the VRP is considered an active research field in hyper-heuristics. In [53] the adaptive very large neighborhood search hyper-heuristic was successful in finding the state-of-the-art results for many benchmark instances of several variants of VRP. Also the CVRP with time windows is one of HyFlex (hyper-heuristic flexible framework) problem domains, for which newly developed approaches in hyper-heuristics are tested [65]. Other VRP problems solved using selection hyper-heuristics include the periodic VRP [17], dial a ride [64], urban transit routing [1], and inventory routing [34].
The previous VeRoLog challenge 2016-2017 tackled a rich VRP problem related to a cattle improvement company that regularly measures the milk quality at a number of farms using specialized tools. These tools have to be delivered to a number of farms (customers) on request and picked up again a few days after delivery. The key challenge is how to schedule the deliveries to satisfy the requests, while at the same time design efficient routes for the pick-ups and deliveries. The second place winner on this challenge used a hyper-heuristic approach based on an online selection method [38].

DESCRIPTION OF THE PROBLEM
The real-world problem from VeRoLog Solver Challenge 2019 can formally be stated as follows (see Figure 1). There are a number of customers Cr = {cr 1 , cr 2 , …, cr |Cr| } geographically spread at different locations, and a depot located at H 0 . The distance between any two locations i and j is given by d i, j . The purpose is to respond to customers' requests by delivering machines and getting them installed by a technician within a defined time horizon T given as a consecutive number of days.
An unlimited number of identical trucks (i.e., vehicles) K = {k 1 , k 2 , …} can be hired to transport the machines to the customers. They are located at the depot each with a maximum capacity C. Also, a number of machines M = {m 1 , m 2 , …, m |M| } are available to be delivered to the customers at their request, and there are different types of machines with different sizes expressed in the same size unit as the truck capacity. The machines are all located in the depot, with enough machines to satisfy all the demand. A set of customer requests R = {r 1 , r 2 , …, r |R| } should be satisfied. The requests are known at the start of the planning period. A single request r i = {cr i , w i , m i , n i } asks for one machine type m i ∈ M, of quantity n i , for exactly one customer cr i , and w i is the associated time window to deliver the request, where w i is specified by the earliest day e i and the latest day l i to deliver request r i . A request of the same type of machines cannot be split and should be delivered by the same truck, and if a customer requires another machine type, a separate request is made. Each truck journey on a day should start and end at the depot location H 0 and can carry different types of machines to satisfy several customers' requests, and each request occupies c i of the truck capacity, and should not exceed its maximum capacity C. The truck can return back to the depot location multiple times during the day to pick up more machines. Also, there is a limit D on the maximum distance a single truck can travel per day. It does not take any time to load a machine at the depot or to unload a machine at a customer.
There is a fixed number of technicians S = {s 1 , s 2 , …, s |S| } who are responsible for installing the delivered machines, at the customer location, at least 1 day after the delivery. Each technician s i ∈ S is located at a certain home location H s i . A technician's daily route starts and ends at his/her home location, and like trucks, there is a maximum distance D s i the technician can travel per day. In addition, there is a maximum limit N s i on the number of requests a technician can carry each day, where carrying out a request means installing all the machines for that request. The technician can maximally work for five consecutive days, and must have 2 days off if he/she has worked for five consecutive days.
Each technician has a skill set for installing certain types of machines. a sm refers to technician s ∈ S installing machine m ∈ M, and is equal 1 if the technician is eligible to install this machine, and zero otherwise. Installing a machine does not take any time. A technician is described with the following entry s = {H s , D s , N s , {a sm 1 , a sm 2 , … , a sm |M| }} referring respectively to the technician home location, maximum travel distance per day, maximum number of installations per day, and which machines they have the skill to install.
The last point to mention, is that the technician should install a delivered machine as soon as possible after the delivery, and for each delayed installation of request r i , a penalty Cl i is added to the cost, and each machine type has a different penalty value.
The main objective is to reduce the overall costs associated with trucks, technicians and idle machines costs. The trucks/technicians total cost is constituted of the following parts: the cost of hiring a truck/technician per day, the cost of hiring a truck/technician within the planning horizon T, and the cost per unit of distance for the traveling of truck/technician. The distance between coordinates (x 1 , y 1 ) and (x 2 , y 2 ) is defined as the ceiling of the Euclidean distance, In addition, there is the cost for penalizing idle machines that remain without installation for more than 1 day. This penalty cost is dependent on which machine it is and the number of days it was idle. The objective function is described by Equation (1) in the following section.
A solution gives, for each day in the planning horizon, the routes followed by each truck/technician (see Figure 2). Assuming that {1,6,7,0,1,2} is the route of a single truck in one of the planning days. The first element in the route "1" is the truck ID, followed by the requests ID's that this truck served. The ID "0" refers to the depot, and it means that truck 1 visited the depot after serving requests "6" and "7" and was loaded to serve requests "1" and "2." Each series of requests before the truck goes back to the depot is named "tour". In this route, there are two tours given as {6, 7} and {1, 2}. The start and end of the truck journey at the depot is not explicitly written in this route format.
The technician routes are very similar, starting with the technician ID, followed by the ID's of the requests that this technician served. Also, the start and end of the technician journey at their home location is not explicitly mentioned in the solution.

Problem instances
We have used two datasets of instances, one of them has been developed for this work and one, referred to as the hidden dataset, was used in the competition to evaluate the participants' algorithms in the VeRoLog Solver Challenge. Each of the instances provides different types of information such as the weights of the objective function components, the maximum truck capacity, the number of days in the planning horizon, and the maximum travel distance allowed by each truck. The details of the requests, locations given as x, y coordinates (i.e., depot, technicians homes, customers) and technicians are also given. The characteristics of these datasets are provided in Table 1.
The small dataset, which includes instances of sizes varying between 6 and 16 requests, is developed specifically for this work. The reason for generating this dataset is to provide an ideal size of instances for testing the mathematical model which can only be applied on instances of such sizes. It is also essential to test our developed hyper-heuristic approach on instances with different characteristics and scales and to compare its performance to the exact model by its ability of finding optimal solutions in a short duration of time.
The hidden dataset was used to assess the performance of the competitors algorithms and rank the finalists in the restricted resource challenge. 2 This dataset contains instances of large sizes up to 900 requests. The number of different types of machines vary between 3 and 7 in each instance, and the number of technicians range from 25 to 125. The highest variation can be found in the costs of using trucks and technicians, distance costs, and the costs per day that trucks and technicians are used. These values range from 10 to 100 000. We refer the reader to [31] for a comprehensive description of the problem and the formal challenge rules. 3

Sets and indices
. S: technicians (s ∈ S). R s : requests that technician s can install and the home location of technician s (R s = {s, i : a si = 1, ∀ i ∈ R} ∪ H s ).

Parameters
T: number of days in the entire planning horizon.

Mathematical modeling
∑ k∈K,j∈R 0 ,i≠j ∑ i∈R s ,j∈R,i≠j ∑ s∈S,j∈R s ,i≠j The exact model for the given problem is formulated by Equations (1)- (27). The objective function (1) is composed of three types of costs: the vehicle (first three summations), technician (next three summations), and idling cost (last summation). The sum of the vehicle hiring cost, the vehicle cost per day, and the vehicle cost per distance is equal to the total vehicle cost. Similarly the sum of the personnel hiring cost, the personnel cost per day, and the personnel cost per distance is equal to the personnel cost. The idling cost is calculated by multiplying the cost of idling all the machines at each request by the difference between the delivery and installation days.
In the mathematical model, constraints (2) Constraints (2) and (11) are balance equations for the vehicles and technicians respectively. They ensure that in any day, in any vehicle and personnel route, the number of arcs entering to a location of a request, depot or home should be equal to the number of arcs exiting from the same location.
Constraints (3) and (12) ensure that both the vehicles and technicians start/end their routes from/at the depot and home respectively. For the technicians, the problem dictates upper bounds for the number of installations (N s ) and travel distance (D s ) for any given day. Therefore, in an optimal solution, if a technician is used on any given day, they should leave and return their home only once. On the other hand, for the vehicles, in addition to the total travel distance on any given day (D), there is an upper bound on the weight of the machines that are being carried by the vehicle, defined as vehicle capacity (C), at any given time of the day. This makes it possible for a given vehicle on any given day to deliver more than its capacity by making multiple visits to the depot. Because of the difference of the restrictions on the vehicles and technicians, constraints (12) are satisfied with an equality to p t s whereas constraints (3) are satisfied with an inequality to v t k . The former constraints force technicians to leave and return to their home only once if they are working on day t, whereas the latter constraints force vehicles to have an equal number of trips that leave from and return to the depot, and these trips can only occur if the truck is operating on day t. M t on the RHS of the constraints (3) is calculated by counting the number of orders that can receive a delivery on day t. Note that, the maximum value that M t can take on different days also sets the upper bound for the maximum number of hired vehicles in the planning horizon.
Constraints (4) and (13) restrict the total travel distance for each day of the vehicles and the technicians, respectively. In addition, constraints (14) set the maximum number of installations for a technician on a single day. Similarly, in the problem definition, there is a limit set on the total weight of machines a truck can carry. This is ensured by constraints (9) and (10). The variable q H 0 k , that represents the weight of the machines right after the vehicle is leaving the depot, is set to be C for any truck by constraints (10). Since constraints (9) ensure that the weight on the truck decreases at every request stop by the weight of the delivery of the same request, and the q ik s are defined as nonnegative integer variables, no truck can carry more than its capacity.
Constraints (5) and (15) ensure that in order to use a vehicle or a technician respectively in any day of the planning horizon, we need to hire them first. Similarly, constraints (6) and (16) ensure that if a vehicle or a technician travel between two requests on any given day, they are already hired for the day.
Constraints (7) and (17) establish the relationship between the routing (x t ijk and z t ijs ) and service, that is, delivery (w t i ) and installation (y t i ) variables for the vehicles and technicians respectively. According to constraints (7), if the location of a request is visited by a vehicle, then the machines ordered by this request are delivered between the first (e i ) and the last (l i ) days the delivery can be done. Similarly, according to constraints (17), if the location of a request is visited by a technician, then the installation that is ordered by this request is done at least 1 day after the earliest day that the request can be delivered.
Constraints (8) and (18) ensure that all the deliveries and installations are done within their predefined times respectively. Constraints (8) ensure that each request is delivered by one of the vehicles. These constraints restrict each request to be delivered between their first (e i ) and the last (l i ) days the delivery can be made. Similarly, constraints (18) ensure that each request is installed by one of the technicians. Since the earliest installation day is 1 day after the delivery of the machines, the constraints only consider the days after the first day that the request can be delivered.
Constraints (9) and (10) prevent subtours in vehicle routes. Constraints (10) assign the maximum weight (C) a vehicle can carry when leaving the depot to variables q ik on any given day. Constraints (9) subtract the weight of the machines that are being delivered from q ik every time a vehicle visits a request. Keeping track of each q ik and forcing q jk less than or equal to q ik if request j is visited immediately after request i or depot by vehicle k (i.e., x t ijk = 1) prevent the formation of subtours. Note that, since a vehicle can make multiple visits to the depot on any given day, these constraints are not being forced for the trips to the depot.
Constraints (19) and (20) prevent subtours in technician routes. Constraints (20) assign the maximum number of visits a technician can make (N s ) to variables g is on any given day. Note that, since the maximum number of installations that a technician can make is an upper bound for the maximum number of visits in a day, we use this constant in our model. Constraints (19) subtract 1 from g is every time a technician visits a node. Similar to vehicle subtour elimination constraints, since g js is forced to be less than or equal to g is if request j is visited immediately after request i or home H s by technician s (i.e., z t ijs = 1), subtours never form in technician routes.
Constraints (21) and (22) ensure that the solution complies with the working day restrictions for the technicians. According to constraints (21) and (22), if a technician works four or fewer consecutive days, that is, LHS is less than or equal to 4, he/she can work either the next day or the day after the next day unless he/she is working more than five consecutive days. If a technician works five consecutive days, that is, LHS is equal to 5, the technician cannot work the next two consecutive days. Constraints (23) ensure that this restriction still holds at the end of the planning horizon and prevents any technician from working more than 5 days in the last 6 days.
Constraints (24) calculate the idling time, the difference between the delivery day and the installation day, for all requests. Since b i is defined as a nonnegative integer, these constraints also ensure that machines are installed at least 1 day after the delivery for every request.

HYPER-HEURISTICS METHODOLOGY OF CVRP FOR DELIVERY AND INSTALLATION OF MACHINES
In this section we describe the hyper-heuristic framework applied to this problem and discuss solution initialization and representation and the low level heuristics set.

Population-based hyper-heuristic framework
Following the description of the selection hyper-heuristic framework in Section 2.2, most of the previously proposed solution methodologies in selection hyper-heuristics utilize a single solution during the search process and iteratively improve it, while some other methodologies adopt the idea of using a population of solutions during the search as a whole, or during some part of it. Our proposed framework is based on a population of solutions from which one of them will be selected and applied to a selection hyper-heuristic at each step in the search. We are motivated in this work to use a population of solutions as we believe that this provides diversity in the search and allows better exploration of new areas in the search space. The process starts by initializing a number of solutions using a generation method to create an initial population pop = {sol 1 , sol 2 , … , sol pop size }. A number of selection hyper-heuristics HH = {hh 1 , hh 2 , …, hh n } combining different selection and move acceptance methods are implemented.
A solution sol i and a selection hyper-heuristic hh j are randomly selected from pop and HH respectively, where sol i will serve as S current to hh j . The selection hyper-heuristic hh j selects a heuristic (or sequence of heuristics) and applies it to S current to create new solution S new , which is checked for feasibility, and rejected if it is not feasible (i.e., violates at least one of the problem constraints). If S new is feasible it will be evaluated and the decision of its acceptance is made by the move acceptance component of hh j . The best found solution S best is replaced by S new if it is better. The iteration between the selection and acceptance components continues until the termination criteria are satisfied, that is until the global time limit is exceeded or when there is no improvement on the best obtained solution for a certain number of iterations.
After hh j terminates, S best is checked against the best found global solution S global and replaces it if it is better. Afterwards, S best is shuffled by randomly selecting and applying a series of low level heuristics. The number of steps to shuffle a solution is tuned by the user, and is constant during the search. This shuffling is necessary in order to avoid the possibility of early convergence and to refresh the population. Next, a solution from the population and a selection hyper-heuristic are randomly selected and the same steps mentioned above are repeated. This iterates until the specified time for running an instance passes. Algorithm 1 outlines our applied framework (Figures 1 and 2).
For this framework we have tested a total of eight selection hyper-heuristics combining the selection methods: simple random (SR), sequence-based selection hyper-heuristic (SS), and the move acceptance methods: record-to-record (RR), Naïve acceptance (Naïve), great deluge (GD), and simulated annealing (SA).
SR is the most basic heuristic selection method. It randomly selects a low level heuristic at each step according to a uniform probability distribution. SS on the other hand applies sequences of heuristics to the solution, instead of single applications. This selection method learns and identifies the sequences of heuristics most likely to improve the current solution. More details of this method can be found in [35,36].
The move acceptance methods applied accept nonworsening solutions, and worsening solutions are accepted using different criteria.
In Naïve and SA, the worsening solutions are accepted with a certain probability. This probability is fixed for Naïve which is predefined by the user, while in SA, the probability varies in time and it is calculated using the following formula: where Δf is the change in the cost at time t current , ΔF is the expected maximum change in the cost, and t limit is the time limit. GD is a threshold move acceptance method which allows worsening solutions if their cost (objective) value is less than or equal to a threshold value referred to as "level" ( t ) which gets updated at each time step (t) during the search. The level is initially set to the initial cost. The update formula for the level is as follows: where f is the final expected cost value, ΔF is the maximum expected change in the cost value, and t current , t limit is the time at the current step, and the time limit respectively. RR is a variant of GD which accepts worsening solutions that are not much worse than the best solution in hand to an extent based on the following formula: where fr is a factor that is updated during the search, starting with a large value and gradually decreasing.

Solution representation scheme
The described hyper-heuristic framework requires initializing a number of solutions to build a population, and this is achieved using an initial generation method. The structure of each of the initialized solutions is demonstrated in Figure 3. Each solution is composed of two main components: truck visits, and technician visits. The truck visits component corresponds to the schedule of trucks during the planning horizon which can be modeled as four levels: days, trucks dispatched on each day, tours performed by each truck, and the requests to deliver on each tour. Similarly, the technician visits correspond to the technicians' schedule composed of three levels: days, technicians scheduled on each day, and visits performed by each technician. The initial generation method randomly produces these schedules, while ensuring the final constructed solution is feasible. The main focus of the initial generation method is the feasibility of the solution and not its quality. The feasibility of the solution must also be maintained during the hyper-heuristic operation. A feasibility test is implemented to ensure that the constraints described in Section 3 still hold after each application of low level heuristic(s). A single violation

Low level heuristics
The hyper-heuristic controls a total of 25 low level heuristics to improve the quality of a given initial solution. These low level heuristics perform swap and insert operations for requests in truck and technician routes (see Figure 4). Low level heuristics are restricted, as needed, to only produce routes that respect some of the constraints. For example, some low level heuristics perform operations between different days in the planning period; in this case if the operation involves delivery requests, the time windows of these requests must be respected and any installations that as a consequence violate the time windows constraints must be rescheduled. If it involves installation requests, the delivery of these requests must be ensured at least the day before.
• LLH0: selects a random day, a random truck route and a random tour, and swaps any two randomly selected requests on this tour.
• LLH1: selects a random day, a random truck route, and two different random tours, and swaps any two randomly selected requests on each tour. • LLH2: selects a random day, two different random truck routes, two random tours from each route, and swaps two randomly selected requests from each tour. • LLH3: selects two different random days, two random truck routes from each day, and two random tours from each route, and swaps two randomly selected requests from each tour. • LLH4: selects a random day, a random technician scheduled on this day, and swaps two randomly selected requests of this technician. • LLH5: selects a random day, two different random technicians scheduled on this day, and swaps two randomly selected requests of these technicians. • LLH6: selects two different random days, and two random technicians, and swaps two randomly selected requests of these technicians. • LLH7: selects a random day, a random truck route, a random tour, and two random positions on this tour. The request on the first position is inserted into the second position. • LLH8: selects a random day, a random truck route, two different random tours on the selected route, and a random position on each tour. The request on the first position of the first tour, is inserted into the second position of the second tour. • LLH9: selects a random day, two different random truck routes, a random tour on each route, and a random position on each tour. The request on the first position is inserted into the second position. • LLH10: selects two different random days, a random truck route on each day, a random tour on each route, and a random position on each tour. The request on the first position is inserted into the second position. • LLH11: selects a random day, a random technician scheduled on this day, and two random positions on the technician route. The request on the first position is inserted into the second position. • LLH12: selects a random day, two different random technicians, and a random position on each technician route, and inserts the request on the first position into the second position. • LLH13: selects two different random days, a random technician on each day, and a random position on each technician route. The request on the first position is inserted into the second position. • LLH14: selects a random day, two different random truck routes, and a random tour on each route, and swaps the two selected tours. • LLH15: selects two different random days, a random truck route on each day, and a random tour on each route, and swaps the two selected tours. • LLH16: selects a random day, two different random truck routes, and a random position on each route. The tour on the first position is inserted into the second position. • LLH17: selects two different random days, a random truck route on each day, and a random position on each route. The tour on the first position is inserted into the second position. • LLH18: selects a random day, two different random truck routes, and two random positions. A block of consecutive requests starting at the first position is swapped with another block of requests starting at the second position. The size of the block is randomly selected. • LLH19: selects two different random days, a random truck route on each day, and two random positions on each route.
A block of visits starting at each of the positions are swapped with each other. • LLH20: selects a random day and two different random technicians, and swaps two blocks of requests for these technicians. • LLH21: selects two random different days and two random technicians from each day, and swaps two blocks of requests of these technicians. • LLH22: selects a random day, and two different random truck routes. A block of requests is moved from the first route to the second route at randomly selected positions. • LLH23: selects two different random days and two random truck routes. A block of requests is moved from the first route to the second route at a randomly selected positions. • LLH24: selects a random day and two different random technicians, and moves a block of requests from the first technician to the second technician. • LLH25: selects two different random days and two random technicians, and moves a block of requests from the first technician to the second technician. Number of iterations to shuffle solutions 10 7

EXPERIMENTAL RESULTS
The experiments were performed on several machines. CPLEX 12.10 and Gurobi 9.0 together with C# and Python were used for the exact technique and Microsoft Visual Studio 2017 C++ for the heuristic method. The experiments for the heuristic method on the hidden dataset were performed on a device with the specifications: Intel Core i5 at 2.3 GHz with memory of 8 GB. On both datasets, the experiments were designed according to the competition rules, which required nine runs per instance with nine different random seeds also determined by the competition rules. The run time for each instance in both datasets was also calculated according to the competition rules, where it has been specified that each instance is run for a limited time on the user machine calculated with the formula: T limit = fb × (10+| R| ), where T limit is the time limit for running an instance according to the user local machine, fb is a factor calculated by a benchmark tool provided by the competition to estimate the equivalent time on any machine compared to the organizers core machine, and |R| is the number of delivery requests in the instance.
To tune these parameters we have followed two approaches: a manual approach where a series of extensive experiments were performed to fine tune the design parameters. We arrived at a combination of parameter values that resulted in a relatively better performance across a subset of public instances. The second approach is using the irace package [44] to automatically tune the parameters on the five instances with the highest variance in the small dataset. Irace performed a maximum of 8000 experiments and ran for about 9 hours to find the top four configurations. The best configuration in the top four was selected to perform another round of experiments on the small dataset with the same experimental design (i.e., nine runs per instance with the random seed values set by the competition rules). The parameter settings of the hyper-heuristic using the two approaches are shown in Table 2. The results of the small dataset reported in the next section are the ones found using the manual tuning. For convenience, the developed approach is denoted as POHH. Table 3 provides the results of the small instances dataset for the exact and the hyper-heuristic method. For the exact model, the upper and lower bounds are provided for each instance. The lower bound indicates that an optimal solution was not found for a particular instance. The results of the hyper-heuristic experiments are reported in terms of the minimum and maximum objective values achieved in the nine runs, the average of the nine runs and the SD. The time in seconds is the time that was required to find the reported results by the exact model, and the time limit for each run of POHH. The time was normalized to its equivalent in the standard machine using the calibration tool provided by the competition. We also note that the execution time of the exact model on any instance was limited to up to 30 minutes. From Table 3 we can directly compare the performances of the two models in terms of finding optimal solutions and the time required to find them. The exact model was able to find optimal solutions for 12 instances out of the 25 with CPLEX and Gurobi, and no feasible solution with CPLEX or Gurobi was found for the instance Small_09. For the rest of the instances the same upper bounds were found by CPLEX and Gurobi, while we notice that for some instances Gurobi was able to find better lower bounds. Comparing the results of the exact model to the minimum value in the nine runs, POHH was able to find the optimal solutions in all 12 instances where the exact model found optimal solutions. In the other cases where no optimal solution was proved by the exact model, the POHH algorithm either found the same upper bound or, in the case of seven of the instances, a better upper bound was discovered. Also, a feasible solution for Small_09 was found by POHH. Although the POHH was able to find the same value for the upper bound as Gurobi and CPLEX in many instances, we cannot yet argue that this upper bound is the optimal solution to these instances. In terms of run time, POHH achieved improved run times in most of the cases. The exact model in many instances required more than 3000 seconds to find a solution, while POHH required less than 30 seconds on the same instances.

Results on the small dataset
We also compared the results on the small dataset using the manually tuned parameters with the best configuration found by the irace package. The results of the two approaches were very similar using the minimum value in the nine runs as comparison Table 3 Summary of the results base. Both found the same minimum in all the instances, except for Small_18, in which the irace configuration found a new best result of "479 817 740", and Small_23, in which the manual configuration was better.

Results on the hidden dataset
As mentioned previously, the hidden dataset was used to assess the competitors' algorithms in the restricted resources challenge, and according to the results of this challenge, eight teams were selected as finalists. We have followed the same competition rules as the other finalists team, with this set of experiments to fairly compare and justify our results. Table 4 summarizes the results of the top six teams, including our hyper-heuristic approach, for each instance ranked based on their best found solution from best to worst. For each instance, the results are reported for the best six teams out of nine using the average of the nine runs and the minimum objective value.
Considering the minimum and average objective values obtained over nine runs for each hidden instance, POHH is relatively competitive between the finalists. The POHH achieved a position in the top six in 15 instances out of 25. For 11 out of 25 instances, including: Hidden_02, Hidden_07, Hidden_09, Hidden_11, Hidden_13, Hidden_16, Hidden_17, Hidden_19, Hid-den_21, Hidden_22, Hidden_25, POHH performs better than at least half of the finalist approaches in terms of average objective value. Except the instances Hidden_13 and Hidden_17, the same phenomena is observed with those instances with respect to the minimum objective values. The best achieved results are found on instances Hidden_07, Hidden_11, Hidden_16, and Hid-den_21, where POHH is ranked the fourth based on both the average and the minimum objective values. These instances are all of different sizes: 150, 450, 600, and 750 requests, reflecting the ability of POHH to perform well on instances with varying characteristics and complexities ( Figure 4).
We have also ranked our approach among the eight finalists using the same method used in the competition. A ranking score is calculated per instance for each submitted solver by removing the two best and worst solutions found by this solver. We then take the average objective value of these five solutions as score for the algorithm, and rank all methods accordingly. The average of all ranking scores for the instances represents the final mean rank of the solver, which was used to order the competitors from the first to the last. Figure 5 displays the ranking of the POHH algorithm among the eights finalists based on this method. It is clearly seen that our method was able to produce results competitive with the finalists by achieving a better final mean rank than three teams, and an insignificant difference from the ranks of the third, fourth, and fifth teams.
Although the proposed algorithm did not succeed in improving any current best known solutions on hidden instances, it performed well (see Section 6.1) on small instances with few requests, and the results on the hidden instances are considered reasonably good. Figure 6 visualizes the six sample instances, including Small_01, 03, 06 and Hidden_01, 03, 06 that we used for the analysis purposes, reflecting the varying characteristics of each instance. These instances were arbitrarily chosen to represent varying sizes. The visualization was obtained with the tools provided by the challenge organizers [51]. The large light blue circle indicates the depot, and this may or may not have technicians on the premises at any given time.  Figure 6 also shows semi-transparent green circles centered on the locations with technicians. The radius of these circles is equal to the half of the maximum daily distance of the corresponding technician. Because these circles are semitransparent, overlap leads to color intensification, making visually clear which customer locations are within reach of few or many technicians. One may notice from the clustering, for example, as for Hidden_06, that the locations are making the shape of the Netherlands. This implies that these are indeed based on real-world data offered by ORTEC.

Performance analysis of POHH
Although Figure 6 gives a rough picture of those instances (e.g., some instances are more limiting in terms of number and action radius of technicians), we must emphasize that it does not describe a given problem instance fully. For example, the importance of violating a given constraint is not depicted and, as we mentioned before, penalties for violating the different constraints can differ substantially as a part of the cost function. Table 1 is particularly useful in this case.
The pie charts in Figure 7 depict the utilization rates of the different selection hyper-heuristics applied in our framework. The utilization is calculated in terms of the ratio between the number of times a selection hyper-heuristic was successful in finding an improved solution over the best global solution to the total number of improvements made by all the selection hyper-heuristics in the duration of run time. Table 4 The performance comparison of the finalists and POHH, based on the average, and minimum of the objective values over nine runs for each instance  There are particular selection hyper-heuristic methods that clearly performed better than the rest in making improvements to the solutions during the search process. The incorporation of the RR-based selection hyper-heuristics in the POHH appears to play a key role of solving the problem in a relatively effective manner, in particular SS-RR which performed equally well in the small and hidden datasets. SR-GD was very successful in the larger size instances, where 50% of the improvement rate was achieved by SR-GD in Hidden_06 that has 900 requests. The least successful selection hyper-heuristics are the ones combined with the simulated annealing acceptance. The utilization of SR-SA was down to 0% in all instances, except for an insignificant improvement rate of 5% in Small_03. Also, SS-SA did make much contribution in terms of improvement for the hidden set. The naïve acceptance is more successful for the small instances than the larger ones, but only when it is combined with the sequence based selection method.
There seems to be a variation in the performance between the selection hyper-heuristics in instances with different complexities, and we cannot generalize that a certain combination of a selection and a move acceptance methods would be successful in every instance in this problem. An interesting idea would be to embed a high-level intelligent control mechanism that can observe these variations, and apply the components of selection hyper-heuristics accordingly during the search time, similar to the online selection methods. The random selection criteria that we apply currently in our framework was able to find "reasonable" results as reported, and we expect that the suggested improvements in the selection mechanism could yield even better results.

Performance comparison to the constituent hyper-heuristics
Another round of experiments have been conducted by applying the eight selection hyper-heuristics employed in our framework independently on the instances displayed in Figure 6. Each selection hyper-heuristic was run for nine times using the same rules to calculate the run time of an instance described in this section. The results are displayed in Table 5 using the best and average objective values from the nine runs, along with the associated SD. The Mann-Whitney-Wilcoxon test is performed with a 95% confidence level in order to compare pairwise performance variations of two given algorithms statistically. The following notations are used: (a) "+" denotes that our algorithm (POHH) is better and this performance variance is statistically significant, (b) "−" denotes that the performance of POHH is worse and this performance variance is statistically significant, and (c) "=" indicates that there is no statistical significant between the two methods. The POHH algorithm performed statistically better than each of the constituent hyper-heuristic for all instances, except Hid-den_03, where the two methods SR-RR and SS-RR found slightly better averages and performed statistically better. Other than those two cases, the POHH algorithm found the best averages and minimum values in all instances, performing exceptionally better in particular on the largest instance of Hidden_06. This provides evidence for the success of two proposals: (a) utilizing multiple solutions allows better exploration and more possibilities for further improvement in new areas of the search space, instead of focusing on a single solution; (b) applying a sequence of selection hyper-heuristics to a solution might be useful in utilizing the varying performances of these selection hyper-heuristics. This can lead the search to different directions that can yield further improvement.
Overall, POHH creates a synergy among the eight selection hyper-heuristics resulting in an improved performance for almost all instances. Although POHH does not utilize learning, various learning mechanisms embedded into each one of the low level selection hyper-heuristics potentially contributes to the overall success of the proposed approach when compared to the performance of each constituent selection hyper-heuristic.

CONCLUSION
In this study, we tackled a complex VRP problem which was the subject of the fourth edition of the VeRoLog solver challenge (2019). The challenge consisted of a novel VRP problem comprising two interdependent stages: a capacitated VRP problem with time windows (CVRPTW) for delivering various equipment to customers on their requests, and a service technician routing and scheduling problem (STRSP) for the installation of the delivered equipment. We propose two approaches, and apply them to the set of instances supplied by the competition organizers, and also to another small set of generated test instances. The first approach is an exact mathematical approach, which is the first attempt to implement an exact model for such version of VRP problem. We show that even with small sized instances, the method requires large amounts of computing time to solve the problem. Moreover we prove that the exact model cannot solve instances of large sizes. In light of this, and due to the large number of constraints and the wide range of differing instances which make it difficult to create problem-specific algorithms, Table 5 The performance comparison of POHH, SS-RR, SS-GD, SS-SA, SS-Naïve, SR-RR, SR-GD, SR-SA and SR-Naïve based on the average (Avg), associated standard deviation (Std), minimum (Min) of the objective values over nine trials and the pairwise average performance comparison of POHH vs (SS-RR, SS-GD, SS-SA, SS-Naïve, SR-RR, SR-GD, SR-SA and SR-Naïve) based on Mann-Whitney-Wilcoxon for each instance produced by each approach