key: cord-0961877-j51gaesp
authors: Amiri, Afsaneh; Salari, Majid
title: Time-constrained maximal covering routing problem
date: 2018-12-04
journal: OR Spectr
DOI: 10.1007/s00291-018-0541-3
sha: dcb2f77299e7339d5a36d3b2aae45c266ab15f37
doc_id: 961877
cord_uid: j51gaesp

We introduce the time-constrained maximal covering routing problem (TCMCRP), as a generalization of the covering salesman problem. In this problem, we are given a central depot, a set of facilities and several customers which are located within a pre-determined coverage distance of available facilities. Each facility can supply the demand of some customers which are within its coverage radius. Starting from the depot, the goal is to maximize the total number of covered customers, by constructing a set of p length constraint Hamiltonian cycles. We have proposed a mixed integer linear programming model and three heuristic algorithms, namely iterated local search (ILS), tabu search (TS) and variable neighborhood search (VNS), to solve the problem. Extensive computational tests on this problem and some of its variants clearly indicate the effectiveness of the developed solution methods.

patients receive medical treatments in the shortest possible time in order to prevent the outbreaks of disease. Therefore, the problems in which researchers are seeking optimal methods to reduce the traveling time or distance as well as increase the number of satisfied customers are considered as effective ways to deal with the spreading of epidemics. For this reason, in this paper, we discuss a problem which simultaneously increases the number of the covered customers and decreases the length of the traveling distance. Indeed, we have combined the routing and the covering problems in order to benefit from the characteristics of both problems. Several combinatorial optimization problems, in the field of routing and covering, are proposed and developed. In particular, we can refer to the vehicle routing, the maximal covering and the traveling salesman problems. In the following, the mentioned problems and their variations are discussed first and then we outline the proposed problem, its characteristics and applications.

The vehicle routing problem (VRP) is a combinatorial problem, seeking to visit a number of customers by utilizing a set of vehicles which are located at a central depot (Dantzig and Ramser 1959) . Most popular variants of this problem are the VRP with time windows (VRPTW), the capacitated VRP (CVRP) and the multi-depot VRP (MDVRP). In the VRPTW, each customer must be visited within a certain time interval. The problem has several applications including the meal delivery, the waste collection and the school bus routing problems (Jabali et al. 2015; Bae and Moon 2016; Buhrkal et al. 2012; Hiermann et al. 2016) . The vehicles are capacitated in the CVRP, and several applications, namely road transportation and network routing problems, have been introduced for this problem (Faulin et al. 2011; Sörensen and Schittekat 2013; Hosseinabadi et al. 2017) . The MDVRP uses multiple depots instead of one depot to construct the routes. Delivery of food and drugs, chemical and oil products is some applications of the MDVRP (Wasner and Zäpfel 2004; Kergosien et al. 2013; Lalla-Ruiz et al. 2016) . For further studies, readers can refer to Toth and Vigo (2001) and Golden et al. (2008) .

The maximal covering problem (MCP) was introduced by Church and Velle (1974) . The goal of the MCP is to maximize the amount of demand that can be covered by a limited number of facilities within a specified distance or time. Generally, in the MCP we are given several customers and facilities; each facility can respond some customers' demand which is located within its pre-determined distance. The goal is maximizing the number of covered customers by using a limited number of facilities. In Pirkul and Schilling (1991) the capacitated MCP was introduced in which each facility has a limited capacity.

The covering salesman problem (CSP) is an extension of the traveling salesman problem (TSP) which is seeking construction of a minimum length tour over a subset of n nodes such that each unvisited node is within a pre-determined covering distance of at least one of the visited nodes. Routing of rural healthcare delivery teams, locating aircraft for overnight delivery systems and locating mail boxes are some applications of the CSP (Current and Schilling 1989) . Gendreau et al. (1997) introduced the covering tour problem in which the vertices are composed of three sets including a set of vertices that must be visited (S 1 ), a set of vertices that can be visited (S 2 ), and a set of vertices that must be covered (S 3 ) by visited ones. The purpose is to find a minimum length cycle over the vertices of S 1 and a subset of those belonging to S 2 in order to cover all the vertices of S 3 . One of the applications is locating the post boxes among potential sites in a way that all of the users located in a reasonable distance from boxes are covered and the cost of collection tour over all selected boxes is minimized. In Golden et al. (2012) the authors provided a generalized version of the CSP. The objective is to find a minimum cost route over a subset of vertices such that every unrouted vertex has to be covered at least k times by those vertices that are visited on the tour. As an application, in routing of rural healthcare delivery teams, in each day, a limited number of people could be served by the team; therefore, some points need to be visited more than once. Another related paper is Current and Schilling (1994) in which two bi-objective problems known as the median tour problem (MTP) and the maximal covering tour problem (MCTP) are introduced. In each problem, the vehicle must visit only p (out of n) nodes on the network. Furthermore, minimization of the total tour length is considered as one of the objectives of both problems. The second objective is maximizing the access of the tours to the nodes that are not directly on the tour. This objective in the MTP is realized by minimizing the total distance traveled from nodes that are not on the tour to the nearest stop on the tour, while in the MCTP it is achieved by minimizing the demand of nodes which are not covered by the stops on the tour. These two problems have numerous real-world applications such as rural healthcare delivery, designing of distributed computer networks and design of hierarchical transportation networks (Current and Schilling 1994) .

The other related problem is the orienteering problem (OP) which is also known as the maximum collection problem and was first introduced in Tsiligirides (1984). This problem is a development of the TSP and also known as the selective TSP. Here, each vertex has an associated profit and instead of visiting all of the vertices, a subset of vertices with the maximum profit is visited, while the tour length is not allowed to exceed from a specified threshold (Laport and Martello 1990) . The team orienteering problem (TOP) is an extension of the previous problem and consists of determining m length constraint circuits. This problem has been proposed in Butt and Cavalier (1994) with the name multiple tour maximum collection problem (MTMCP) . Several variations of the TOP are presented in the literature including TOP with time windows (Vansteenwegen et al. 2009; Gunawan et al. 2015; Souffriau et al. 2013 ) and the capacitated orienteering problem in which vehicles have limited capacities (Archetti et al. 2010; Aras et al. 2011; Zachariadis and Kiranoudis 2011; Bock and Sanità 2015) . For a comprehensive survey of OP and TOP, the reader is referred to Vansteenwegen et al. (2011) and Gunawan et al. (2016) .

In another relevant problem which is referred to as the ring star problem, the goal is to minimize the length of a cycle which is constructed over a subset of vertices, as well as to minimize the allocation cost which is imposed to the problem by the remaining vertices that are not located in the cycle and served by the closest vertex in the cycle. This problem is an extension of the location allocation problem and has been introduced firstly by Lee et al. (1996) and developed and solved in Labbé et al. (2004 , Baldacci et al. (2007 , and Sundar and Rathinam (2017) .

In this paper, we present the time-constrained maximal covering routing problem (TCMCRP) which is a generalization of the time-constrained maximal covering salesman problem (TCMCSP) (Naji-Azimi and Salari 2014). In this paper, we develop the TCMCSP by considering multiple vehicles instead of just a single vehicle due to the fact that in real world, we have to deal with the problems with several vehicles. In the TCMCRP, we are given a set of vertices including a central depot, a set of customers Central Depot Customer Facility Coverage Radius Fig. 1 An illustrative example of the TCMCRP and a set of facilities. Each facility can supply the demand of some customers within its pre-specified coverage distance. The goal is maximizing the total number of covered customers by constructing a limited number of length constraint Hamiltonian cycles over a subset of facilities.

Since the TCMCRP is the combination of the CSP and the TOP, it can be applied to all applications already proposed for these two problems. For example, in rural healthcare delivery system, teams should travel to numerous areas in a limited time to provide general care services like vaccination, emergency services or such activities for prevention of epidemics. Having a limited amount of resources, the teams have to select some areas for visiting in order to maximize the number of individuals that will be covered. Another application is discussed when natural or human disasters like earthquakes, wars or typhoons occur. In these situations, the problem includes the assignment of some locations for establishing field hospitals among the potential areas. The locations should be chosen in a way that the number of injured people who can receive emergency assistance is maximized and the traveled distance from main centers to these locations is minimized. In above examples, maximizing the covered demand subject to the availability of limited resources is the main issue which is the basic characteristic of the TCMCRP model.

An illustrative example is shown in Fig. 1 . In this figure, we are given 41 nodes including 10 facilities, 30 customers and a central depot. In this example, 6 facilities are visited by two available vehicles and totally 22 customers have been covered.

The remainder of this paper is organized as follows: Formal description of the problem is provided in Sect. 2. Section 3 describes the characteristics of the heuristic algorithms. Computational results are presented in Sect. 4, followed by conclusions which are reported in Sect. 5. 2006; Montané and Galvao 2006) . In our developed model, the sub-tour elimination constraints are a generalization of those proposed in Karaoglan et al. (2009) and Kara (2011) for the single-vehicle version of the problem. In particular, the proposed model is a generalization of the already proposed model in Naji-Azimi and Salari (2014) for the TCMCSP.

This TCMCRP is defined on a directed graph G (V , A) where V {0} ∪ T ∪ W is the set of vertices. {0} represents the central depot at which a set of p homogeneous vehicles P {1, 2,…, p} are located. W is the set of customers, and T is the set of facilities. A {(i, j)|i, j ∈ T ∪ {0}} is the arc set where each arc is associated with a weight t i j (travel distance or time). The maximum travel length of each vehicle k ∈ P cannot exceed the given limit L k .

Covering matrix D d i j is given in which d uv 1 when customer u is located within a pre-specified coverage radius of facility v; otherwise, d uv 0.

The set of variables are defined as follows:

f i jk : the total traveled distance from depot to node j ∈ T ∪ {0} when traversing arc(i, j) by vehicle k.

i∈T ∪{0}

x i jk i∈T ∪{0}

x jik y jk ∀ j ∈ T , ∀k ∈ P,

j∈T ∪{0}

x j0k 1 ∀k ∈ P,

j∈T ∪{0}

x 0 jk 1 ∀k ∈ P,

k∈P j∈T

The objective function (1) maximizes the total number of covered customers. Constraint (2) is to limit the maximum travel time of each vehicle k ∈ P. Constraint (3) implies the connectivity of each route. Constraint sets (4) and (5) assure that all of the routes originate from and terminate at the central depot. Constraint (6) shows that customer i can only be assigned to facility j visited by vehicle k, if y jk 1. By Constraint (7), each customer can be assigned to at most one facility. Constraint (8) indicates that each facility can be visited at most by one vehicle. Constraint sets (9)-(13) prohibit sub-tours which are adopted from (Kara 2011) . Constraint (9) is the flow conservation constraint. By using this set of constraints, the total distance traveled from the central depot to node j ∈ T ∪ {0} minus the distance traveled from central depot to node i ∈ T , when traveling (i, j) by vehicle k, is equal to the distance from node i to node j. Constraints (10)-(15) are valid inequalities. In particular, Constraint (10) indicates that the total distance from node i ∈ T to the central depot, when traversing arc (0, i) by vehicle k ∈ P, is equal to the distance between them. Constraint (11) stipulates that the distance from the central depot to the facility i ∈ T which is traveled by vehicle k ∈ P cannot exceed L k . Constraint (12) sets upper bound on f i jk which guarantees the maximum distance traveled to j must not exceed L k minus the shortest distance that the vehicle k can travel from j to the depot. Constraint (13) sets lower bound on f i jk which shows that the distance traveled to j must be at least equal to the sum of the distance from the depot to i and the distance between i and j. Constraint (14) does not allow sub-tours with only two facilities. Constraint (15) represents that vehicle k must visit facility j when there is an arc entering to j. Finally, Constraints (16)-(19) specify the type of decision variables.

Due to the nature of the problem which makes it NP-hard and because of the applications in emergency situations, it is important to reach qualified solutions in a reasonable computing time. Exact solvers such as CPLEX could not reach optimal solution in a short time, especially for large size instances. Therefore, three heuristic methods, namely ILS, TS and VNS algorithms, are proposed for solving the introduced problem. In addition, we have developed different procedures which are used to generate initial solutions. In the following subsections, the proposed procedures for the initialization phase, the local search procedures which are used in the structure of the three metaheuristic algorithms and the properties of the ILS, TS and VNS are described, respectively.

In this section, we develop two procedures, namely "Sequential" and "Parallel" algorithms, to generate initial solution. In particular, the Sequential algorithm constructs one route at a time, while the Parallel one builds several routes at the same time (Solomon 1987; Potvin and Rousseau 1993; Pisinger and Ropke 2007) . In the following, we introduce the characteristics of these algorithms.

In this method, before constructing the routes, facilities are sorted in a descending order of their corresponding scores that will be defined later in this section. Starting from the central depot, the algorithm constructs the first route by visiting the first facility from the ordered set. Upon this step, the corresponding score of each facility is updated and the algorithm selects the next facility to be inserted into its best feasible position of the working route. Essentially, the best position is the place having the minimum extra insertion cost which is imposed to the solution by visiting the new facility. The algorithm keeps adding facilities to the route until it is not possible to visit a new facility because of the route's length threshold. After completing the first route, the procedure is continued by constructing a new route and the whole procedure stops whenever we have constructed p routes or covered the demand of all customers.

Associated scores for preference of facilities are calculated by using two methods. In the first method, score S j 1 of each facility j is equal to the number of uncovered customers within its coverage distance. In the second method, score S j 2 is computed by applying the following relation.

In this equation, d j indicates the distance imposed to the solution by inserting facility j into its best position on the tour and F is the set of unrouted facilities. It is worth mentioning that for the first facility, the value of d j corresponds to its distance from the depot.

In this section, we introduce two types of algorithms, namely "Parallel I" and "Parallel II." In Parallel I algorithm, facilities are sorted by utilizing one of the methods proposed in Sect. 3.1.1. At the first step, the algorithm selects the p first facilities from ordered set and p routes are constructed at the same time. At each step, the corresponding scores of the remaining facilities are updated and a facility is selected from the ordered set and assigned to the best feasible position generating the minimum extra insertion cost.

In Parallel II algorithm, at the first step, the algorithm selects the p facilities located as far as possible from each other. Following this step, p routes are generated by connecting each facility to the central depot. Routing of the remaining facilities is the same as Parallel I algorithm.

In this section, we provide the details of the local search procedures which are used within the structure of the heuristic algorithms. In particular, four procedures have been developed, for which the details are provided in the following.

Two types of swap move including intra-route and inter-route are used. The goal of the intra-route swap, called "Swap I," is to minimize the length of the each route k, i.e., i∈S k ∪{0} j∈S k ∪{0} t i j x i jk . In this definition, S k is the subset of facilities visited by vehicle k ∈ P. To accelerate the search process, swapping facility j ∈ S k with other facilities is considered only for those facilities located within a pre-determined neighborhood radius of j. Neighborhood radius (R) is achieved by applying Eq. (21).

In the above equation, θ is an input parameter taken from [0, 1] for which the corresponding value will be set in Sect. 4.4.

The goal of the second type of swap, called "Swap II," is to minimize the total length of the routes, i.e., k∈P i∈S k ∪{0} j∈S k ∪{0} t i j x i jk , by considering the replacement of each facility j ∈ S k with all of the other visited facilities within its neighborhood radius.

This move is seeking the maximization of the objective function value, i.e., i∈W j∈T k∈P d i j z i jk , by visiting a set of facilities from F which can cover at least one uncovered customer. For these facilities, the best insertion position is specified and the facility with the largest value of S j 2 (see, Eq. 20) is added to its best feasible location. By inserting the first facility, the set F and the value of the objective function are updated. This procedure continues until the possibility of inserting all members of F is examined.

Since some facilities cover common customers, by deleting a routed facility, all of its covered customers may be served by other routed facilities. So, by removing a routed facility, without any changes in the objective function value, i.e., i∈W j∈T k∈P d i j z i jk , the length of the routes may be decreased. Therefore, the possibility to insert other facilities may increase and the objective function will increase as well. This move is done for all the routed facilities, and the facility with the most reduction in the total length of the routes is selected.

This procedure includes two types. In the first one "Extraction-Insertion I," the goal is to decrease the tour's length by repositioning the routed facilities. In particular, a facility is extracted from its location and reinserted to the place on the same route with the minimum reinsertion cost. Applying this procedure is tested for all of the routed facilities, and a facility with the most reduction cost is selected.

In the second type of this procedure "Extraction-Insertion II," the goal is to increase the value of the objective function, or in the case that the objective value could not be improved anymore, it is to decrease the tours' length without reduction in the objective value. In particular, it is done by removing routed facilities from their locations and adding unrouted facilities into the best position on the available routes. Starting from the first constructed route, the first visited facility is extracted from its location and a set of unrouted facilities which covers at least one of the uncovered customers is inserted into the solution. At each step, the added facility is the one having the minimum insertion cost. If this substitution leads to an improvement in the objective function or it can decrease the total length without any reduction in the objective function, the corresponding move is saved. Applying this procedure is checked for the rest of the routed facilities. Finally, the substitution with the best improvement is selected for executing the procedure.

In this section, we describe the ILS procedure which is used to solve the proposed problem (Algorithm 1). The ILS is a metaheuristic algorithm that combines local search and perturbation procedures to search the feasible region more efficiently (Lourenço et al. 2003) . In particular, this algorithm has been successfully applied to solve similar problems such as TOP (Vansteenwegen et al. 2009; Gunawan et al. 2015; Zachariadis and Kiranoudis 2011) . For this reason, we use the ILS as one of the proposed methods for solving our problem. In this algorithm, a solution is created by following one of the proposed procedures in Sect. 3.1. The algorithm is composed of two loops. The inner loop is to improve the solution's cost by applying local search procedures. Essentially, at each iteration, a local search procedure is selected by using the roulette wheel technique (see Sect. 3.3.1). Before entering the outer loop, as long as the solution can be improved, the LS algorithm is applied. The pseudocode of the LS procedure is depicted in Algorithm 2. In fact, we tested several orders of the moves and the current order is the best combination. In this combination, at first we utilize the swap procedures followed by the extraction-insertion procedures. The reason of using the deletion and insertion procedures is that they can further improve the solutions obtained by applying the two mentioned procedures (i.e., swap and extraction-insertion), and this is why we call them after each execution of the swap and extraction-insertion procedures. The whole procedure is iterated for "iter1" iterations. In addition, during the execution of the outer loop, the goal is to escape from local optima by utilizing the perturbation procedure. The outer loop is repeated "iter0" iterations. The corresponding values for iter0 and iter1 will be set in Sect. 4.4.

In each iteration of the ILS algorithm, one of the local search procedures is selected based on roulette wheel (see line 11 in Algorithm 1). In this method, in the first iteration the possibility of choosing a move is equal to 1 n where n is the total number of moves. The possibility of selecting different procedures in later iterations is proportional to the effectiveness of the procedures in the previous iterations. Essentially, if the best solution is improved by applying a specific procedure, the possibility of selecting the corresponding procedure is increased in the subsequent iteration.

To escape from local optima, after every "iter1" iterations of the inner loop, the perturbation procedure is applied on the best found solution (see line 31 in Algorithm 1). In this procedure at first, ρ percent of the visited facilities are selected, randomly, and removed from their corresponding routes and added to the set F (set of unrouted facilities). Each time a facility is removed from the route, the assignment of customers is updated. In the second step, a facility is selected from F, randomly, and is added into its best feasible insertion position. The preliminary results show in case that we use the same set of customers in the structure of the perturbed solution, there will be a high probability to get stuck in a local optima solution by applying the local search procedure. In fact, by inserting the facilities randomly, we try to escape from local optima and expand the search space. After inserting the facility, the objective function and the subset F are updated. This step is iterated as long as we are not able to insert a new facility.

TS is a metaheuristic algorithm introduced by Glover (1989) . To escape from local optimum solution, this algorithm always moves to the best neighbor solution, even if the corresponding neighbor solution does not lead to an improvement in the cost of the current solution.

The proposed TS algorithm starts from the best solution obtained by utilizing the LS procedure (see Algorithm 2) and in each iteration selects a local search procedure by using the roulette wheel technique to explore the current solution's neighborhoods. For each selected neighborhood, the "best neighbor" with the most improvement in the cost is selected. If the implemented move is not tabu, it is selected as the new incumbent solution. Otherwise, the move is applied only if it has improved the best solution that has been found. In addition, we execute the diversification phase after "iter4" iterations without any improvements in the cost of the best solution. Finally, the intensification phase is applied after each improvement of the best solution which is discussed in Sect. 3.4.1. The pseudocode of the TS algorithm is provided in Algorithm 3 for which the details are given in the following subsections.

Upon applying each procedure, the attributes of the obtained solution must be considered as tabu in order to prevent choosing them for a specified number of iterations. In this paper, we have two types of tabu list as the short memory. In the first tabu list, the attributes of the swap I and the swap II are stored, while the second list contains the attributes of the Extraction-Insertion I and the Extraction-Insertion II. Attributes of the first tabu list are stored as (i,j) which indicates swapping facilities i and j are forbidden for a fixed number of iterations. Attributes of the second list are the indices of the facilities extracted from or inserted to one of the routes through the corresponding extraction-insertion procedure. So, if the facility has inserted (extracted), it will not be extracted (inserted) for the following τ iterations.

Long-term memory is based on the frequency of edges which appear in different solutions. During search process, whenever we have an improvement in the cost of the best solution, the corresponding edges of the improved solution are kept in longterm memory to be used in the diversification phase. In particular, diversification is implemented after "iter4" iterations without any improvements. In this phase, the objective function of solution s, i.e., f (s), is decreased by imposing a penalty which is proportional to the frequency of edges in s. To do so, a new objective function is calculated based on Eq. (22). In this equation, f (s) is the number of covered customers by visited facilities, i.e., i∈W j∈T k∈P d i j z i jk .

In addition, r i j is the number of improved solutions that contain edge (i, j) and σ denotes the total number of improved solutions. In the diversification phase, selection of the local search procedures is based on the roulette wheel technique and evaluation criterion for selecting neighborhoods is g(s) instead of f (s). In particular, by using g(s), solutions consisting of edges with less frequencies during the search history will have a lower penalty. Diversification is iterated for "iter5" iterations, and in each iteration, the neighborhood which has the best value of g(s) is selected. Finally, in each iteration of the TS, if the selected move improves the best solution, local search phase is applied as the intensification strategy.

Proposed by Mladenović and Hansen (1997) , the VNS is a metaheuristic algorithm that explicitly applies a strategy based on dynamically changing the neighborhood structures. The general framework of the VNS algorithm for the introduced problem is given in Algorithm 4. The algorithm starts to construct an initial solution (Cur-rentSolution) by applying the initialization procedure (see Sect. 3.1). It then improves upon this initial solution by applying the LS procedure. The improvement of the solution continues in a loop until the termination criterion is met. The loop contains the perturbation procedure and the LS procedures. If a solution is improved by applying the perturbation and local search procedures it will be accepted as the incumbent solution. At each iteration of the VNS algorithm, we apply the perturbation procedure in order to dynamically expand the neighborhood structure. In particular, the algorithm generates a solution which is in the neighborhood size ρ of the CurrentSolution, i.e., Perturbation (Current Solution, ρ) . In this definition, ρ is the percentage of customers which is extracted by utilizing the perturbation procedure (see Sect. 3.3.2).

Several instances have been generated to evaluate the effectiveness of the proposed algorithms. These samples are randomly generated and classified into three sets including small, medium and large size data. In particular, we use the instances from the TSP library having 52, 76, 100, 150, 200, 318, 417, 575, 657 and 724 nodes (Reinelt 1991) . In generated instances each facility can cover up to 5 nearest customers, randomly. By using each instance, we have generated 3 new instances in which for each sample, 50, 60 and 70% of the nodes are considered as customers (|W |). Three different scenarios, namely 2, 3 and 4, are considered for the number of vehicles (|P|). Finally, three values are proposed for L k as follows:

In this equation α ∈ {1, 0.9, 0.8}. As a result, totally we have 10 × 3 × 3 × 3 270 instances. The instances are divided into the three groups. In particular, the groups with small and medium instances contain instances having up to 100 and 200 nodes, respectively, while the other instances belong to the group with large instances.

ILOG CPLEX has been used to run the mathematical model in Sect. 2. In particular, for each instance 5 h of computing time is put on the maximum run-time of CPLEX and the results are reported in Tables 1, 2 and 3. In these tables, number of nodes, number of facilities, number of customers, number of vehicles and the maximum travel length are provided in columns 2-6, respectively. The corresponding objective function is provided in the column "Obj." and the corresponding run-time of CPLEX (in seconds) is given in the column labeled by "Time". Finally, the optimality gap (in percentage) is provided in the column "Gap". Totally, 78 small instances (out of 81), 34 medium instances (out of 54) and 32 large instances (out of 135) have been solved to optimality within the given time limit. For those instances that CPLEX is not able to reach the optimal solution, the average optimality gap is 2.28% for the small instances, while this is 3.84% for the medium instances. In addition, the average computing time for the small and medium size instances is 905.75 and 7803.84 s, respectively. As it is shown in Table 3 , CPLEX is able to reach the feasible solution only in 38 instances (out of 135) of the large size data and it fails to reach a feasible solution for the other 97 instances. The average computing time for 38 instances for which CPLEX can achieve a feasible solution is 1079.38 s. There are some instances for which the optimality gap is a positive value, while CPLEX has not reached the time limit. In these instances, the computer runs out of memory before finding the optimal solution and the best found solution is reported in the corresponding position in Table 3 .

According to the applications of this problem in critical situations and emergencies and due to the fact that in real-world situations we may face the larger size problems, we should use heuristic methods to accelerate obtaining high-quality solutions in a reasonable computing time.

As discussed in Sect. 3.1, we have proposed three methods, namely Sequential, Parallel I and Parallel II, for generating initial solutions. In addition, we developed two scoring methods, i.e., S j 1 and S j 2 , for preference of facilities. Therefore, we have 6 different scenarios to generate initial solutions.

Since CPLEX fails to reach a feasible solution for several large instances, we have only used the results obtained by CPLEX for the small and medium size data in order to analyze the performance of the initialization procedures. Results by utilizing each of the 6 scenarios on small and medium instances are given in Table 4 . In particular, for each scenario the overall average gap for available instances, with respect to the results obtained by CPLEX, is given in the "Avg. gap" row. Based on the results, the Parallel I algorithm with S j 1 scoring method has the minimum average gap, while the worst gap is obtained by Parallel II algorithm with S j 2 scoring method. Finally, the row labeled by "Avg. time" shows, on average, different methods need almost the same computing time to achieve an initial solution.

The proposed heuristic algorithms contain a set of parameters to be tuned. For each algorithm, Table 5 gives the assigned values for the involved parameters and the final selected values. To set the values, we have chosen one-third of the small and medium size problems as the set of benchmark instances to run the experiments. In particular, all the combinations of the parameters are tested and the best combination leading to the best average performance is selected for the corresponding algorithm. To analyze the results, we have used Wilcoxon nonparametric test with 95% of confidence. This test is presented in Wilcoxon (1945) and used instead of t test when the probability distribution of samples is unknown. For each algorithm, to select the best combination of parameters we select the combinations having the best and worst overall performance for the set of benchmark instances. Wilcoxon test shows that the difference between the results of these two combinations is not statistically significant. As a result, for each algorithm the combination having the minimum computing time is chosen to run the final experiments.

Each of the initialization procedures is used to run the ILS algorithm. Essentially, for each of the initial solutions we run the ILS algorithm five times with five different random seeds and the average and the best performance of the ILS algorithm are provided in Table 6 . In this table, for each initialization procedure, the "Avg. gap" row gives the overall average gap of the CPLEX results with the solutions obtained by utilizing the ILS algorithm on all the small and medium instances over the five runs. In addition, the best gap over five different runs is given in the row labeled by "Best gap". Finally, the average computing time for all the small and medium instances over the five different random seeds is given in the row represented by "Avg. time". According to the results, on average, all combinations have better performance with respect to CPLEX. Considering the best performance of the algorithm, best methods for generating initial solution are Sequential algorithm by using the S j 2 as the scoring method, and Parallel I and Parallel II algorithms by using S j 1 as the scoring method. Finally, for each row of the table, the column labeled by "Total" gives the average of values over the corresponding row.

For small and medium instances, the results of applying the VNS algorithm are given in Table 7 . In particular, each row of this table reports the same information as that already described for Table 6 . Taking into account the average performance of the algorithm over five runs, Parallel II algorithm and Sequential algorithm with S j 1 as the scoring method have the best performance. In addition, Sequential method with S j 1 scoring method, Parallel I with S j 1 scoring method and Parallel II have the best performance, when considering the best performance of the algorithm. Based on the results reported in the last column (i.e., Total), on average, the ILS algorithm outperforms the VNS algorithm. 

In this section, we analyze the results obtained by the TS algorithm for small and medium instances. According to the results reported in Table 8 and taking into account the average performance of the algorithm over 5 independent runs, Parallel II algorithm has the best performance. In addition, Sequential method with S j 2 as the scoring method and Parallel II are the best scenarios when considering the best performance of the algorithm. Finally, comparing the results obtained by all three heuristic algorithms indicates the superiority of the TS algorithm for the small and medium instances. In particular, when considering the average performance of the algorithms, the overall average gap is − 0.49 for the TS algorithms, while this is − 0.29 and − 0.46 for the VNS and ILS algorithms, respectively.

In this section we provide the results of the CPLEX and heuristic algorithms over the set of large instances. In particular, according to the results reported in Tables 6,  7 and 8, for each of the heuristic algorithms the best scenario is chosen to run the experiments.

Each heuristic algorithm is run five times, and the average results are reported. The results show that for 38 instances in which CPLEX is able to reach a feasible solution, the proposed TS algorithm can improve the result by 2.64%, while this is 1.66% and 1.94% for the VNS and ILS algorithms, respectively. In addition, there are 133 instances for which TS has obtained the best solution, while VNS and ILS have achieved the best performance in 109 and 112 cases, respectively. Detailed results of the heuristic algorithms over the large size instances are provided in Appendix 1.

As discussed in the definition of the problem, the TCMCSP is a special variant of our proposed problem. So, we have adopted our TS algorithm to be applied for this problem. For generating initial solutions, we have used Parallel II algorithm with S j 1 scoring method. The results are reported in Table 9 , in which columns labeled by |W | and L k have the same meaning as those already introduced in Table 1 . Column labeled by "Best Obj. E " gives the best results achieved by the exact algorithm, proposed in Naji-Azimi and Salari (2014). For each instance "Avg.Obj. H " and "Best Obj. H " give, respectively, the best and the average performance of the proposed heuristic algorithm in Naji-Azimi and Salari (2014) over 5 different runs. Columns "Avg.Obj. T " and "Best Obj. T " give the average and the best results of the TS algorithm over five different runs, respectively. "Best Gap T E " column represents the gap between the best performance of TS algorithm (i.e., column Best Obj. T ) and the results reported in column "Best Obj. E ". Finally, the column "Best Gap T H " shows the gap between the best performance of our TS algorithm and the best performance of the heuristic algorithm proposed in Naji-Azimi and Salari (2014). As it is given in the last row of the table, on average, the TS algorithm improves the best known heuristic results in Naji-Azimi and Salari (2014) by 0.15%.

Taking into account the best performance of the TS algorithm, the results show that this algorithm is able to improve the solutions obtained by CPLEX by 1.27%. Average running time for the TS is 1.66 s which is a great improvement in comparison with 1680.24 s in CPLEX. We should mention that the average time of the proposed heuristic algorithm in Naji-Azimi and Salari (2014) is 1.54 s.

As already reported in Sect. 1, the TCMCRP is a generalization of the TOP. For this reason, in order to validate our proposed algorithm, we have modified the TS algorithm to be used for the TOP instances. In particular, we use 212 benchmark instances developed by Chao et al., which are divided into 4 datasets (Chao et al. 1996) . In each dataset, the position and the profits of customers are the same for all instances, while the number of vehicles varies between 2 and 4. In addition, the travel limit is different for available instances.

We compare the performance of the TS with three heuristic algorithms, namely MA, PSOiA and PSOMA, and the results are reported in Table 10 . In particular, the MA is a memetic algorithm proposed in (Bouly et al. 2010) . PSOiA is a heuristic algorithm that uses the particle swarm optimization (PSO) to solve the problem (Dang et al. 2011) . Finally, the PSOMA is a PSO-based memetic algorithm that combines the PSO and the memetic algorithms (Dang et al. 2013 ).

For each algorithm tested on an instance, the algorithm is run 10 times. Following this step, for each instance the best obtained solution (Z best ) by utilizing 4 different algorithms is used to analyze the effectiveness of the proposed TS algorithm and the results are reported in Table 10 . In this table, for each instance set, the average gap 

Comparing the algorithms for the TOP instances and the best gap of different heuristics with respect to Z best are reported in rows "Avg. gap" and "Best gap", respectively.

The results show the efficiency of the TS algorithm. In particular, taking into account the average performance of the TS algorithm the worst results have been achieved in set p.4 for which the overall average gap with respect to the best solution is 0.57%. Considering the best performance of the TS algorithm, the results show that in the worst case the average gap is 0.27% for the set p.4. Detailed results of available algorithms are provided in Appendix 2.

In order to evaluate the effectiveness of the local search procedures which have been used in the structure of the proposed algorithms, we calculate the percentage of improvement made by each procedure as follows. In this regard, we consider the proportion of the iterations in which a specific procedure can improve the current solution, to the total number of the iterations that the procedure had been called. The obtained results are summarized in Table 11 . As it can be concluded from this table, all the presented procedures are useful and the most effective procedure is Extraction-Insertion I, which has been successful in 37.50% of its execution.

In this section, we do some tests to analyze the impact of the parameters on the behavior of the TS algorithm. To do so, we fix all parameters and change one parameter at a time. Essentially, ρ and iter2 parameters have been selected. ρ is the percent of the facilities which has been removed from the solution in perturbation phase parameter (see Sect. 3.3.2). In our experiments ρ ∈ {0.6, 0.7, 0.8, 0.9} and the impact of this parameter for the TCMCSP is given in Fig. 2 . In this figure, the vertical column shows the overall gap between the TS and the best solutions proposed for the TCMCSP instances in Naji-Azimi and Salari (2014). In addition, Fig. 3 gives the same information for the iter2 parameter. We have selected five different scenarios, namely 90, 100, 110, 120 and 130, for running this test. Results show that the best result has been achieved when iter2 is equal to 120. Figures 4 and 5 show the impact of iter2 and ρ parameters over the 135 instances of small and medium size for TCMCRP. Vertical axis shows the improvement in CPLEX' results by applying the TS algorithm. The iter2 and ρ are set to 120 and 0.8, respectively. 

In this section, we investigate the relation between L k and the total percentage of covered customers. In particular, by increasing the corresponding value of α which concludes increasing L k parameter, the percentage of covered customers increases as well. According to Fig. 6 , the decision maker is able to select the best value for α by considering the number of customers whom he wants to cover. Figure 7 shows the relation between the number of customers that can be covered by each facility and the percentage of customers that are covered by routed facilities. According to this figure, decision maker can determine the coverage radius of each facility by considering the distance between facilities and customers.

In this paper, the time-constrained maximal covering routing problem (TCMCRP) has been introduced, which is the combination of the covering salesman problem and the team orienteering problem. In this problem, a central depot, a set of facilities and a set of customers are given. The TCMCRP objective is to maximize the number of customers covered by facilities which are visited by a set of length constraint vehicles. We proposed a flow-based mathematical model to solve the model. Since the TCMCRP is a NP-hard problem, exact model is not effective for real size instances and heuristic algorithms including ILS, VNS and TS are presented for solving the problem. We tested the performance of these algorithms on 270 instances having 52-724 nodes and dividing into small, medium and large size sets. The numerical results showed that on average all the proposed algorithms outperform CPLEX especially in large size instances. Regarding the performance of the heuristic algorithms, the tabu search algorithm has the best performance in comparison with the other two algorithms. In addition, to evaluate the effectiveness of the developed heuristic algorithms we have adapted our algorithms to be used for the time-constrained maximal covering salesman problem and the team orienteering problem. The results indicate that developed heuristic algorithms can achieve promising results for the similar problems. In particular, by applying the TS algorithm on the TCMCSP instances, the results show the superiority of the TS with respect to the best results available in the literature. The results also indicate that the TS algorithm performs efficiently when applying on the TOP benchmark instances.

For the future research, we propose to develop a variation of the problem in which we are given a threshold on the minimum percentage of customers that have to be covered while minimizing the total travel time. This problem happens in a situation that we have a limited amount of resources to be delivered to the customers within the shortest possible time.

See Tables 12 .  Table 12 Performance of the heuristic algorithms for large instances of the TCMCRP Instance Obj. Table 12 continued Instance Obj. 103.00 Table 12 continued Instance Obj. Table 12 continued Instance Obj. 

See Tables 13, 14 , 15, 16. 

Selective multi-depot vehicle routing problem with pricing

The undirected capacitated arc routing problem with profits

Multi-depot vehicle routing problem with time windows considering delivery and installation vehicles

The capacitated m-ring-star problem

The capacitated orienteering problem

A memetic algorithm for the team orienteering problem

The waste collection vehicle routing problem with time windows in a city logistics context

A heuristic for the multiple tour maximum collection problem

The team orienteering problem

The maximal covering location problem

The covering salesman problem

The median tour and maximal covering tour problems: formulations and heuristics

A pso-based memetic algorithm for the team orienteering problem

An effective PSO-inspired algorithm for the team orienteering problem

The truck dispatching problem

A branch-and-price approach to the vehicle routing problem with simultaneous distribution and collection

Solving the capacitated vehicle routing problem with environmental criteria based on real estimations in road transportation: a case study

The covering tour problem

Tabu search-part I

The vehicle routing problem: latest advances and new challenges

An iterated local search algorithm for solving the orienteering problem with time windows

Orienteering problem: a survey of recent variants, solution approaches and applications

The electric fleet size and mix vehicle routing problem with time windows and recharging stations

A new efficient approach for solving the capacitated vehicle routing problem using the gravitational emulation local search algorithm

Self-imposed time windows in vehicle routing problems

Arc based integer programming formulations for the distance constrained vehicle routing problem

Vehicle routing problem with simultaneous pickup and delivery: mixed integer programming formulations and comparative analyses

Metaheuristic algorithms for solving two interconnected vehicle routing problems in a hospital complex

The ring star problem: polyhedral analysis and exact algorithm

An improved formulation for the multidepot open vehicle routing problem

The selective travelling salesman problem

A branch and cut algorithm for a Steiner tree-star problem

Iterated local search

Variable neighborhood search

A tabu search algorithm for the vehicle routing problem with simultaneous pick-up and delivery service

The time constrained maximal covering salesman problem

The maximal covering location problem with capacities on total workload

A general heuristic for vehicle routing problems

A parallel route building algorithm for the vehicle routing and scheduling problem with time windows

A traveling salesman problem library

Algorithms for the vehicle routing and scheduling problems with time window constraints

Statistical analysis of distance-based path relinking for the capacitated vehicle routing problem

The multiconstraint team orienteering problem with multiple time windows

Multiple depot ring star problem: a polyhedral study and an exact algorithm

Philadelphia Tsiligirides T (1984) Heuristic methods applied to orienteering

Iterated local search for the team orienteering problem with time windows

The orienteering problem: a survey

An integrated multi-depot hub-location vehicle routing model for network planning of parcel service

Individual comparisons by ranking methods

Local search for the undirected capacitated arc routing problem with profits