A Novel Spatial-Temporal Voronoi Diagram-Based Heuristic Approach for Large-Scale Vehicle Routing Optimization with Time Constraints

: Vehicle routing optimization (VRO) designs the best routes to reduce travel cost, energy consumption


Introduction
Flexible transportation service not only reduces travel cost, but also alleviates related energy and environmental concerns, such as traffic congestion, energy consumption, and carbon emission [1,2].With the advancement of geographic information science (GIS), space-related problem decision-making has been embedded in spatial decision support systems (SDSSs) to provide flexible transportation service in both public and private sectors [3], such as the design of school bus routes [4], collection of solid waste [5,6], distribution of goods in chain business [7], design of police patrol routes [8], and planning express deliveries [9].Spatial function tools are used for intelligent transportation operation in many services including management of geo-referenced transport network data, geocoding massive human demands, positioning moving vehicles, and navigating routes for drivers [5][6][7][8][9].
Vehicle routing optimization (VRO) aims to design least cost routes to satisfy geographically-distributed human demand through spatial intelligence of SDSSs in the field of transportation [10], including the shortest path problem (SPP) [11][12][13], traveling salesman problem (TSP) [14], and the vehicle routing problem (VRP) [5][6][7].One type of VRO is modeled as the vehicle routing problem with time window (VRPTW), which imposes time constraints on customers, vehicles, and depots to accept or distribute high quality service [4,6,14].It specifies that a customer only accepts certain transportation-related service in a given time interval [15][16][17].The VRPTW is computationally intractable in many real-life applications [18].The heuristic algorithm is the only justified approach as it can find good solutions in a reasonable time [16,17].There are two types of heuristic algorithms for the VRPTW, which include local search and evolutionary optimization.Local search iteratively modifies local route structures to improve solution quality and follows a single trajectory to explore the solution space [3,6,7].Evolutionary optimization simultaneously searches multiple parts of the whole solution space to converge on the best solution [5,19,20].Recently, state-of-the-art heuristic algorithms have been reported to yield a high quality solution of the VRPTW with hundreds of customers [5,9,18].More efficient heuristics are needed for even larger VRPTW applications.
Spatial principle has exhibited good performance in space-related problem decision-making [21][22][23] in problems such as bus stop location [3], hazmat route planning [20], and path covering [24].To efficiently solve VRO, from a spatial perspective, a few speed-up strategies have been developed to handle spatial issues and extend solvable VRO instances up to thousands of customers, which include the granular neighborhood [25], the kth nearest neighbors [26], and the k-ring-shaped Voronoi neighbors [27].The common denominator of these effective strategies is to search only a limited number of spatial neighbors to save computing efforts [28].However, such useful strategies become ineffective for the VRPTW since time constraint imposes an exceptional influence on the proximity of consecutively served customers in a route.A long wait will occur if two spatially-close customers with quite different time windows (for example, one is 8:00-9:00 a.m. and the other is 3:00-4:00 p.m.) are visited consecutively in a route, which challenges the motivation of spatial proximity-based heuristics approaches.Hence, spatial proximity measures must be extended to address the additional temporal constraints [18].Moreover, how to accelerate the VRO with the help of spatial-temporal thinking should be investigated.
Voronoi diagrams describe both spatial proximity and topological proximity between spatial objects [29][30][31][32].By extending it to the space-time context, a spatial-temporal Voronoi diagram (STVD) is well suited for measuring the proximity of customers with time constraints.Therefore, this article proposes a spatial-temporal Voronoi diagram-based heuristic approach to solve very large-scale VRPTWs.A derived spatial-temporal Voronoi neighbors from the Voronoi diagram is used to select promising candidates for a local search.Furthermore, motivated by the truth that nearer neighbors have a higher possibility to be consecutively served in a route [28], a Voronoi distance decay strategy is proposed to assign reasonable search efforts on neighbors with different Voronoi distances.A spatial-temporal feature guided search is used to reconstruct unpromising micro route structures.Intensive experiments on benchmark datasets (25-1000 customers) and real-world large scale VRPTWs and its multi-depot variants (2000-10,000 customers) have been implemented to assess the performance of the proposed approach.This efficient and effective approach enables current local search heuristics to solve super large-scale VRPTWs and will contribute to improving spatial intelligence for transportation-oriented SDSSs.
The remainder of this paper is organized as follows.Section 2 reviews related literature.Section 3 introduces the VRPTW, the STVD, the derived spatial-temporal Voronoi distance and neighbors.The proposed spatial-temporal Voronoi diagram-based heuristic for the VRPTW and its multi-depot variants are described in Section 4, and details of the experiments, results, and comparisons are reported in Section 5. Section 6 discusses the impact of the STVD and spatial-temporal proximity behind the best found solutions.The last section draws the conclusions.

Literature Review
This section reviews heuristic algorithms for the VRPTW and the integration of vehicle routing optimization and GIS.

Heuristic Algorithms for the VRPTW
Heuristic algorithms for VRPTWs are divided into two categories: local search and evolutionary optimization.Starting from an initial solution, local search shifts from a current solution to a better neighborhood solution with a simple change of local routes [33].The change is made iteratively by moving nodes or exchanging arcs within a route or between routes such as 1-0 exchange, 1-1 exchange, 2-Opt, and λ-opt.Usually, the change process becomes trapped in a local optimum.To overcome unpromising results, many intelligent strategies that accept some worse solutions have been developed to further improve the solution quality, such as simulated annealing [34], iterated local search [35], large neighborhood search [36], variable neighborhood search [37], and tabu search [38].Numerous tests on small-and medium-size VRPTWs with 25-100 customers have demonstrated good performance of local search heuristics [34,36,37].
Evolutionary optimization concurrently improves many solutions and results in high quality solutions.Four major steps are involved: representation, selection, combination, and mutation.Some of the most successful evolutionary optimization algorithms for the VRPTW are provided by Repoussis et al. [39] and Vidal et al. [40].Mester and Brä ysy [41] used the standard evolutionary optimization framework to guide exploration in the VRPTW solution space with solution initialization and evolution.Gehring and Homberger [42] parallelized the genetic algorithm and first solved VRPTW instances up to 1000 customers.Usually, local search is used as a component in the evolutionary optimization approach.Fast and simplified node or arc exchange-based local search have been used for offspring education in the genetic algorithm [41].Therefore, more efficient local search heuristics also help to improve evolutionary optimization.
For larger VRPTW instances, recent advances have proposed several accelerating techniques to save computing efforts while still achieving high quality solutions.Following the local search principle, taking spatial distance into account, Cordeau [38] limited the local search candidates with kth nearest neighbors to shorten computing time.Toth and Vigo [25] selected the local search evaluated nodes with a fixed distance threshold.Considering time constraints, Ibaraki [43] proposed a time-oriented neighbor list-based strategy to reduce the local search evaluation.Some unpromising search candidates had been excluded with the unmatched time windows.But the spatial proximity has been ignored.Nagata [44] developed a time warp operation for customers with a late service to avoid complex time window analysis.An additional penalty function was added to the original objectives to justify the local search heuristics.However, the spatial dimension and the temporal dimension are still separated in these local search heuristics.
Another effective accelerating approach is to decompose large-scale VRPTW instances into many smaller subproblems.Using spatial clustering, Dondo and Cerdá [45] divided the large scale VRPTW into a set of small-sized problems.Routes were generated for each small problem to provide a high quality initial solution.Qi et al. [18] developed an effective spatial-temporal distance measure to deal with both space and time issues.A large-scale VRPTW was partitioned with a k-medoid clustering method to generate a good initial solution.A standard genetic algorithm was used to improve the initial solution.Good performance of these effective approaches has been verified by VRPTW instances with no more than 1000 customers.
With regard to the Voronoi diagram, the only work has been provided by Milthers [46] who used two-dimensional Voronoi diagrams to split the VRPTW into many subproblems and then solve them with large neighborhood search heuristics.It showed the effectiveness of Voronoi diagrams in guiding the search process.However, this study dealt only with decomposition of the VRPTW in the solution construction stage.Use of the spatial-temporal Voronoi diagram in speeding up the local search for the VRPTW should still be further considered.

Integration of Vehicle Routing Optimization and GIS
Effective approaches to find high quality solutions for VRPTWs have been embedded in a GIS-based SDSS to deal with many real world applications.Generally, the integration between GIS and vehicle routing optimization is tightly coupled.Spatial data management, processing, and visualization tools are used to collect customer orders, georeference related data, activate the solving process, and display routes for VROs, as in GIS software such as ArcGIS and TransCAD.In line with local search, Weigel and Cao [7] first reported their efforts in the implementation of a tabu heuristics approach in a GIS environment to deal with VRO for a major American retailer.Following the evolutionary optimization line, Mendoza et al. [5] integrated a customized routing module, which improved solution quality with commercial solutions such as SAP/R3 and ArcGIS to cope with VROs in public utilities.Experiments on a real-world case in Bogotá , Colombia with 323-601 customers verified the effectiveness of evolutionary optimization and GIS software.
Recently, progress of GIS related to cloud computing transfers the integration of GIS and VRO to a loose coupling way.Using Google maps, Santos et al. [47] developed user-friendly web-based SDSSs that embed VRO to deal with urban trash collection in Coimbra, Portugal.TU et al. [48] presented a cloud GIS-based spatial decision support framework with variable neighborhood search heuristics for dynamic vehicle routing using historical traffic information.These practices have proven the dominance of VROs in real-life transportation applications.However, they also indicate that current spatial intelligence should be improved to enhance the ability of SDSSs in transportation departments to cope with increasing numbers of customers.
In summary, state-of-the-art heuristics, including local search and evolutionary optimization, solve small-and medium-size VRPTWs in a reasonable time.They have been verified to be effective in a GIS environment for real-life transportation cases with no more than 1000 customers.For larger-size real-life applications, more efficient heuristic algorithms are needed to reduce computing time.Using the spatial-temporal Voronoi diagram, this paper proposes an efficient local search-based approach to address VRPTWs up to 10,000 customers in a reasonable time.Differing from the neighborhood in the spatial context [25][26][27]48,49], the spatial-temporal Voronoi neighborhood considering both spatial distance and time windows is proposed to limit the search space.A Voronoi distance decay strategy is developed to assign more searching efforts to nearer neighbors.The spatial-temporal feature guided search is used to escape from local minima.Details of the STVD, Voronoi distance, Voronoi neighbors, and the proposed approach for the VRPTW are described in the next sections.

The VRPTW and the Spatial-Temporal Voronoi Diagram
This section introduces the VRPTW, the STVD, and the spatial-temporal Voronoi distance.

The VRPTW
The VRPTW has five components: customers, depot, vehicle, routes, and the solution.Customers: a set of customers   = { 1 ,  2 , … ,   } expects to be served.A customer vi (I > 0) located at (xi, yi) has a unique demand, qi, a service duration, si, and a service time window <ei, li>.The service is required to arrive after start time, ei, and before end time, li.Such a time window is termed "soft" when it could be violated with a penalty cost, and termed "hard" when no violation is permitted.
Depot: The depot  0 located at (x0, y0) provides service in a time window <e0, l0>.Vehicles: A fleet of K identical vehicles with the capacity Q is located at the depot.Each vehicle departs from the depot, serves multiple customers, and returns to the depot.Travel distance and time between nodes (includes customers and the depot) i, j are denoted as dij and tij, respectively, where  ≠ .

Routes:
A route  =<  0 ,  1 ,  2 , . . .,   ,  +1 > is formed by a sequence of visited customers, which represents the service process of a vehicle, where R denotes the number of visited customers. +1 =  0 , so that the vehicle must return to the depot.The earliest departure time at the depot  0 is   , the earliest service time at a customer   is    , and the earliest return time at the depot is   +1 , which are defined recursively as follows: A route is specified following two restrictions: (1) should be no more than the vehicle capacity Q.The total travel distance of the route is defined as ∑     +1  =0 . A route is defined as feasible if the customer's service time window, the route duration time, and the vehicle capacity are satisfied.
The VRPTW solution: A solution for a VRPTW instance is a set of feasible routes  1 ,  2 ,  3 , . . .,   such that a single vehicle visits each customer once.A hierarchy of objectives of the VRPTW is followed in this paper.The number of used vehicles is first minimized and then the total travel distance is reduced.

The Spatial-Temporal Voronoi Diagram for VRPTWs
Components of the VRPTW, such as customers, depot, vehicles, and routes, have both spatial and temporal attributes.Spatial attributes refer to customers' locations, the depot location, and the vehicle trajectories.Temporal attributes refer to the time window, service time, travel time, and route duration time.To represent them simultaneously, a spatial-temporal coordinate system is built by taking the temporal dimension to be perpendicular to the planar XY plane.Therefore, time windows of the customers and the depot are modeled as vertical lines rooted in their locations.
Figure 1 illustrates a small-size VRPTW with five customers (v1 to v5), one depot (v0), and two routes (R1, R2).Customers (v1-v5) and the depot (v0) are modeled as vertical lines from its start time point to end time point, as line segments L1-L5 and L0.Projections on the XY plane denote their spatial locations, whereas projections on the time dimension denote their time windows and actual service time.The transportation between customers denotes the service process in the routes (R1 or R2).
Distance in this coordinate system is defined as Equation ( 2), where (xi, yi) and (xj, yj) denote the spatial location of i and j, and ti and tj denote the time at i and j, ̅ is the mean travel speed: Voronoi diagrams partition space into adjacent Voronoi cells with given points [30].Following this principle, a spatial-temporal Voronoi diagram divides the spatial-temporal space into many cells.To generate a Voronoi diagram for the VRPTW, the spatial-temporal central point of a customer's spatial temporal line, (xi, yi, (ei + li)/2.0),which is the red point in Figure 1, is used as the seed point.The STVD for VRPTWs is built by the fast quickhull algorithm [50].
In the STVD, a pair of points whose cells share a boundary (Voronoi face or Voronoi edge) are Voronoi neighbors [30].The spatial-temporal Voronoi distance from a given point i to the point j, vd(i, j), is defined as the minimum number of crossings of the spatial-temporal Voronoi boundary in any route from i to j [27,31].Figure 2 provides an example that illustrates the Voronoi distance.The Voronoi distance from P1 to P3 is 2 for the two crossings of Voronoi boundaries (C1-C2).The Voronoi distance from P1 to P13 is also 2 for the crossings C3-C4.The Voronoi distance from P1 to P12 is 3 for the crossing C5-C7.Thus, Voronoi neighbors with different spatial-temporal proximity can be measured.As both spatial distance and time windows are considered, it is well suited for accelerating local search for the solving of VRPTW.Using these useful definitions, this paper proposes an efficient local search heuristic to address large-scale VRPTWs.

The Spatial-Temporal Voronoi Diagram-Based Heuristic
This section introduces the spatial-temporal Voronoi diagram-based heuristic (STVDH) for large-scale VRPTWs.The principle of the STVDH is to speed up the local search procedure by using useful spatial thinking, including the spatial-temporal Voronoi neighbors, distance decay search strategy, and local route features.The overall framework of the proposed heuristic is summarized in Figure 3.After the construction of an initial solution (step 1-1), a local search heuristic based on spatial-temporal Voronoi neighborhood is followed to improve solution quality (step 1-1 to step 1-4).A spatial-temporal Voronoi distance decay strategy is developed to assign reasonable searching efforts on neighbors with different distances (step 1-2).Spatial-temporal route features are identified and guides the searching process escape from local minima (step 1-3).The stopping condition is a limitation on the number of iterations (step 1-4).Finally, the best found solution is reported (step 1-5).

The Construction Algorithm
The construction algorithm iteratively inserts an unrouted node into a proper position until all customers are routed.During the insertion, instead of evaluating all unrouted nodes, neighbor customers within a given Voronoi distance are selected as candidates of the next inserted nodes.The details of the construction algorithm are described in Figure 4. First, customers are sorted by the start time window and stored in a queue L (step 2-0).Then, L's front unrouted node is popped out as the operated node  (step 2-1).An empty route r is initialized with u (step 2-2).Another unrouted customer v is randomly selected from ' Voronoi neighbors.Node v will be inserted before or after a routed node i in the route r (step 2-3).If the insertion is performed, u is updated with v.The best insertion place is defined as the one that minimizes the added total travel length and the total value of the end time window violation as Equation ( 3), where j denotes the next visited node of i.  is a random value in [0.5, 2].It should be noted that the arrival time at v, j should not be before the start time, as Equations ( 4) and ( 5 During the insertion, several constraints including the start time window, route duration time, vehicle capacity, and the maximum route length should be satisfied.If no place is found for insertion, we pop L's front node as u, insert it into a new route, and label it as routed (step 2-2 and step 2-3).The insertion will end if all customers have been labeled as routed (step 2-5).Finally, obtaining all routes is the initial solution.

Local Search
Local search iteratively explores the solution space to improve solution quality.Usually, the exploration removes or adds edges from the current solution by using local search operators.Three typical local search operators, including the 1-0 exchange move, 1-1 exchange move, and 2-Opt move [51], are used to explore the neighborhood solution in this article.
 1-0 exchange move injects a node from the original place and inserts it after another.As Figure 5a shows, node b is removed from its current position and inserted before node i.  1-1 exchange move swaps positions of a pair of nodes.As Figure 5b shows, the places of nodes b and i are swapped.
 2-Opt move swaps the ends of two routes after the positions of a pair of nodes.As Figure 5c illustrates, the partial routes after b and i are swapped.The computing complexity of these three operators depends on two aspects: operated nodes and local search evaluation.Two operated nodes, named b and i, are involved as shown in Figure 5. Therefore, the number of operated candidates is n(n − 1), where n is the number of customers.The time complexity is ( 2 ).For each candidate pair, local search evaluation checks the time window for each customer involved, total duration time, and vehicle capacity for each route.Due to the cascade effect of the arrival time, it will, on average, take (/2) time, where R is the number of routes.To utilize limited computing effort more efficiently, we propose a spatial-temporal Voronoi distance-decay strategy to assign searching efforts on selected candidates.A time warp operation is also used to shorten computing time for each candidate.

The Spatial-Temporal Distance Decay Strategy
Differing from the equal searching intensity of the k nearest neighbors [26] and the k-ring Voronoi neighbors [27,48], the spatial-temporal Voronoi distance decay strategy searches more on nearer neighbors than further neighbors to balance the solution quality and computing time.For spatial-temporal Voronoi neighbors with Voronoi distance k, the searching probability Pk is defined as Equation (6).As the distance k increases, searching probability Pk will decrease.Therefore, most computing effort will focus on near neighbors, but some necessary efforts are left for further neighbors.
where k = 1, …, kmax, and kmax is the maximum searchable Voronoi distance,  ∈ [1, ∞) is a value that controls the spatial-temporal distance decay speed.As the value of  increases, the decay speed becomes faster.
To implement the Voronoi distance decay strategy, local search randomly generates a value of δ between (0, 1), finds the probability range in which δ is located, and then selects Voronoi neighbors with the corresponding distance k as candidates to evaluate.As it avoids evaluating all neighbors, this useful strategy puts more searching efforts on closer nodes and therefore significantly accelerates the local search.As the number of Voronoi neighbors is fixed, the computing complexity reduces from ( 2 ) to ().

Time Warp Operation
The local search evaluates distance and time constraints for customers and routes, which cost the most computing effort [43].Nagata [45] proposed a new penalty function in which change can be computed in O(1) time for the VRPTW.When a later service happens for a customer, a time warp operation is performed such that the arrival time is adjusted to the end of the time window, and a penalty is added to the original objective.Figure 6 shows an example of the time warp.The additional penalized cost is defined as the time warp used in Equations ( 8) and ( 9), where twi is the cost value at customer i, and (ai − ei) + is a non-negative value of the service late time at i. Therefore, the objective is changed to minimize the total route number, total route length, and total penalized cost TW(s).For O(1) time consumed, the time warp operation is very fast.

Acceptance Criterion
Local search will naturally drop into local minima.A threshold acceptance criterion is used to escape.If the found solution is not 0.5% worse than the current solution, it will be accepted as the new solution.
Otherwise, it will be rejected.If the best solution does not improve in Zmax iterations, the spatial-temporal features-guided search will be started.

Spatial-Temporal Features-Guided Search
Due to the used spatial-temporal Voronoi neighbors, local search explores a relatively small proportion of the neighborhood solutions.Some potential neighborhood solutions will be rejected due to spatial and temporal constraints, such as route duration, time window, and vehicle capacity.Different from traditional random perturbations [37], spatial-temporal features guided search destroys unpromising structures defined by spatial-temporal rules and rebuilds a part of the current solution to escape local minima.It also increases the searchable probability for farther neighborhood solutions.Three components are the spatial-temporal features, the removal algorithm, and the reinsertion algorithm.

Spatial-Temporal Features
Three types of spatial-temporal features are defined in this paper including time window violation nodes  , longest distance nodes   , and smallest route's nodes   .Their definitions are described below.

 Time-window violation nodes
The time warp operation will introduce time window lateness in which the arrival time at a customer may be later than the end time of its time window.However, this route feature should not belong to the best solution.We remove such unpromising nodes from the current solution and reinsert them later.Formally, time window violation nodes   are defined as Equation (10):

 Longest-distance nodes
To cooperate with the time window, some long segments may be kept in routes and difficult to replace in local search, which leads to the increase of total route length [27].Longest-distance nodes   are end points of the long segment, defined as Equation (11), where S denotes the current solution, eij or eji denotes a segment in S. ̅ denotes the mean spatial distance of the Voronoi neighbors with Voronoi distance 1.The length of the long segment is more than three times ̅ .

 The smallest route's nodes
In some solutions, small routes serve only one, two, or three customers near the depot, which requires more vehicles [48].Usually, this feature is kept for the suboptimal property.The small route's nodes are defined as   , as Equation (12), where |r| denotes the number of visited customers in the route r:

Removal Algorithm
To increase diversification of the guided search, more nodes are removed from the current solution in this study.Additional nodes   are Voronoi neighbors of defined features, just as: Therefore, removed nodes  are spatial-temporal features nodes and their Voronoi neighbors   , as in Equation ( 14).The removal algorithm sequentially ejects all nodes belonging to  and leaves a partial solution for the insertion.

Insertion Algorithm
This algorithm inserts all removed nodes into new places in the partial solution.First, all removed nodes are randomly pushed into a queue L.Then, the front node of L is popped out and inserted at the best place before or after its Voronoi neighbors to keep spatial-temporal proximity, which minimizes insertion cost as Equation ( 3), under the limitation of total duration time, time window, and vehicle capacity.It is of note that the end time windows should not be violated.If there is no proper place, a new route is created, and the front node of L is popped out and initialized as the first-served customer.The insertion operation is repeated until all removed nodes are located in proper places.Due to the best choice at each insertion, the final solution is up to the node sequence in the queue.In this article, the sequence is randomly arranged 50 times to obtain different rebuild solutions.Finally, the obtained result of the permutation that minimizes the total objective is accepted as the final solution.

An Extension of the Solving Algorithm for MDVRPTW
MDVRPTW is a variant of VRPTW, which has more than one depot to serve geographically scattered customers.Compared with the VRPTW, MDVRPTW has an additional problem of how to allocate customers to the right depot.Voronoi partition that divides a study area into many regions with given objects is another typical function of the Voronoi diagram.For the effectiveness, it has been widely used in service area analysis [32].By using Voronoi partition characteristics of the depots' spatial-temporal Voronoi diagram, we assign customers to be served by the depot of its located Voronoi regions.By that, the presented STVDH is easy to be extended to deal with a large-scale MDVRPTW by decomposing it to many smaller VRPTWs.
Figure 7 illustrates the workflow of the STVDH's adaption for the MDVRPTW.The depots' spatial-temporal Voronoi diagram is built to allocate customers to the located depot so that the MDVRPTW is divided into several VRPTWs, each of which has a depot and many assigned customers.Routes for each VRPTW are first initialized using the construction algorithm in Section 4.1 and then separately improved by the spatial-temporal Voronoi diagram-based local search heuristic in Sections 4.2 and 4.3.To overcome the border effect of Voronoi diagram, an additional improvement between depots is used to move or exchange nodes between routes served by different depots.This is done by using 1-0 and 1-1 exchanges to move nodes between routes served by different depots after the spatial-temporal features guided search.Details of this process are referred to TU et al. [49].The stopping condition is the same as the STDVH.Finally, the best found solution is reported.

Experiment and Comparison
To assess the performance of the STVDH, experiments with benchmark problems and large-scale real-world VRPTW and MDVRPTW problems in Shanghai, China were conducted.The STDVH algorithm was implemented with C++, running on a Windows 7 64-bit system with an Intel Core i7-3.4Gprocessor and 16 GB memory.This section reports test datasets, parameter tuning, computing results, and comparisons with other heuristic algorithms.Computing times of compared algorithms are transformed to the computing platform of STVDH with Dongarra factor [52].

Test VRPTW and MDVRPTW Problem Dataset
Two types of VRPTW and MDVRPTW problem datasets are tested.The first dataset includes the VRPTW benchmarks of Solomon [15] and Gehring and Homberger [53].The benchmarks of Solomon have 56 problems with 100 customers.They are divided into six groups with eight to twelve problems, named R1, R2, C1, C2, RC1, and RC2.In the Euclidean plane, customers are randomly distributed (R1 and R2), cluster distributed (C1 and C2) or semi-cluster distributed (RC1 and RC2).Every instance has a single depot within the spatial domain of the customers.The travel time between nodes is equal to the corresponding Euclidean distance.R1, C1, and RC1 have a short time scheduling horizon, whereas R2, C2 and RC2 have a long time scheduling horizon.Similar to Solomon's VRPTW problems, the VRPTW benchmarks of Gehring and Homberger [53] have 300 VRPTW instances with 200, 400, 600, 800 or 1000 customers [54].Customers have the same distributions.
The second dataset is a large-scale VRPTW and MDVRPTW problem dataset in Shanghai, China [55].Five VRPTW instances (Sh1a-Sh5a) and 5 MDVRPTW instances (Sh1b-Sh5b) with 2000 to 10,000 customers and one to six depots are generated to simulate the daily parcel delivery service of a logistics company in Shanghai, China.Customers are randomly selected from the commercial points of interest (POI) using a professional digital navigation map to obtain real-world spatial locations.The time window of each customer is a random value in min bounded by 30 and 60, and the service duration is a random value from 1 to 4 min.Demand is a random value in the interval [1, 100].Depots are located in professional logistic districts in the city.All vehicles have the same capacity, 2000, and maximum route duration, 480 min.The mean travel speed is set to be 45 km/h.

Parameter Tuning
There are four parameters to be tuned in the presented STVDH algorithm as shown in Figure 3.Following Coy's calibration approach [56], a preliminary experiment is conducted on an instance in the problem dataset to find the best parameter settings.The tuning process initially identifies the first two parameters related to the stopping condition and then the last two parameters related to the Voronoi distance decay strategy.The first two parameters, the maximum number of iterations, Nmax, and the maximum number of iterations to start the spatial-temporal features guided local search, Zmax, control the exit of the STVDH.After intensive experiments with different stopping conditions, Nmax should be set to 5000N to converge on a stable high-quality solution, where N is the number of customers.Zmax is set to 100N to start the spatial-temporal feature guided search.
The last two parameters control the Voronoi diagram-based speedup strategy.Experiments on some selected VRPTW instances are conducted to set parameter values for one type of VRPTW dataset.Taking the second VRPTW dataset as an example, we solved Sh2a with all combinations of the last two parameters.Finally, we set   = 3 and = 2 for the best performance in the experiment.Without loss of generality, identical parameter values were set for each instance.The impact of the Voronoi diagram will be discussed in Section 6.A summary of parameter settings is listed in Table 1.Taking Sh2a as an example, Figure 8 displays the convergence of the number of routes and the total length of routes with the parameter settings used in solving Sh2a.

The Results of the Solomon [15] and Gehring and Homberger [53] VRPTW Benchmarks
The parameter setting for these VRPTW benchmarks is as: Nmax = 5000N, Zmax = 100N,   = 2 and  = 1.5.For each problem, the STDVH was run 10 times.The total objective, total routes length, and computing time were recorded.Results of three state-of-the-art methods, including the arc-guided evolutionary algorithm (AGEA) [39], the hybrid genetic algorithm with adaptive diversity management (HGASDC) [40], and the local search based heuristics VRPEJ [57], which limits the searching space with the kth nearest neighbor strategy, are compared with the STVDH.Tables 2 and 3 present the comparison of the results on the Solomom and Gehring and Homberger VRPTW benchmarks respectively.Each entry reports the average of computing results of a group (average computing time | the average total routes length).They indicate that the presented STVDH is able to achieve high quality solutions of the VRPTWs with no more than 1000 customers in a relatively short computing time.In terms of the number of routes, the STVDH achieves the fewest routes in both the Solomon (total of 405 routes) and Gehring and Homberger (total of 10296 routes) VRPTW benchmarks.In terms of the total route length, the STVDH performs better than the VRPEJ for all size VRPTW instances, only inferiorly to AGEA on Solomon benchmarks, but superiorly to AGEA on the Gehring and Homberger benchmarks.However, the STVDH's results are still worse than that of the HGSADC.For the computing efficiency, the STVDH consumes the least computing time in all six size VRPTW instances (100~1000).Therefore, the proposed STDVH exhibits a good performance on the VRPTW benchmarks of Solomon and Gehring and Homberger.

The Results of Large Scale VRPTW and MDVRPTW Problem in Shanghai, China
The parameter setting for this large scale VRPTW and MDVRPTW problem is as follow: Nmax = 5000N, Zmax = 100N,   =3 and  = 2.For each problem, the STVDH algorithm was also run 10 times.Total travel distance, number of used vehicles, and computing time were recorded.The obtained results are summarized in Table 4 (https://github.com/spatialsmart/VRPTW/tree/master/Shanghai/Solution).
As Table 4 indicates, the STVDH algorithm reports the best found solution for large-scale VRPTWs up to 10,000 customers in 139.5 min (about 2.32 hours).As the number of customers increases from 2000 to 10,000, the computing time of the STVDH increases from 15.7 min to 120.8 min, whereas for the MDVRPTW, the time increases from 16.1 min to 139.5 min.Compared with the VRPTW, more computing efforts are required for the MDVRPTW due to additional customer allocation between depots.In terms of solution quality, the STVDH utilizes an average of 847.1 total routes that travelled 30,184.725km to serve all customers in the five VRPTW instances.For the five MDVRPTW instances, 906.4 routes that travelled 26,353.572km are required to satisfy customer needs.It should be noted that the performance of the STVDH is robust, as indicated by the gap between Savg and Sbest (−1.60% for the VRPTW and −1.73% for the MDVRPTW).Therefore, the proposed STVDH provides promising results for large-scale VRPTW and VRPTW within reasonable computing times.
Taking problem Sh1a as an example, Figure 9 displays the best found solution and the spatial-temporal vehicle route with ArcScene.It demonstrates the best route compromise between spatial proximity and time constraints.The same vehicle may not yet serve spatially-near customers as Figure 9b shows.
To assess the solution quality of the obtained results, we compare them with solutions reported by two heuristic algorithms.One solution is obtained by solving the simulated instances with the network analysis module in ArcGIS TM 10.2, which uses a tabu heuristics algorithm to solve complex real-world VRPTW problems.Another compared solution is calculated using VRPEJ [55].Other metaheuristics such as AGEA and HGSDAC are not compared due to unavailability.Each algorithm was run 10 times to find the best result.

Discussion
Taking the large scale VRPTW and MDVRPTW in Shanghai, China as an example, this section discusses the impact of the spatial-temporal Voronoi diagram and spatial-temporal proximity behind the best found results.

Impact of the Spatial-Temporal Voronoi Diagram
To evaluate the impact of the Voronoi diagram, this study investigated the balance of solution quality and computing time with different values for parameters   and  .As Figure 10 illustrates, considering the spatial-temporal Voronoi neighbors, an increase in   improves solution quality but decreases computing efficiency because more spatial-temporal Voronoi neighbors are involved in local search procedures.With regard to the effect of the Voronoi distance decay strategy, as Figure 10a shows, both the slow decaying speed ( = 1) and the fast decaying speed ( = 3) generate worse results.However, as Figure 10b shows, as a fast decaying speed requires less evaluation on far Voronoi neighbors, the STVDH consumes less computing efforts as the parameter  increases.

Spatial-Temporal Proximity Analysis on the Best Found Results
To understand spatial-temporal proximity behind the obtained results, we conducted an analysis on the Voronoi distance between consecutively-visited customers in the route of best found solutions in Section 5.2. Figure 11 displays the average percentage of Voronoi distance in the best found solutions.Two typical features are indicated below.
 The number of larger Voronoi distances is only a small proportion of the best found results.As Figure 11 displays, the percentages of Voronoi distances greater than three in the VRPTW and MDRPTW solutions are only 1.77% and 1.58%, respectively.Such a distribution agrees with the setting of parameter   (   = 3) in Section 5.2.Therefore, the spatial-temporal Voronoi neighborhood is very typical in the found solution. The percentage decreases sharply as the Voronoi distance increases.For the large-scale VRPTW dataset (Sh1a-Sh5a), the percentages for Voronoi distances 1, 2, 3, and >= 4 are 59.44%, 30.13%, 8.66%, and 1.77%, respectively, which is similar to the distribution in the MDVRPTW solution.This verifies that there is a spatial-temporal local compact structure in the routes of best found solutions.Compared with the theoretical searching intensity (as Equation (7) in Section 4.2.1 where   = 3,  = 2) in the Voronoi distance decay strategy, these percentages systematically slightly shift to small Voronoi distances.This result demonstrates the effectiveness of the Voronoi distance-decay strategy, which searches more on near neighbors in the local search but still spends necessary efforts on far neighbors.
In summary, analysis of best found results demonstrates reasons why the spatial-temporal Voronoi diagram is effective in the STVDH algorithm.

Conclusions
Vehicle routing designs the least costly routes to satisfy the geographical distribution of human needs.Embedded in the GIS environment, it not only benefits transportation decision-making for both public and private sectors, but also enriches the value of geoinformation with the assistance of location-based service.Inspired by spatial-temporal proximity, this article presents a novel spatial-temporal Voronoi diagram-based heuristic approach to solve large-scale VRPTWs quickly.The derived spatial-temporal Voronoi neighborhood measures proximity considering both spatial and temporal issues.The used Voronoi neighbors limits searching space in the construction algorithm, local search, and spatial-temporal feature-guided search.Moreover, the presented Voronoi distance-decay strategy assigns more searching efforts on spatial-temporal near neighbors.
Experiments on the VRPTW benchmarks and large-scale VRPTW instances in Shanghai, China with 2000-10,000 customers have demonstrated the good performance on VRPTW and MDVRPTW problems of different sizes.The obtained results indicate that the STVDH presented in this study can provide high quality solutions for large-scale VRPTWs in a reasonable time with the help of a spatial-temporal Voronoi diagram.It also verifies that the spatial-temporal proximity is very typical in best found solutions as 98% of segments had a Voronoi distance less than three.
The main contributions of this article are twofold.First, considering both time windows and spatial locations, the spatial-temporal Voronoi diagram is proposed to accelerate the solving process for VROs.In contrast to the spatial proximity based speeding up strategies [25][26][27], it makes use of spatial-temporal Voronoi neighbor to further reduce computing effort for local search based heuristics.The idea behind the presented approach that makes uses of spatial principles to accelerate complex optimizing processes can be extended to facilitate other space-related decision-making problems such as near real-time emergency response and large-scale facility location.Second, a spatial-temporal Voronoi diagram-based heuristic was developed to solve large-scale VRPTWs.This novel approach exhibits good performance on super large-scale VRPTWs up to 10,000 customers, which can cope with challenges from many complex real-world transportation and logistics applications.
Among future developments that we intend to undertake, we plan to integrate the presented approach with outdoor/indoor ubiquitous positioning and friendly navigating technologies for vehicles with cloud GIS.The presented effective approach will be used to upgrade the network analyst module in SDSSs to provide a more flexible transportation service for daily large-scale logistics and distribution activities.It will further enhance the spatial intelligence of modern transportation and logistics applications.

Figure 2 .
Figure 2. The spatial-temporal Voronoi diagram and the spatial-temporal Voronoi distance.

Figure 3 .
Figure 3.The framework of the spatial-temporal Voronoi diagram-based heuristic algorithm.

9 )Figure 6 .
Figure 6.Wait time and time warp at a customer.

Figure 7 .
Figure 7.The adaption of the STVDH for the MDVRPTW.

Figure 8 .
Figure 8.The convergence of the total number of routes and the total routes length (Sh2a).

Instance:Figure 9 .
Figure 9.The result for large scale VRPTW problem Sh1a.(a) the final best solution; and (b) a spatial-temporal route.

Figure 10 .
Figure 10.The impact of the spatial-temporal Voronoi diagram on the STVDH.(a) Deviation to the best solution; and (b) computing time.

Figure 11 .
Figure 11.Average percentage of the Voronoi distance in the best found solutions.
The route duration time   +1 −   should be no more than the maximum duration Tmax; (2) The total served demand ∑    =1

Table 1 .
The third parameter,   , indicates the maximum searchable Voronoi neighbors.As the value of   increases, local search will evaluate more neighbors, which requires more computing efforts.According to the distribution of the Voronoi neighbors,   should be no more than 5 [49].The last parameter, , determines the distance decaying speed.A high value indicates a tendency for fast decay.The proper value of  is within the bound [1,4].The value space of parameters are summarized in Table 1.Parameter settings of the STDVH algorithm.
N is the number of customers.
Instance: group name.N: number of customers.CNV: cumulative number of vehicles.CND: cumulative travel distance.T: time reported by metaheuristics.T2: time on the computing platform of the STVDH.

Table 3 .
Cont.Instance: group name.N: number of customers.CNV: cumulative number of vehicles.CND: cumulative travel distance.T: time reported by metaheuristics.T2: time on the computing platform of the STVDH.

Table 4 .
Results from the large-scale VRPTW and MDVRPTW instances in Shanghai, China.

Table 5
compares the obtained results.It indicates that the STVDH presented in this study outperforms both ArcGIS and VRPEJ.In terms of efficiency, comparison with ArcGIS (total of 2141.1 min) and VRPEJ (total of 2355.3 min) indicates that the STVDH costs the least computing effort (total of 703 min).In terms of solution quality, the STVDH requires 1729 routes travelling 55,990.55km to serve all customers in 10 instances.ArcGIS requires more vehicles (1764 routes) that travel a greater distance (62,047.809km) to provide the same service.The VRPEJ, which utilizes a spatially kth nearest neighbors strategy, requires the most routes (1893 routes) that travelled 63,659.365km.Hence, it confirms that the spatial-temporal Voronoi neighbor strategy in the presented STVDH is much better than the spatial neighbor strategy for local search to solve large scale VRPTWs.

Table 5 .
Comparison of the results of the STVDH with other heuristic algorithms.

.365 Gap to STVDH 10.82% 13.70% Total Computing Time (/min) 703.0 2141.1 2355.3 Instance
: problem name.Type: VRPTW or MDVRPTW; N: number of customers.M: number of depots.STVDH: best of 10 runs (number of routes| total routes length).ArcGIS: best result of 10 runs (number of routes| total routes length).VRPEJ: best of 10 runs (number of routes| total routes length).CNV: cumulative number of vehicles.CND: cumulative travel distance.