An iterative local search based hybrid algorithm for the service area problem

This article presents a hybrid algorithm for the service area problem. The design of service areas is one of the essential issues in providing efficient services in both the public and private sectors. For a geographical region with a number of small spatial units, the service area problem is to assign the service-demand units to the service-supply units such that each facility has a service area. The basic criteria for the service areas are the highest service accessibility, the contiguous service areas, and that the service demand does not exceed the service supply in each service area. A hybrid algorithm for the service area problem is proposed by extending iterative local search (ILS) algorithm with three schemes: population-based ILS, variable neighborhood descent (VND) search, and set partitioning. The performance of the algorithm was tested using 60 well-designed instances. Experimentation showed that the instances could be solved effectively and efficiently. The solutions found by the hybrid algorithm approximate optimal solutions or the lower bounds with an average gap of 0.15%.


Introduction
This article deals with the service area problem (SAP). The design of service areas is one of the essential issues in providing efficient services in both the public and private sectors (Daskin, 2011). The delineation of service areas for schools (Caro et al., 2004;Ferland & Guenette, 1990;Schoepfle & Church, 1991), hospitals and healthcare centers (Emiliano et al., 2017;Hu et al., 2018;Jacobs et al., 1996;Pezzella et al., 1981), disaster shelters (Hu et al., 2014;Li et al., 2008), green energy resources (Yanik et al., 2016) and many other facilities can be generalized as a contiguity-constrained capacitated facility SAP. For a geographic area with small spatial units, the SAP should assign the service-demand units to servicesupply units such that each facility has a service area and some criteria are satisfied. The basic criteria for the design of service areas are the highest service accessibility, the contiguous service areas, and that the service demand does not exceed the service supply in each service area (Daskin, 2011;Ferland & Guenette, 1990;Kalcsics et al., 2005;Kong et al., 2017;Xiao, 2008). Service accessibility can be evaluated by total travel distance from demand units to their supply units. The shortest travel distance is usually preferred when using the service. Contiguous service area is also a necessity for satisfying policy or management purposes. In addition, the total service demand in each service area should be no greater than its service capacity.
The SAP can be defined as a contiguity-constrained generalized assignment problem (GAP). It aims to minimize the total travel distance when using the service while satisfying constraints on facility capacity and service area contiguity. The GAP is known to be nondeterministic polynomial time hard (NP-hard). The constraints on spatial contiguity also pose great obstacles in modeling and solving the geographic problems (Duque et al., 2011). Consequently, various exact and heuristic methods have been proposed for solving the SAP.
Mixed integer linear programming (MILP) has been widely utilized to solve the districting problems. Owing to the computational complexity of the districting problems, many of the MILP models developed since the 1960s were simplified by ignoring some districting constraints such as district contiguity and unit integrity. The model solutions must be adjusted by repairing split units and fragmented districts. Such methods have been suggested for school districting (Caro et al., 2004;Franklin & Koenigsberg, 1973;Koenigsberg, 1968), political districting (George et al., 1997;Hess et al., 1965;Hojati, 1996;Li et al., 2007), and the territory design problem (Kalcsics et al., 2005). Modeling and repairing approaches could be adapted to solving the SAP. However, in case of a spatial mismatch between service demands and supplies in a geographic region, it may be difficult to find satisfactory solutions.
The second class of MILP models have considered district contiguity. Three types of formulations on district contiguity-tree-based, order-based, and flow-basedcan be embedded into an assignment model or a location-allocation model for solving districting problems (Duque et al., 2011;Shirabe, 2005). Nemoto and Hotta (2003) proposed a MILP model for political districting. The model considers district contiguity, but it does not guarantee the compactness of the districts. Salazar-Aguilar et al. (2011) proposed a bi-objective model for designing commercial territories. Constraints such as balance of sales volume and district contiguity are formulated in the model. However, 4 h of computation time were used to solve the problem instances with 150 units and 6 territories. Plane et al. (2019) proposed a minimum inter-person separation model for political districting. For instances with 6 × 6 grids can be solved optimally or near-optimally in 6.30-53,632.22 s. Kong et al. (2019) formulated a center-based allocation model with flow-based constraints on district contiguity to solve the districting problem. Given district centers from a K-medoids algorithm, the districting instances could be solved efficiently. Since the service-supply units are already provided in the SAP, the model could be modified to solve the problem.
The third exact approach has been suggested based on the set partitioning problem (SPP). Given a large number of feasible districts generated by the construction or heuristic methods, the SPP model attempts to select an optimal set of districts from the candidate districts as a districting solution (Garfinkel & Nemhauser, 1970;Mehrotra et al., 1998;Nygreen, 1988). Solutions from the problem instances indicate that the SPP model for the small-sized instances can be solved exactly. The model results also seem to be satisfactory for large problems. To obtain satisfying solutions, they suggested different techniques to generate promising candidate districts. Kong et al. (2017) proposed an iterated local search algorithm with set partitioning to solve the school districting problem. The SPP model could choose a much better solution from the historical districts identified by the local search. However, the algorithm was not fully tested on general SAP instances.
Local-search-based and population-based heuristics are the mainstream methods used for solving the SAP. Starting from an initial solution, the local-search based algorithms could iteratively improve the incumbent solution using neighborhood search (Kong et al., 2017;Li et al., 2008;Ricca & Simeone, 2008;Ríos-Mercado & Fernandez, 2009;Wang & Kong, 2021). The neighborhood space of the incumbent solution is explored by one or several search operators, and the incumbent solution is updated according to the acceptance criteria adopted in a specific algorithm. The one-unit shift is a widely used operator that moves a boundary unit to its neighboring district while maintaining the connectivity of its original district (Li et al., 2008;Liberatore & Camachocollados, 2016;Ricca & Simeone, 2008;Tavarespereira et al., 2007;Xiao, 2008). Butsch et al. (2014) introduced three operators for districting problems: shift, double shift, and swap. Kong et al. (2017) suggested much more complicated local search operators in iterative local search (ILS) algorithm for school districting. However, algorithms with expensive local search operators have not fully tested in the existing literature. Wang and Kong (2021) proposed a generic algorithm framework to implement the exact and hybrid algorithms for facility service districting.
Population-based heuristics have become popular for solving districting problems in recent literature. The evolutionary algorithms (Bacao et al., 2005;Bergey et al., 2003;Chou, 2011;Chou et al., 2007;Forman & Yue, 2003;Hu et al., 2014;Liu et al., 2016;Tavarespereira et al., 2007;Xiao, 2008) maintain and improve multiple candidate solutions using mechanisms such as crossover, mutation, and selection. The crossover operation for districting problems is usually implemented by overlaying two individuals and then repairing the fragmented overlaid areas (Chou et al., 2007;Datta et al., 2008;Liu et al., 2016;Xiao, 2008). The mutation operation usually moves every boundary unit to its neighboring districts with a small mutation possibility. To increase the speed of the algorithm convergence and the solution quality, Tavarespereira et al. (2007) proposed an evolutionary algorithm with local search for the districting problem. The population diversity in evolutionary algorithms and the search intensity in local search heuristics are balanced in the hybrid algorithm. In addition, natureinspired metaheuristics such as scatter search (Salazar-Aguilar et al., 2012) and artificial bee colony (Rincón-García et al., 2015) have also been applied to solve districting problems.
These are some limitations to the solution methods for the districting problems. First, the existing algorithms are capable of solving artificial or real problem instances; however, the performance of the algorithms has not been sufficiently tested on benchmark instances, so the optimality of solutions from the algorithms is still unknown. Second, some studies report that districting problems are difficult to directly solve using exact methods. However, along with progress in MILP, especially the embedded heuristics in mixed-integer programming (MIP) optimizers, state-of-the-art MIP solvers could solve increasingly difficult problems (Lodi, 2017). Designing new algorithms based on the framework of modern MIP solvers may be a promising approach to the districting problems.
This article presents a hybrid algorithm for the SAP. It was proposed by extending iterative local search (ILS) algorithm with three schemes: population-based ILS, variable neighborhood descent (VND) search, and set partitioning. The algorithm was tested using 60 welldesigned problem instances. The effectiveness and efficiency of the proposed hybrid algorithm for the SAP were verified using the problem instances. Compared with exact solutions, the solutions from the hybrid algorithm approximate optimal or the lower bounds with an average gap of 0.15%. To the best of authors' knowledge, it is the first time to introduce a population-based ILS with VND search and set partitioning for the SAP.

Problem formulation
For a geographic area, let V = {1, 2 … n} be a set of n small units. Each unit i has service demand p i . Let S = {s 1 , s 2 … s K } (S ⊂ V) be a set of K service-supply units, and each unit s k has service capacity q k . Let d ik be the distance between demand unit i and supply unit k. Let a ij indicate whether units i and j share a border and N i be a set of units that are adjacent to unit i (N i = {j| a ij = 1}).
A MILP model can be formulated by defining the SAP as a contiguity-constrained GAP. The contiguity of service areas can be ensured by a network flow model (Duque et al., 2011;Shirabe, 2005). Let x ik ∈ {0, 1} denote whether unit i is assigned to supply unit s k , H k (H k ≥ 0) be the service overload of supply unit s k , and f ijk (f ijk ≥ 0) be the flow from unit i to unit j in service area k. The model for the SAP is formulated as follows: Subject to : The objective function (1) is to minimize the total travel distance from service demand units to their supply units. It also penalizes service overloads by using a large enough coefficient α. Constraints (2) confirm that each unit i is assigned to only one service-supply unit. Constraints (3) are soft constraints on service capacities. Constraints (4), (5), and (6) are the flow-based formulations on contiguity. For two adjacent units i and j, if both of them are assigned to service-supply unit s k (x ik = 1 and x jk = 1), there might be a flow from unit i to j in basin k (f ijk ≥ 0); otherwise, these is no any flow between them (f ijk = 0). This is confirmed by constraints (4) and (5). Constraints (6) ensure that one unit of flow must be created at non-service-supply unit i if it is assigned to service-supply unit s k . Constraints (4), (5) and (6) enforce all the flows created in service area of k must run to the sink unit s k . This flow-based formulations were proposed for p-regions problem (Duque et al., 2011), and were rewritten for political districting (Kong et al., 2019) and facility service districting (Wang & Kong, 2021). Note that the redundant constraints on district sinks in Kong et al. (2019) and Wang and Kong (2021) is deleted from the SAP model. Constraints (7) and (8) impose restrictions on decision variables.
The soft constraints (3) on facility capacities are used in the model. The service overloads (H k ) in the constraints are penalized by the objective function. There are two advantages of using soft constraints. First, the soft constraints could help the MIP optimizer to find a feasible solution efficiently and guide the search toward preferred solutions. Second, using hard capacity constraints, there is no feasible solution for some instances, e.g. the total demand is greater than the total supply, or both the assignment constraints and capacity constraints cannot be satisfied. Such infeasible instances can be handled by the soft constraints.

Iterative local search
A hybrid algorithm for the SAP was proposed based on ILS algorithm. ILS is a simple, easy to implement, and quite effective metaheuristic for discrete optimization problems (Lourenço et al., 2010). ILS starts from an initial solution and iteratively improves the solution by local search and perturbation. Local search heuristic is used to intensively explore the solution space. However, the iterative search may easily get trapped in local optima that are far away from the global optimum. ILS escapes from local optima by applying perturbations to the current local minimum. The basic ILS procedure is illustrated in Algorithm 1. Note that f(s) in step 6 refers to the objective function of solution s.

Initial solution
The transportation problem (TP) model was revised to generate initial solutions (Wang & Kong, 2021). The model assigns the demand in each unit to one or multiple supply units with minimum traveling cost. The TP can be easily solved by a MIP solver owing to its bipartite structure. The small random coefficients ϵ ik (|ϵ ik | < 0.02) are used in the objective function to obtain a random solution.

Minimize
Subject to : X i∈V p i x ik ≤ q k; ∀k∈S ð11Þ The solution from the TP model needs to be repaired. First, some units in a solution may be split into two or more parts. For each split unit, the algorithm assigns it to a service area according to the largest portion of split. Second, some service areas in a solution may be noncontiguous. The non-contiguous areas can be repaired by deleting the fragmented units and then reassigning them to their neighboring areas. Note that the solution from region growth may violate the constraints on service capacities. The infeasible solutions are allowed in the algorithm. However, a solution with service overloads will be penalized in the following local search by the objective function (1) and tend to be a feasible solution.

Local search
The design of local search operators in optimization algorithm is critical for repeatedly improving the incumbent solution. The one-unit shift is a widely used operator that moves a boundary unit to its neighboring district while maintaining the connectivity of its original district (Kong et al., 2017;Li et al., 2008;Liberatore & Camachocollados, 2016;Ricca & Simeone, 2008;Tavarespereira et al., 2007;Xiao, 2008). Butsch et al. (2014) introduced three operators for districting problems: shift, double shift, and swap.
Two local search operators are used in the algorithm to improve the solutions. The operators attempted to move one or more units located on the boundary to their neighboring service areas. Note that only the feasible moves are allowed, because when moving a boundary unit from its original area to a destination area, the original area may be non-contiguous. The one-unit shift and two-unit shift are illustrated in Fig. 1.
The one-unit shift moves boundary unit i to one of its neighboring area k as outlined in Algorithm 2. In the algorithm, all the boundary units are selected (step 1) and then shuffled randomly (step 2). For each boundary unit i, the algorithm tries to move the boundary unit from its original area to one of its neighboring areas (step 3-7). The move is accepted in cases that the original area is connective, and the new solution is better than the incumbent solution (step 6 and 7).
The two-unit shift moves boundary unit i to one of its neighboring areas k, and at the same time moves boundary unit j in area k to one of its neighboring areas (Algorithm 3). In other words, one unit is moved into area k, and another unit in area k is moved out. Differing from the one-unit shift that involves two areas, a two-unit shift usually involves three areas. Swapping two units between two adjacent areas is a special case that involves two areas. The two local search operators have different search space and computational complexity: O(Kn) and O(Kn 2 ). However, the number of possible moves is generally much fewer, because only the boundary units can be moved to their neighboring areas, and there are only a few neighboring areas available for a boundary unit.

Perturbation
Perturbation is one of the key components of ILS algorithm. Using local search operators, the solution may easily reach to local optima. A ruin and recreate procedure is useful to perturb the SAP solution from local optima. There are multiple methods to ruin the solution such as deleting some boundary units, deleting all the units in some areas, and deleting some units in a connective region (Kong et al., 2017). Then, the deleted units can be reassigned to their nearest facilities. A repair procedure on the new solution is usually necessary since some areas in the perturbed solution may be non-contiguous. On the other hand, the objective of perturbed solution may become worse. However, the following local search could improve the solution quality significantly.
The strength of perturbation is critical for ILS algorithm. A very strong perturbation tends to generate a much worse solution, such that better solutions will be found with a very low probability. On the other hand, for a very small perturbation, the solution will often back to its original state of local optimum by the following local search. An appropriate strength of perturbation will maintain the diversification of ILS search, and thus benefits to find better solutions.
Four perturbation schemes are randomly used in our ILS algorithm: ruin boundary units and recreate, ruin service areas and recreate, ruin a connective region and recreate, and move boundary units to their neighboring areas. Let strength be the strength of perturbation, i.e. a percentage of the solution components, for example, strength% of boundary units, strength% of service areas, or strength% of all units. In all perturbation schemes, the perturbed units/areas are randomly selected.

Contiguity of service areas
To keep the contiguity of service areas in SAP solutions, several algorithms are frequently used in local search and perturbation. The first algorithm is to check the connectivity of a service area. If a spanning tree can be established using all the units in an area, the area is connective (Liu et al., 2016;Xiao, 2008). The second algorithm is to find the fragmented units in a solution. The spanning tree method can also identify the fragmented units. For each service area, a spanning tree is gradually constructed from the supply unit. If some units cannot be assigned to the tree, they are the fragmented units. The third algorithm is to repair a solution with noncontiguous service areas. The fragmented service areas are repaired by deleting the fragmented units and then gradually reassigning them to their neighboring areas.

A hybrid algorithm
We introduced a hybrid algorithm based on ILS to solve the SAP. The standard ILS algorithm was enhanced by three schemes: population-based search, variable neighborhood descent (VND) search, and set partitioning. The proposed algorithm is outlined in Algorithm 4.
The hybrid algorithm maintains a number of elite and diverse solutions. It is a population-based algorithm rather than the standard single solution-based ILS algorithm. The diversification of ILS search is guaranteed by the perturbation scheme. The population-based extension of ILS may enhance the diversification of local search. In step 1, a population of psize initial solutions are generated by the model (9)-(12). In step 6, a solution is randomly selected for perturbation and local search. In step 12, the population is updated according to the solution objectives and the similarities among the solutions. To maintain solution diversity, only elitist and dissimilar solutions are selected to be the new population. The updated population should be dissimilar. For any two solutions, their dissimilarity can be defined as the percentage of the demand units that are assigned to Kong Computational Urban Science (2021)  different facilities. Based on this definition, a dissimilar population means that the dissimilarity between any two solutions in the population is greater than a minimum dissimilarity. This percentage can be predefined as an algorithm parameter (mindiff). In some cases, the size of the new population may be smaller than the parameter psize; however, the following process of perturbation could generate new dissimilar solutions. The second scheme is to use VND local search instead of simple local search in ILS. A simple local search is to improve the solution by performing one-unit shift and two-unit shift respectively. VND is repeatedly exploring different neighborhood structures by deterministic change of neighborhoods until the solution cannot be improved. Using VND search, ILS algorithm can easily find solutions in local optimum, and then moves the solution to new search point by perturbation. The VND search is outlined as follows.
The third scheme is to improve ILS solution by set partitioning. The set partitioning has been successfully coupled with ILS to solve districting problems (Kong et al., 2017;Wang & Kong, 2021). Let Ω denote the set of historical service areas identified in ILS algorithm. Each area i has an objective o i and a set of units U i in the area. Let subset Ω j = {i| i ∈ Ω, j ∈ U i }, the SPP model (13)-(15) could be used to select a subset service areas of Ω, which gives a minimal objective and also covers all the units in set V. The decision variable x i indicates whether the candidate area i is selected.

Minimize
Subject to : x i ∈ 0; 1 f g; ∀i∈Ω ð15Þ In each ILS loop, all the service areas in solution s* are recorded in the list pool (step 11). Each record includes the units in the service area and its objective value. At the end of ILS, thousands or more areas will be recorded in the pool. The SPP model could be used to select a better solution from the candidate service areas (step 14).
In algorithm implementation of the hybrid algorithm, the ruin methods described in section 3.4 are randomly used in step 7. The one-unit move and two-unit move local search operators are used in VND search (step 8).

Experiment settings
Three study areas, ZY, GY, and GY2 (Kong et al., 2019) were adapted for delineating the service areas. In urban area ZY, there are 324 spatial units and 3783 primary school students. In rural area GY, there are 297 spatial units and 36,824 junior middle school students. In area GY2, the same area of GY, there are 1276 residential communities /natural villages and 819,812 residents. In this article, the number of students or the number of residents in each spatial unit is assumed to be the quantity of demand. Figure 2 shows the spatial units and the service demand in each unit.
Using the three study areas, 60 instances for the SAP were prepared to test the hybrid algorithm. For each study area, 20 instances were designed as follows. (1) The schools (or township centers) are supposed to be the candidate service-supply locations, and the number of students in each school (or the number of residents in each township) is assumed to be its service capacity. In addition, more units are manually selected as candidate service-supply locations; their capacities are randomly set to be a number between the minimum and maximum capacities of schools (or township centers). Consequently, 36, 44 and 24 units in the three areas were selected to be the candidates of service-supply units with total supplies of 9195, 42,326 and 1,324,763, respectively.
(2) Select a number of facilities from candidates randomly. 13-17, 37-41, and 18-22 service-supply units were selected for the three areas, respectively. (3) Adjust the facility capacities in each instance proportionally to ensure that the supply-demand ratio is around 1.15 and 1.03, respectively. (4) Select a number of facilities from candidates by solving a capacitated p-median problem, and then adjust the facility capacities with the same supply-demand ratio in step (3). Table 1 shows the attributes of the instances. The 20 instances for each study area can be classified into 4 sets: random facility locations with higher supply-demand ratio (A), p-median facility locations with higher supplydemand ratio (B), random facility locations with lower supply-demand ratio (C), and p-median facility locations with lower supply-demand ratio (D). As a result, the instances are diverse in terms of the number of demand units, the number of facilities, the facility locations, the supply-demand ratio and the average number of units served by one facility.
The proposed algorithm was implemented by using the Python programming language. The Python script were run in PyPy 6.0, a fast and compliant implementation of the Python language (see http://pypy.org). Each instance was repeatedly solved for ten times. Since random mechanism was used in initial solution generation, local search and perturbation, different solutions will be obtained by repeatedly executing the algorithm. To verify the optimality of the solutions, the instances were also solved by the MILP model formulated in Section 2. The algorithm ran on a desktop computer with Intel Core I7-6700 CPU 3.40-GHz, 8-GB RAM and the Windows 10 operating system.
In the experiment, algorithm parameters are shown in Table 2. The coefficient α in objective function (1) was set to 10,000. It is large enough to penalize the service overload in some service areas. To analyze the sensitivity of the parameter settings, the instances were also solved by changing one of the algorithm parameters. The SAP, TP and SPP model were solved by the IBM ILOG CPLEX Optimizer 12.6.3. Two parameters were set for the CPLEX optimizer for solving the SAP model: the relative MIP gap tolerance MIPGap = 10 − 10 and optimizer time limit Timelimit = 7200 (seconds).

Solution results
The solutions for 60 instances are shown in Table 3. The instance name in the table consists of the area name, the number of facilities and the type of facility configurations. For each instance, the lower bound, upper bound, and computation time obtained from the CPLEX optimizer are shown in columns LB, UB, and Time. The upper bounds labeled with asterisk are optimal objectives. For each instance, there are 10 solutions obtained from the hybrid algorithm. The best, average and worst objective gaps are shown in columns Gap min , Gap avg and Gap max respectively. The gaps are calculated by the formula (obj-LB)/LB*100%. Note that the gap value is 0 in case of the solution is optimal. Column Dev shows the standard devotion of the objective values. Column Time indicates the average computation time in seconds. Note that all the solutions are feasible in terms of capacity constraint and contiguity constraint. Table 3 shows that all the SAP instances can be solved optimally or near-optimally by CPLEX optimizer. Among the model solutions, 46 are optimal and 14 are near-optimal with MIPGap less than 0.32%. Different from the early reports that only very small districting instances could be exactly solved, large instances for the SAP could be solved in our experiment. It is also found that the instance complexity depends on both the instance size and the supply-demand ratio. The instances with high supply-demand ratio (instance sets A and B) would be solved more efficiently than those with low supply-demand ratio (instance sets C and D).
The results in Table 3 also show that the proposed hybrid algorithm is effective and efficient. First, the solutions from the hybrid algorithm approximate to the optimal solutions or the lower bounds with an average gap of 0.15% and ranging from 0.00% to 0.82%. Among the best solutions in column Gap min , 28 (46.7%) are optimal. Table 4 summarizes the number of optimal solution found by CPLEX and the hybrid algorithm. Second, the objective deviations are small, ranging from 0.00% to 0.38%. Third, all the instances were solved by the hybrid Fig. 2 The spatial units in areas ZY, GY and GY2 (Kong et al., 2019). The grey circles indicate the size of the service demand in each unit Set partitioning in the hybrid algorithm is capable of improving the ILS solutions. The average improvement on all instances is 0.12%. Table 5 illustrates the improvements for each set of instances. For area GY, the instances of type C and D are significantly improved by set partitioning: 0.60% and 0.43%, respectively. Occasionally, the ILS heuristic may be difficult to find a feasible solution; however, the set partitioning procedure could select a feasible solutions with high quality. The computation times for solving SP models are shown in Table 6. The SP models for SAP instances could be solved very efficiently. Generally, the SAP solutions could be improved by set partitioning in a short computation time.

Analysis of algorithm parameters
To evaluate the sensitivity of the parameter settings, all the instances were solved by the proposed algorithm using 9 additional sets of parameters. The solution gaps are summarized in Table 7. Column Default shows the solution gaps from the algorithm with parameters listed in Table 2. Other columns shows the solution gaps from the algorithm by setting different parameters: P1: default parameters, but set population size to 1 (psize = 1); P5: default parameters, but set population size to 5 (psize = 5); P20: default parameters, but set population size to 20 (psize = 20); R10: default parameters, but set perturbation strength to 10% (strength = 10%); R15: default parameters, but set perturbation strength to 15% (strength = 15%); OP1: default parameters, but use one-unit shift search operator only; OP2: default parameters, but use two-unit shift search operator only; LS: default parameters, but use simple local search in ILS loops; RG: default parameters, but use region growth method to generate initial solutions.
There are several findings from the solution results in Table 7. First, solution gaps in columns Default, OP1 and OP2 show that the two-unit shift search operator is more effective than the one-unit shift; and the combination of one-unit shift and two-unit shift operators is slightly better than the two-unit shift operator. The oneunit shift is a fast local search operator, while the twounit shift is slower but more effective. For designing the local search procedure in ILS, an issue is that whether it is best to have a fast operator or an effective one. The experimentation suggest that the best choice is to use both of the one-unit shift and two-unit shift in ILS search.
Second, solution gaps in columns Default, P1, P5 and P20 indicate that the population-based search (Psize = 5, 10 or 20) is better than the single-solution-based search (Psize = 1). Maintaining a population of elite and diverse solutions, local search in ILS could explore much larger neighborhood space, and thus have more possibility to find better solutions. It is observed that the populationbased ILS performs better than the single-solution-based ILS, and the population size is not sensitive to solution quality.
Third, other parameters have some effects on the performance of the proposed algorithm, but are not very sensitive to the solution quality. The perturbation strength is a key parameter to balance the intensification and diversification of ILS search. In our algorithm, the perturbation strength between 5%~15% is appropriate for solving most instances. For ILS local search, both the simple local search and VND search could guide the algorithm to find good solutions. Since the VND search is capable to reach local optimum in solution space, the ILS with VND search might be better than the ILS with simple local search. The final solutions for the SAP   slightly depend on the initial solution method. For most instances, the procedure of solving a TP model and repairing the model solution is the best mothed to generate initial solutions.

Conclusion
A new hybrid algorithm was proposed to solve the contiguity-constrained capacitated facility SAP. To the best of the author's knowledge, it is the first time to introduce a population-based ILS with VND search and set partitioning for delineating service areas. The performance of the proposed method was intensively tested on 60 well-designed instances. Experimentation shows that the algorithm is capable of finding high-quality near-optimal solutions in a reasonable computation time. The proposed hybrid algorithm is different from existing algorithms in some aspects. The algorithm was designed based on the standard ILS algorithm, and further enhanced by three schemes: VND search, population-based ILS, and set partitioning. VND search in ILS could extensively explore the neighborhood space of the incumbent solution, and thus increases the convergence speed of ILS. The population-based ILS maintains a set of elite and diverse solutions. Starting from the multiple initial solutions, the perturbation and local search in ILS have a higher possibility of finding better solutions. Additionally, the service areas explored by the local search are reused for selecting a better solution by set partitioning. In traditional set-partitioning-based heuristics, even if different techniques are suggested to generate a large number of service areas, it is still difficult to identify the promising candidates for large problem instances. In local-search-based or population-based heuristics, a large number of solutions are discovered; however, most of them are abandoned. Historical solutions in metaheuristics should be the high-quality candidates for the set-partitioning-based approaches. In our algorithm, all the service areas identified in each iteration were recorded and finally reselected by a SPP model.   Twelve sets of benchmark instances with different sizes and complexity for SAP were introduces based on three geographic areas. The instances are diverse in terms of the number of demand units, the number of facilities, the facility locations, the supply-demand ratio and the average number of units served by one facility. The optimal or near optimal solutions for the instances were obtained from CPLEX optimizer. The instances and all the solution results shown in this article can be downloaded from https://github.com/yfkong. The authors believe that the benchmark instances are valuable to evaluate existing and newly proposed algorithms for the SAP.
Finally, some general issues regarding the hybrid algorithm should be investigated in further studies. The parameter settings and the adaptive adjustment of the parameters for the algorithm remain as open issues. Both the theoretical and practical analysis on this topic may inform the algorithm to solve a specific instance more efficiently. Considering that districting problems have some common criteria, it could be possible to extend the algorithm to solve other problems such as the political districting problem and location-districting problems.

Code availability
The code is available to download at https://github.com/yfkong/AFSAP.

Author's contributions
The author wrote the whole manuscript and read and approved the submitted version.

Funding
This work is supported by the National Natural Science Foundation of China (41871307).

Availability of data and materials
The data is available to download at https://github.com/yfkong/AFSAP.