A GRASP-Tabu Heuristic Approach to Territory Design for Pickup and Delivery Operations for Large-Scale Instances

1Faculty of Engineering and Applied Sciences, Universidad de Los Andes Chile, Monseñor Álvaro Portillo 12455, Las Condes, Santiago, Chile 2Tecnologico de Monterrey, Escuela de Ingenieŕıa y Ciencias, Eugenio Garza Sada 2501, 64849 Monterrey, NL, Mexico 3School of Computing, Informatics and Decision Systems Engineering, Arizona State University, Tempe, AZ 85287-8809, USA 4Universidad Autónoma de Nuevo León, Facultad de Ciencias Fı́sico-Matemáticas, Av. Universidad s/n, 66450 San Nicolás de los Garza, NL, Mexico


Introduction
In this study, we consider a logistic districting problem faced by a parcel company that operates in the metropolitan area of Monterrey, Mexico.This company is concerned with the delivery and pickup of packages within Monterrey.The latter is divided into districts for logistic purposes.Each district is served by a single vehicle that departs from a central depot, in which the packages are received.At the end of the day, the vehicle comes back to the depot.Currently, drivers take the route decisions in a dynamic manner during the day due to the variability of the demand.That is, it is assumed that the driver may receive new dispatch orders along the day.
The company's current practice is to redesign the districting configuration every year and a half in order to account for demand variations.Although the day to day locations of the customers and the daily volume of demand are stochastic, the logistics planning managers consider that a reasonable approach is to use the data of a high workload day which is representative of the current and growing demand.The number of pickups and deliveries may vary, but, at the time, the company was considering data sets of about 1,100 customers, with about equal numbers of pickups and deliveries.In general, the days chosen by the company to serve as representative days have above-average demand and a geographical distribution of demand that is considered to be "typical," that is, not unusually high or low in any area with respect to an average behavior.
Once a representative day is selected, the logistics manager identifies the locations of the customers on a map.Then, the districts are readjusted considering their estimated capacity.But this readjustment needs to consider the compact Routes Problem (FRP) studied in [10], in which the service region is divided into subareas but each route does not vary from day to day.
Daganzo [11] proposes an approximation method for the design of multiple-vehicle delivery zones seeking tours of minimum total length.The objective in that paper is to explore the impact that the zone shape has on the expected length of each route.Later, in [12], a methodology in which the region is partitioned into zones of nearly rectangular shape elongated toward the source is presented.The number of points considered in the instances is large compared to the capacity of the vehicles; this fact complicates the problem's resolution.Newell and Daganzo [13] analyze the districting of a region in which the underlying network of roads forms a dense ring-radial network.The authors proposed an approximation method for the design of multiple-vehicle delivery tours that minimizes the total traveled distance.The design of delivery zones for distributing perishable freight without transshipment is studied in [14].Lei et al. [15] propose the multiple traveling salesperson and districting problem, which considers multiple periods and depots.The authors propose an adaptive large neighborhood search metaheuristic.
Regarding problems for designing delivery tours within districts, it can be mentioned [16], in which the design of multiple-vehicle delivery tours that satisfy time constraints for letter and parcel pickup is studied.The authors propose a methodology that involves partitioning the region into approximately rectangular delivery zones that are arranged into concentric rings around the depot using a continuous approximation of the model.Rosenfield et al. [17] study the problem of planning service districts with a time constraint.One of their main contributions is the derivation of analytical expressions to determine the optimal number of service districts for the US postal system.A methodology to design multidelivery tours associated with the servicing of an urban region of irregular shape is presented in [15].The authors' methodology is based on a sweep approach and assumes a rectangular grid structure.Novaes et al. [18] present a methodology for solving the same problem as [19], but they used a continuous approach to represent the region.Previous works are extended by [20] introducing some improvements to the ring-radial model.The authors present a special case of a Voronoi diagram and model the problem with a continuous approximation.
Muyldermans et al. [21,22] address the problem of districting for salt spreading operations on roads.They assume that each district is served by a single facility aiming to minimize the deadheading distances and the number of vehicles required to service the region.The authors present a novel model based on a graph with demand at the nodes.Another districting problem modelled on a graph but considering the demand in the arcs is presented in [23], in which the concept of Eulerian districts is introduced.
Haugland et al. [24] consider the problem of designing districts for vehicle routing problems with stochastic demands.They presented a Tabu Search (TS) and multistart heuristics to solve the problem.Also, a multiobjective dynamic stochastic districting and routing problem is considered in [15].In this problem, the customers of a territory stochastically evolve over the planning periods.Therefore, several objectives must be optimized via a coevolutionary algorithm.Another nondeterministic districting problem is studied in Sheu [25], in which an integrated fuzzy-optimization customer grouping based logistics distribution methodology is presented.The proposed method groups customers' orders primarily based on the multiple attributes of customer demands.Then, the customer group-based delivery service priority is determined, and finally the routes for deliveries are established.
Concerning multiple objectives to be considered in a districting model, Tavares-Pereira et al. [26] consider a districting problem with multiple criteria and propose an evolutionary algorithm with local search to approximate the Pareto front.Also, a problem related to a bottled beverage distribution company is studied in [27].Hence, a biobjective programming model is introduced considering dispersion and balancing with respect to the number of customers in each district as performance criteria.Two continuous location-districting models applied to transportation and logistics problems are developed in [28].A Voronoi diagram is combined with an optimization algorithm for solving both proposed models.
In [29,30], the same districting design problem for the operations of a parcel company is studied.In both references, the same mathematical model is proposed but in [30] a twostage approach is considered.The second stage is related to a mathematical model that considers the solution of the first stage and, then, aims to minimize the dispersion of the workload in the districts.Small and medium size instances are solved via a hybrid algorithm.That algorithm is unable to efficiently solve large size instances.In [31], there is a review of the state of the art regarding modelling techniques and solution methods for problems in supply chain planning that handle districting or customer clustering.
The problem herein addressed is closely related to [29,30].However, there are significant differences among them.Particularly, in this paper we are proposing an adjusted heuristic that is capable of solving large-scale instances with low computational effort.In previous papers, only limited size instances are solved, but herein the proposed heuristic efficiently solves larger instances.The main differences between the algorithms is that the procedure for constructing feasible solutions, the neighborhoods considered, and the movements explored during the local search are modified.The latter changes allow the heuristic to solve large size instances.Moreover, instances with different structure are tested, and the heuristic algorithm shows to be robust.Specifically, a set of symmetric instances is considered in which the optimal solution is known due to their particular structure.These instances allow properly measuring the efficiency of the proposed heuristic.
Furthermore, other important characteristics of this research are that we aim to optimize the balance of workload content and compactness in a single objective function instead of a single objective (as usually).Also, parcel applications are scarce in the literature.Hence, this research is worthy due to parcel companies handling many customers and the proposed heuristic can efficiently solve these situations.
The modeling structure proposed in this work also differs from others found in the literature, mainly in the fact that the districting problem is modeled as a graph in which demand occurs at the nodes.
In our situation under study, districting is a strategic decision that is not defined day to day and demand points are not fixed.We also distinguish between pickup and delivery operations.The only reviewed paper to do so is [16], but the objective is to minimize the expected length of the tours and the workload balance is addressed as a constraint.So, it differs from the model presented in this paper, in which we aim to balance workload content and also achieve districts of compact shape.

Mathematical Formulation
Consider a connected, undirected graph (, ), where  is the set of vertices and  is the set of edges.In general, the graph is not complete.All the edges   = (V  , V  ),   ∈ , (V  , V  ) ∈  × , have a positive length and represent a real path between adjacent vertices V  and V  .Distances between vertices are edge lengths for adjacent pairs and shortest path distances for nonadjacent pairs.A district is defined as a subset of the vertices and  = {1, . . ., } is the district set.Each vertex represents a customer that may require either a pickup or a delivery.The depot is defined as the vertex {0}.
We define   and   as the number of pickups/ deliveries requested by each demand point ,  ∈ ; and   and   as the stopping time per pickups/delivery in each demand point ,  ∈ .As mentioned previously, the distance between a pair of demand points is denoted as   . represents the average vehicle speed. and  are continuous variables that represent the maximum workload content and compactness metrics, respectively.  is a continuous auxiliary variable that takes the value of the travel time from the depot to the farthest point in district ,  ∈ .
The districting procedure seeks to optimize two criteria simultaneously: balancing workload among the districts and compactness of their shapes.Compactness is not defined uniquely for all the districting problems in the literature and it is generally defined according to the application context.For this problem, we define it as the distance between the two furthest apart vertices in a district, that is, max{  }, considering all vertices V  and V  belonging to a district.We employ a minimax objective in which the maximum compactness metric among the districts is minimized.Although the compactness metric that we employ could have the same value for districts of different shape (e.g., a long thin district), when combined with the objective to balance workload content, in practice it produces districts of approximately circular or square shape.
For each district, the workload content is defined as the time required to perform all required pickups and deliveries and also includes an estimation of the time required to travel from the depot to a vertex in the district.The latter is approximated by the time required to travel from the depot to the farthest point in the district.Despite the fact that the closest vertex in the district or the centroid could also be employed to approximate the line haul distance from the depot, considering that the vehicle will visit all vertices in a district during the route, the farthest vertex is a reasonable metric.Additionally, the use of the farthest vertex simplifies the mathematical model.In order to balance the workload content among districts, the maximum workload among all districts is minimized.The following expression shows how the workload of a district is calculated: We propose a linear single objective mixed integer model in which a weighted sum of the compactness metric and the maximum workload content assigned to a district is minimized.Each metric is normalized with respect to an estimation of the optimal values of compactness and workload content ( and , resp.).Since a feasible number of vehicles considered in the model are known a priori, the capacity of the vehicles is not explicitly considered in the model.To ensure that a feasible solution with respect to vehicle capacity is found, we balance district workloads and we keep the number of pickups and deliveries within certain limits.These limits will not allow pickups and deliveries to be intermixed, and this will release vehicle capacity by deliveries to allow subsequent pickups.The limits on the number of pickups and deliveries for each district are denoted by  and , respectively.
The following decision variables are defined: The mathematical model is as follows: ∈ {1, 0} ; Equation ( 3) is the objective function that minimizes a weighted average of the maximum workload and compactness.Both terms in (3) are normalized and the relative weighting is given by , 0 ≤  ≤ 1. Normalization parameters are estimated based on their optimal solution values as follows:  is determined by estimating an average workload per district when it is approximately balanced by dividing the total workload of the region. is estimated as the traveling time between the two furthest points in a district, assuming that the service region is equally divided into  districts of equal shape (assuming a circular shape for normalization parameters estimation).Constraints (4) guarantee that each demand point is assigned to exactly one district.Constraints ( 5) and ( 6) ensure that each district has a maximum number of  pickups and  deliveries, respectively.Constraints (7) guarantee that   takes the value of the time from the depot to the farthest vertex of each district .Constraints (8) require that  take the value of the maximum travel time between the vertices assigned to a district.Constraints (9) guarantee that  takes the value of the maximum workload from among the districts.Constraints (10) are the binary and nonnegativity requirements.

Solution Methodology
The redistricting problem has been shown to be NP-Complete by Altman (1997).Moreover, since constraints ( 5) and ( 6) impose a limit in the amount of pickups and deliveries assigned to a district, it turns out that even finding a feasible solution is a difficult task.Motivated by the difficulty of the problem, we propose a multistart heuristic that hybridizes GRASP and Tabu Search.The procedure seeks a feasible solution (partition of the customers into districts) of the model proposed in Section 3 aiming to obtain good partition quality in terms of the values of the objective function in which the weighted sum of compactness and work balance metrics are optimized.We will further describe the steps of the algorithm.The general procedure consists of two phases: construction of a feasible initial partition and improvement by a local search.These phases are named as FIP-Construction and Local-Search, respectively.In the first phase, districts are conformed by partitioning all the customers into the districts in set . Hence, a solution of the proposed hybrid algorithm corresponds to a districting decision.The Local-Search procedure incorporates a Tabu Search (TS) short term memory, in which the recently visited partitions are labeled as Tabu active during a predefined number of iterations.The aspiration criterion allows a Tabu active move only if the resulting partition is better than the current best one.Among all the partitions created and improved, the overall best one is reported as the final output of the algorithm.
We propose two local search procedures that may be applied independently or in two different combinations.This results in four strategies that attempt to improve the current partition (set of districts) constructed during each iteration of the algorithm.A key concept is the adjacency of vertices and districts.A vertex is considered to be adjacent to a district if there exists at least one edge connecting the vertex with another one that already belongs to the district.Then, we restrict the allocation of the vertices only to its adjacent districts.Knowledge of the adjacency helps to avoid unnecessary evaluations that may result in long computational times in the local search and also to enhance compactness of the constructed partition.
The encoding of a partition  is as follows: there is a list associated with each district, in which the index of the customers associated with it are included.Moreover, the quality () of each partition  will be measured as the weighted sum of both objectives considered in (3).Algorithm 1 shows the pseudocode of the procedure in which  * represents the best value of the objective function given by (3) and  * is the corresponding best partition found.The term seed will be used hereafter to denote a vertex chosen to initiate the construction of a district.NSeeds is the number of seed selection procedures used in the FIP-Construction method.We employ five different procedures to select a set of seeds.The procedures are used in sequential order,  = 1, . . ., , during the iterations of the hybrid algorithm.In each iteration a feasible partition is attempted to be constructed and, if found, improved.

Phase I: Feasible Initial Partition Construction (FIP-Construction).
We apply a multistart approach, in which a determined number of initial solutions are attempted to be constructed based on a greedy randomized approach.For the construction of  districts, the same number of seeds would be chosen.Algorithm 2 shows the pseudocode of the construction procedure, for which the two main steps are selection of a set of  seeds (Seeds Selection) which employs a single seed selection procedure in each iteration () and allocation of vertices to the districts associated with each seed (Vertices Allocation).Set Ψ contains the unassigned vertices.To enhance compactness, vertices are attempted to be assigned to the closest seed as long as adjacency conditions are fulfilled.If an initial feasible partition has been constructed, a local search procedure is applied.Otherwise the partition is discarded.

Seeds Selection.
Five different approaches are used in this step ( = 1, . . ., ).Each time the algorithm performs an iteration to construct and improve a feasible partition, the five approaches are sequentially employed.Algorithm 3 shows the pseudocode of the general procedure, in which  is the set of potential vertices from which seeds,   ∈ , are selected.() is the greedy function employed in each procedure.The greedy and semirandom functions are further described.

P-Dispersion Greedy Function.
The objective is to select the seeds that are most dispersed relative to each other.The greedy function computes the sum of the distances between the already selected seeds and a potential vertex.Vertices associated with larger sums are considered to be more desirable.
Neighborhood Greedy Function.Vertices are evaluated according to the number of vertices that are located within a neighborhood defined by a threshold distance.The vertices with more neighbors are considered to be better.Once a seed is selected, the neighbors of the seed (adjacent vertices) are discarded as potential candidates.
Angle Greedy Function.Based on the location of the vertices.The service region is partitioned into  wedge-shaped radial sectors of equal size.The RCL is formed by the vertices whose locations are closest to the boundary between two sectors.Semirandom Selection.A vertex is randomly selected as a seed at each iteration and the adjacent vertices are discarded as potential seeds.If all the points have been eliminated before all seeds are selected, the remaining seeds are selected randomly.

Vertices Allocation.
Once the set of seeds has been selected, districts are formed around the seeds.This procedure not only attempts to create a feasible partition but also strives for compactness by assigning vertices to their closest seed subject to satisfying adjacency requirements.Vertices are allocated to districts in three consecutive steps.In each step, we attempt to construct a feasible partition, and the next step is applied only if a feasible partition was not found in the previous step.
The procedure is as follows: (1) a seed,  * , is randomly selected and the  closest vertices to the seed are inspected with the aim of assigning the closest adjacent vertex that satisfies the capacity of the district. is a control parameter and its value is determined based on preliminary experimentation (see Section 5.2).The procedure is repeated until a stopping condition is met.(2) The remaining unassigned vertices are attempted to be assigned to an adjacent district, respecting capacity, and each vertex is explored only once.
(3) The remaining unassigned vertices are attempted to be assigned to an adjacent districting, even if capacity is violated, but always aiming to assign the vertex to the district that least increases its number of pickups or deliveries beyond its capacity limits.If an infeasible assignment of vertices is generated, a last procedure is applied based on local search (exchange of vertices between adjacent districts).The aim of the latter is to find a feasible partition; in case that this is not possible, the partition is discarded.

Phase II: Local Search Phase (LS).
After a feasible partition has been constructed, its quality can be measured by computing the corresponding objective function value according to (3).Hence, the improvement of the current partition is sought.The LS phase contains four neighborhood structures, two of them are referred to as "1-Step LS" (1-S) and "K-Steps/Pairs LS" (K-S/P).Two additional neighborhood structures result from a combination of the above-mentioned neighborhood structures: "Hyperheuristic LS" (HypLS) and "2-Iterations LS" (2-IterLS).In an iteration of the LS phase, the evaluation of the objective function of all the partitions over a search space determined by one of the four neighborhood structures proposed is performed.In case of ties, we select the partition that presents the least dispersion of the workload content among the districts.
Each neighborhood structure implements a Tabu Search short term memory: recently visited solutions are labeled as Tabu active during  iterations.Tabu permanence  was defined as a static value.The aspiration criterion allows a Tabu active move only if the resulting partition is better than the current best one.Local search procedure is applied until stopping conditions are met.
For this phase, to simplify the exploration, a partition  is codified as the matrix   , with   = 1 if vertex  is assigned to district  and   = 0 otherwise.A movement in partition , consists in changing the vertex  adjacent to district   , that is, from   = 0 and    = 1, for all  and  ̸ =   ∈ .An exchange neighborhood of partition  consists of the set of partitions that, starting with   = 1 and    = 1, with point  adjacent to district   and point   adjacent to district , results from setting   = 0,    = 1,    = 0, and   = 1 for all ,  ∈ and  ̸ =   ∈ .The neighborhood structures are the following.

1-Step LS (1-S), 𝑁1(𝑥).
Each iteration is concluded by performing the best move found and setting the partition with   = 1 as Tabu active during  iterations.A list of the three best partitions is also maintained so that a final attempt is made to improve the three best ones found in hopes of finding a better partition.
-Steps/Pair LS (k-S/P), 2().This neighborhood is an extension of the algorithm proposed by Kernighan and Lin [32] for graph partitioning.It consists of the set of partitions that results during a number of steps,  = 1, . . ., .In each step , all moves and exchanges are explored for a pair of adjacent districts which are randomly selected.Partitions with   = 1 and/or    = 1 are set as Tabu active during  iterations.During each step , the best partition,   , is selected, and the search is repeated starting with that partition   .Each iteration is concluded by selecting the best overall partition among  1 , . . .,   .The suggested value of  is ⌊/2⌋, where  is the maximum number of vertices allocated to one of the two districts under consideration.To enhance diversity, probabilities of the districts being selected are adjusted during the iterations in order to enhance the search over spaces that have not been considered.

Hyperheuristic-LS (HypLS). Consider a codified partition 𝑥,
where   = 1 if vertex  is assigned to district .The neighborhood of  is a combination of 1() and 2() neighborhoods.The procedure randomly selects one of the two neighborhood structures and search over the space defined by the corresponding neighborhood.This procedure is repeated until a stopping condition is met, also maintaining a list with the three best partitions to perform a final attempt to improve them based on 1().

2-Iterations LS (2-IterLS
).Consider a partition , where   = 1 if vertex  is assigned to district .The neighborhood structure is the union of the 1() and 2() neighborhoods.Given a current partition, the procedure consists of evaluating all the possible ones according to each search, selecting the best partition over all the space.This procedure is repeated a number of iterations until a stopping condition is met.The three best partitions are also maintained so that a final attempt to improve them is done based on 1().
The best overall partition is reported as the output of the proposed algorithm.Notice that this corresponds to a feasible partition of the model proposed in Section 3; that is, it is a districting configuration that satisfies all the constraints and with best weighted sum of workload balancing and compactness metrics.

Numerical Experimentation
Numerical experimentation is performed to test the performance of all the procedures described above.The proposed algorithm is implemented in C programming language.We generate a vast set of instances of different types and sizes.For comparison purposes, the smallest size instances are solved with CPLEX 11.1.

Instance Generation.
We generate several types of instances in order to test the different components of the proposed algorithm under different conditions.The characteristics of each type of instance are summarized in Table 1.Instance sizes considered are small instances (50 vertices/5 districts and 200 vertices/10 districts), medium instances (450 vertices/15 districts), and large instances (1000 vertices/20 districts and 1500 vertices/30 districts).The Parcel instances have 1109 vertices and 28 districts.For all the instance types, the time required to perform a pickup and a delivery was set to a value of 10 and 5 minutes, respectively.For the test instances, the travel time to an adjacent vertex is assumed to be included in these values.
To generate the first three types of instances, the limits on the number of pickups and deliveries were defined based on two levels of capacity: tight () and less restricted (LR).The relative weighting factor varies over three values:  = 0.25, 0.5, and 0.75.The speed values vary over three values: 25, 30, and 35 kilometers/hour.Three replicates for each of the five instance sizes were generated.For each replicate, additional instances were generated by varying three factors: the value of the relative weighting factor, the two levels of capacity limits, and the three values of speed, resulting in a total of 3 × 2 × 3 × 3 × 3 × 5 = 810 instances.
To generate the urban type instances, we define the location of the vertices over the entire metropolitan region of Monterrey by generating a set of what we call base vertices, from which vertices are randomly selected to form each problem instance, for a large size instances (2211 vertices) and medium and small instances (1185 vertices).For each of the five instance sizes, three replicates were generated, and each of them was tested varying over the two capacity limits, the three values of the relative weighting factor, and the three values of speed, resulting in a total of 5 × 3 × (2 × 3 × 3) = 270 instances.
Parcel instances were generated considering GPS data provided by the parcel company that was collected during

Symmetric
Unrealistic instance but used to measure the performance of the heuristic for larger instances.Optimal solutions consist of "symmetric" districts: districts with the same workload content and the same shape.To generate the instance, a district is created and rotated to form the rest of the districts symmetrically.Workload content is also equal for all the districts.Euclidean distances are computed even if the points are not directly connected to guarantee symmetry.

Semisymmetric
Very similar to the symmetric instance, but Euclidean distances are computed only for those points connected by an edge.For the rest a shortest path distance is computed, and, hence, there is no guarantee that an optimal solution corresponds to a symmetric configuration.

Asymmetric
More general instances consist of vertices uniformly distributed over a plane using Euclidean distances for those vertices connected by an edge and shortest path distances for the rest.

Urban
Instances that consist of vertices distributed over a region that resembles the structure of the metropolitan area of Monterrey, Mexico.

Parcel
Instances generated with GPS data of a representative workday provided by the parcel company.This data set is a full size real instance.some representative days.The instance has a size of 1109 vertices and 28 districts.Each edge represents a real distance between a pair of adjacent vertices.The adjacency is determined respecting the real street structure of the city.Workload content (pickup or delivery) is randomly assigned to each vertex, and a single replicate is tested with each of the two capacity limits defined and over the three values of lambda and the average speed, resulting in total 2 × 3 × 3 = 18 instances.In total 810 + 270 + 18 = 1,098 instances are solved.

Stopping Rules and Control Parameters.
A limit time of 7200 seconds is set for CPLEX.For the heuristic and all the neighborhood structures explored we set a limit time of 3600 seconds.An iteration of the heuristic is allowed to proceed as long as it begins before the 3600-second time limit.The last iteration itself may take a variable amount of time to complete.Based on preliminary experimentation we define a maximum number of iterations for the FIS-Construction and for the local search procedures; the value of the Tabu permanence is set to four iterations, and the size of the RCL is set to a value of three.

Results and Discussion
. For the symmetric instances, the symmetric optimal solution is compared to the solution found by the heuristic and CPLEX.For the rest of the instances, a gap with respect to the best integer solution reported by CPLEX is computed as described by (11).Positive gaps are obtained when CPLEX finds better solutions: Tables 2-8 present results summarized either by type and instance size or by weighting factor.The tables present the maximum, average, and minimum value over the instances that were solved.Table 2 shows the gaps with respect to CPLEX results for the asymmetric, semisymmetric, and urban instances, considering those instance sizes in which CPLEX was able to find at least an integer solution (50 50 and 200 10).CPLEX found the optimal solution for the smallest size instances of 50 5.For the semisymmetric instances, all the tested solution methods found the optimal solution.For the asymmetric instances of 50 5, the 1-S method found the smallest of the maximum gaps (8.7%) but the 2IterLS found the best average gap (1.20%).For the urban instances of 50 5, the 2-IterLS method found the optimal solutions.Good results are observed in general with a maximum gap of less than 14%.
Table 3 shows the results for the asymmetric and urban instances by weighting factor and speed factors, for the instance size of 50 5, which CPLEX solved to optimality.We do not present results of the semisymmetric instances in the table because, for all the instances, all the tested solution procedures found the optimal solution.
As can be observed, the results do not vary too much according to the values of both factors and no significant pattern in the performance of the heuristics was found with respect to the values of the factors.Observe that, regarding the weighting factor, in general, maximum gaps were worse for an equal weighting factor for both criteria (0.5), but that, for the average gaps, some heuristics performed better with a low value of lambda (0.25) and others with a high value of lambda (0.75).Table 4 presents results for the same instance sizes and types, but with respect to the capacity level.As can be observed, the tight instances were more difficult for the procedures, and it was also more difficult to find feasible solutions for this instance type.
For the symmetric instances, the optimal solution (the symmetric solution) corresponds to symmetric districts with the same workload content and shape.Gaps are estimated similarly as in (11) but now considering the value of the optimal solution instead of the solution found by CPLEX.Table 5 reports the results found based on the instance size.We estimated gaps between CPLEX and the heuristics with respect to the optimal solution.CPLEX found the optimal solution only for the instances of 50 vertices.The positive gaps reported for CPLEX for the instances of 200 vertices indicate that it did not find an optimal solution within the time limit.Table 6 presents results for the symmetric instances based on the capacity level.Observe that the K-S/P procedure did not perform as well on 50 5 less restricted instances than on the tight instances.Also for the 200 10 instances, most of the heuristics, and CPLEX found better results for the tight instances.Overall, the 1500 30 instances were the most difficult to solve.
Table 7 shows the computational times for the rest of the instance types, by instance size.We present the maximum and average computational times in seconds for the instances that CPLEX could solve and for each of the heuristics.Low computational times are observed for all the instance types, with the heuristics reaching the maximum computational time only for the asymmetric and urban largest size instances of 1500 30, which are the most realistic and difficult to solve.Table 8 displays the results obtained with each method for the parcel instances.A gap is computed with respect to the current solution similarly to the gap presented in (11).The heuristic solution and the current districting configuration are compared evaluated with (3).Compared to the current districting design, all the procedures are able to find a better districting configuration, with the 2-IterLS heuristic finding the best results.Computational times were generally low.
Figures 1 and 2 each depict an example of a districting configuration.Figure 1 presents a graph in which districts overlap, a configuration that is not desirable for the company.In contrast, Figure 2 shows a configuration in which the districts can be distinguished from each other.The geography of the territory also determines the shape of the districts as can be seen in Figure 2, in which the districts conform to the urban geography of the region which is characterized by having several hills within and around the city.

Conclusions and Recommendations for Further Research
We present a mathematical formulation of the districting problem faced by a parcel company that serves a determined region.The model seeks to balance the workload content among the districts and to form districts of compact shape.The difficulty of the problem motivated us to develop a heuristic approach based on a hybrid heuristic that combines elements of GRASP and Tabu Search.To test the performance of the heuristic, we generate instances of several types and sizes, including a real instance using data from a parcel company.The results show good performance by the proposed approach, with good quality solutions and low computational times.The capacity factor is found to be the most significant factor because it conditions the difficulty of the instance to even find a feasible solution.
Further research includes formulating a stochastic version of the problem and analyzing different demand scenarios.The problem may also be solved as a biobjective optimization problem to find an efficient frontier instead of a single solution.Different mathematical programming solution approaches, such as Lagrangian Relaxation and decomposition, may also be explored.Although the heuristic is fast, it is not practical for the company to redesign the districts on a daily basis, because of managerial and human issues.Changing the districting configuration very frequently  could make it very difficult for the drivers to relearn the limits of their assigned districts, so a reasonable frequency should be defined to allow the drivers to adapt to their districts and avoid making the sorting process more subject to mistakes.It is important to consider that despite the speed of the heuristic, gathering data from the delivery and pickup scanners and the GPS devices requires considerable time because the system is not yet fully automated.In addition, not every day is a suitable representative day, which must meet several criteria and must be chosen with care.For these reasons, another future research topic could be to determine the optimal redesign frequency.

Figure 1 :
Figure 1: Illustration of districts that overlap.

Figure 2 :
Figure 2: Illustration of districts configuration for the parcel instance.
Function.It is similar to previous method but the region is divided into sectors of approximately equal workload content, measured by the number of vertices on it. fl instance of the problem (, (  ,   ),   ,   ,   , , ,   ,   ) Output:  = ( 11 ,  12 , . . .,  1 , . . .,  ||, ) fl A feasible partition of graph , or empty set if no feasible partition was constructed Set  = ⌀; set Seed = ⌀, and set  ∈ Ψ, ∀ ∈   fl instance of the problem (,   ,   ,   , , ,   ,   ) Output: Set Seed.Greedy Evaluation () RCL ← { best elements} Randomly select   ∈  among the elements of the RCL For the Neighborhood greedy function, discard from  the neighbors of   Else,  ← Semi-Random procedure, (neighbors of   are discarded from ) Else, randomly select   ∈  from among the discarded vertices End Return {} Algorithm 3: Pseudocode of the greedy randomized seed selection procedures.

Table 3 :
Asymmetric and Urban 50 5 results by weighting factor and speed factors.

Table 4 :
Asymmetric, semisymmetric, and urban results by capacity level.

Table 5 :
Symmetric instances results with respect to the optimal solution.

Table 6 :
Symmetric instances results by capacity level.