A Multi-Stage Algorithm for a Capacitated Vehicle Routing Problem with Time Constraints

The Vehicle Routing Problem (VRP) is one of the most optimized tasks studied and it is implemented in a huge variety of industrial applications. The objective is to design a set of minimum cost paths for each vehicle in order to serve a given set of customers. Our attention is focused on a variant of VRP, the capacitated vehicle routing problem when applied to natural gas distribution networks. Managing natural gas distribution networks includes facing a variety of decisions ranging from human resources and material resources to facilities, infrastructures, and carriers. Despite the numerous papers available on vehicle routing problem, there are only a few that study and analyze the problems occurring in capillary distribution operations such as those found in a metropolitan area. Therefore, this work introduces a new algorithm based on the Saving Algorithm heuristic approach which aims to solve a Capacitated Vehicle Routing Problem with time and distance constraints. This joint algorithm minimizes the transportation costs and maximizes the workload according to customer demand within the constraints of a time window. Results from a real case study in a natural gas distribution network demonstrates the effectiveness of the approach.


Introduction
In the current competitive environment, manufacturing and service companies are facing different challenges: global markets, increased competition, customers behavior, accelerated lead and delivery times as well as process technology capabilities, which are contributing to an increase in uncertainty and variability [1,2]. Therefore, choosing the shortest path is mandatory to minimize costs and time and allow companies to focus on significant sustainability issues, which have become a core strategic imperative pushing companies towards identifying new opportunities [3].
The Vehicle Routing Problem (VRP) is one of the most optimized tasks studied and it is implemented in a huge variety of industrial applications. The VRP is defined as the optimal delivery route from one or multi depots to customers, taking into account distance, capacity and time constrains. VRP is used when a series of routes must be identified to optimize the distribution of goods to various cities using the most economical method with each vehicle leaving and returning from the central depot (or from several depots). The objective is to design a set of minimum cost paths for each vehicle in order to serve a given set of customers [4].
The VRP manages the design and optimization of vehicle fleets' routes and is applied to distribution, logistics and supply chain issues. The VRP presents several variants on mathematical modelling such as capacitated VRP, multi-depot VRP, periodic VRP, split delivery VRP, stochastic VRP, VRP with backhauls, VRP with pick-up and delivery and VRP with time windows. Additionally, under • Designing of new NGDN or expansion of existing ones [9]; • Routing a fleet of vehicles from a depot to service a set of customers that are geographically dispersed [10]; • Determining the pipelines in the network, the compressor stations and their capacities [11]; • Minimizing the transportation costs and maximizing the transported volume [12]; • Optimizing gas network inventories in the face of system uncertainties [13]; • Synchronizing the minimum risk loss and total cost of natural gas pipeline networks [14].
Despite the numerous papers available on VRP, only a few study and analyze the problems occurring in capillary distribution operations such as those found in a metropolitan area. However, there are some exceptions such as [15], where a set of field service operations are considered in order to manage the assets of a NGDN and maintain high levels of customer satisfaction and gas provider profitability, particularly regarding activation, deactivation and management of utility meters with respect to time, distance and capacity, and where an operational VRP characterized by the presence of both mandatory and not mandatory customers, and partially overlapping time windows were presented. The authors adopted a VRP with time windows (VRP-TW) combined with Simulated Annealing (SA), which allows for real constraints optimization in reasonable computational time and represents an alternative approach with respect to existing contributions to the considered subject. This algorithm interface operates as a decision support system tool that allows decision makers to both visualize and to modify the routes while allowing interesting optimization capabilities.
The implementation of the established algorithm to a real case study resulted in a far more complex system, where both running time and usability were not effective, resulting in a pilot approach which did not allow the company to adapt the changes needed or take into consideration the additional paths and real time constraints.
Accordingly, this research proposes a different algorithm which takes into account these drawbacks and suggests a new solution, that is, simplifying the problem by reducing the number of constraints and varying some variables. Hence, a true comparison of the two studies is not possible.
This work introduces a new algorithm based on the Saving Algorithm (SA) heuristic approach which aims to solve a Capacitated Vehicle Routing Problem (CVRP) with time and distance constraints. This joint algorithm minimizes the transportation costs and maximizes the workload according to customer demand within the constraints of a time window and it allows for increased effectiveness in real-time scenarios.
The paper is organized as follows: in Section 2, the problem is described. Section 3 formally presents the CVRP and the joined algorithms. Experimentation and results are shown in Section 4 and finally, in Section 5 some conclusions are drawn.

Problem Definition
This problem derives from a real case study concerning a small medium service enterprise, which manages domestic gas utilities in the metropolitan area of Genoa (Italy). The problem it faces is quite complex and is related to its field service operations of utility networks in the metropolitan area. The heuristic approach proposed in this paper is based on the CVRP with time constraints, and it exploits the well-known cluster-first-route-second problem decomposition technique to minimize the distance within the cluster by generalized assignment problem.
The study deals with the planning of service visits to a set of customers, performed by a set of company operators, each using a vehicle from the company fleet.
The service company provides two types of operations: Activation and Deactivation. Activation is when a new customer opens an account and deactivation is when a customer closes an account. Activations and deactivations are fixed for a specific day by the gas company provider who arranges an appointment with the customers within a specific time window.
The fleet is a set of vehicles each associated with a single driver whose workload should not exceed two-time based restrictions: • A maximum daily workload (visits per day); • A maximum workload for the time window (time allotted per visit).
Hence, the problem consists of finding a method to optimize company logistic resources, while maintaining a maximum level of customer satisfaction in the most economical way.
This algorithm improves the distribution efficiency and increases customer service quality using the aforementioned heuristics, the Saving Algorithm and 2-Opt.
The algorithm assigns a sequence of interventions to the various operators, which are carried out during the working day (divided into time slots) therefore minimizing the total distance traveled in compliance with the following operational constraints: The algorithm calculates the distances between the various interventions and provides a solution in terms of:

•
A list of interventions for each operator consisting of an account number ID, location and time window for each customer is assigned; • The total time of each intervention (including transport time).

Proposed Algorithm and Formulation
The CVRP is one of the most examined VRP applications. In CVRP, the objective is to minimize the costs of the routes (similar to the generic VRP) but at the same time not exceed the load capacity of the vehicle assigned to those paths [16].
In CVRP, customer demands, locations and depots are predetermined. The depots store all the products or services ordered by customers and each vehicle leaves and returns to the same depot. Therefore, a solution is feasible if it satisfies these aspects. Note that a route must be assigned to only one vehicle.

Basic Components of CVRP
Let us give a complete graph G = (I, J, K), the notation used and the mathematical model are as follows: customer index k route index Parameters: N number of vehicles C ij distance between point i and j, i, j ∈ I ∪ J d j demand operation of customer j (in terms of time) Q k time limit for operator k L j service time at customer j Decision variables: 1, i f i preceeds j on route k 0, otherwise z ijk = 1, i f customer j is matched to a depot i 0, otherwise y j = cumulative travel and service time up through customer j, j ∈ J U lk = auxiliary variable for sub − tour elimination constrains in route k The objective function is to minimize the total distance of all vehicles given by: A single route has to be given to each customer: Each vehicle has a time capacity constraint: The following equation defines the sub-tour elimination constraint: While Equation (5) expresses the flow conservation: Each route is covered by one vehicle: A customer is assigned to a depot only if there is a route from the aforementioned depot going to that customer: Equations (8) and (9) give the binary requirements on the decision variable: Finally, Equation (10) defines the positive values of the auxiliary variable: To summarize, the objective function (1) minimizes the total distance of all vehicles with respect to some constraints: (i) each customer is satisfied with only one route (2); (ii) the total daily time capacity of each operator (3); (iii) the sub-tour avoidance (4); (iv) the flow limits (5) and (v) route to be served (6).
CVRP is NP-hard, therefore an efficient algorithm for solving the problem to optimality is generally unavailable with respect to real scenarios. To deal with the problem efficiently and effectively, a heuristic approach such as the Saving Algorithm and 2-Opt is applied.

Proposed Algorithm
The proposed algorithm starts from the classic problem of CVRP with a time window constraint to define routes which minimize the costs and maximize the number of activations or deactivations per worker with respect to (i) a maximum daily workload and (ii) a maximum workload for each time window. Specifically, the NGDN considered in this paper provides a service, therefore there is no load capacity for the vehicles, but there is a time constraint for the number of operations per worker.
The proposed algorithm includes five phases, shown in Figure 1, that are detailed in the following.

Phase 1: Initialization
The first step of the algorithm is based on the initialization phase, in which all data on operation requests are collected and organized with respect to: • Account number ID of the operation; • Type of operation Location of the operation in terms of latitude and longitude; • Time window for the requested intervention.

Phase 2: Solving Distance and Saving Matrix
The second phase of the algorithm is based on the Saving Matrix calculation, in which the distances between the depot and the customer locations are defined. At this point, it is important to underline that the case considered is based on a single depot, but equations could be extended to multi-depot problems. The Saving Matrix is calculated using the Saving Algorithm heuristic. The Saving Algorithm is a well-known heuristic for solving CVRP [17]. In the Saving Algorithm, distance and travel time savings are calculated when two routes merge together. Based on the saving order, routes are merged together sequentially and a better, more feasible solution will be generated [18].
As Figure 2 shows, the basic idea of this algorithm is to combine routes together based on the saving matrix.  In the initial step, all nodes are connected to the depot (v 0 ) and N routes are generated as an initial solution (Phase 2.1). As shown in Figure 4, two different routes and policies can be generated and calculated through the Haversine formulation. It determines the great-circle distance between two points on a sphere given their longitudes and latitudes:

•
Policy a: the vehicle starts from the depot, carries out the operation in v 1 and then returns to the depot, heads to v 2 , carries out the operation and returns to the depot. The cost of this route is: • Policy b: the vehicle starts from the depot, carries out the operation in v 1 and heads on to v 2 , carries out the operation and returns to the depot. The cost of this route is: Policy b allows a saving of: We further assume that: Therefore, we can calculate the distances between all nodes of V and the savings related to them (Phase 2.2). The saving matrix is generated by calculating the saving distance and therefore the travel time for merging any two possible routes (i, j) together.

Phase 3: Clustering
In the clustering phase, each route is generated by grouping operations which allow maximizing the saving considering the time constraint. Based on the saving matrix, a list L(i, j) is created and sorted in descending order by merging savings. With list L(i, j), road i and road j are selected to be tested with the feasibility of mergence. If the mergence is feasible, these two roads are merged together and a new merged route is created and the saving matrix is updated with new routes. If no feasible mergence is found or it exceeds the given calculation time, the process is terminated.
In the example described above (Figure 4), we need to consider: T 01 + T 1 TW 1 (17) Otherwise, if the time constrain expressed by Equation (16) is not verified, two routes are created and assigned to different operators (Policy a): The output of this phase is the definition of the first solution, which can be improved through the next phase.

Phase 4: Route Sequencing
The first solution found in the second and third phase can be improved by route sequencing (maximizing savings). This means that for each route, we need to:

1.
Evaluate the initial solution cost, and 2.
Improve the current solution by a 2-Opt procedure.
This paper applies a 2-Opt local search algorithm. In each cluster, a new route will be generated by solving the 2-Opt algorithm, subdividing the original problem into sub-problems that are faced in a sequence of successive phases. In a 2-Opt algorithm, two tours are neighbors if one can be obtained from the other by deleting two edges, reversing one of the resulting two paths, and reconnecting [19]. The 2-Opt neighbors of the tour are pair-wise part exchanges within the current solution sequence [20]. An exhaustive search visits all n! tours, while the total number of tours visited by 2-Opt, including the original data set, is given by the following formula: The adoption of the 2-Opt procedure allows the scalability to be optimized, as a 'linear' increase of computational time is expected rather than a factorial. The 2-Opt method improves the local search ability of the algorithm and its overall performance. The output of this phase is the definition of a second, improved solution.

Phase 5: Final Optimization
The final optimization phase aims at improving the overall solution starting from the second improved solution, attempting to maximize the available time of an operator taking into account the corresponding time window. In fact, after sequence optimization provided by the 2-Opt algorithm, there is the possibility that a route can be operated with less time than the time identified by the initial solution. Therefore, an additional operation may be added to the optimized route. Taking a look at the previous simple example ( Figure 4) where: An additional operation can be manually introduced in order to fill AT k . Therefore, a cluster is composed of the route v 0 -v 1 -v 2 -v 3 -v 0 ( Figure 5) if: where: Time necessary to go from node v 2 and v 3 (T 23 ) expressed in (min) Time of the operation performed in v 3 (T 3 ) expressed in (min) Figure 5. Improved path.

A Real-Case Scenario Application of the Proposed Algorithm
In this section, a real-case scenario application of the proposed algorithm is provided. The problem is composed of 54 activations and 51 deactivations (for a total of 105 service operations) which need to be assigned to workers; each operation is described with an account number ID, the location and the time in which the operation has to be completed, an example of the sample is reported in Table 2.
The algorithm was implemented on a Windows 10 × 64, Intel Core i7 CPU 2.40 GHZ, 4 GB RAM. After all the data from the activations and deactivations have been collected (Phase 1), the algorithm provides a solution taking into account all constraints and then allocates operations to the workers, solving the aforementioned Phases 2, 3, 4 as shown in Table 3. The solution is calculated thanks to the definition of distance and saving matrix for each pair of destinations. Table 3 shows the account number IDs assigned to an operator in time window 1. The two last columns represent the total time of intervention (including transport time) and the not allocated time for each operator, respectively. Table 3 shows that the operator 1 carries out the interventions for IDs 70, 17 and 72 with a total intervention time of 77.2 min, consisting of 60 min of operating time (2 deactivations and 1 activation) and 17 min of travel time.
Then, Table 4 reports the workload that has been assigned to each operator throughout the day and information related to the operating time (including travel time) for the whole day. Finally, in order to optimize routes (Phase 5), Table 5 shows the distance from one location to all those located in the same time window. According to Table 5 the position of intervention ID 12 would  be close to interventions ID 2, 1, 4, 3, 62 and ID 18, quite close to interventions ID 70, 17, 72, 64, 13, 56, 57, 55, and ID 5 and distant from the remaining ones. Through this color scale the user can see at a simple glance any useful movements and make the appropriate changes in order to optimize the routes once again. In the example presented here, the intervention ID 12 is moved from operator 9 to operator 8 respecting the AT 8 . Indeed, AT 8 before the Phase 5 was 49.1 min, by adding intervention ID 12, it totals106 min, which is less than the maximum available time established at 120 min. Before Phase 5 Distance from intervention 12 to all those located in the TW1 After Phase 5 Green: distance less than 3 km; Yellow: distance between 3 and 6 km; Red: distance greater than 6 km. Table 6 and Figure 6 show the information relating to the itineraries defined after Phase 5 for Operator 8. It contains information on each intervention assigned to him, such as the type and time window in which it is located and the visual path.
Finally, a comparison between the algorithm provided by [13] and the one proposed in this paper is shown. It takes into account:

•
The running time necessary for the algorithm to provide the solution; • The imbalance; • The total workload for a set of operators.   62  TW1  18  TW1  12  TW1  67  TW2  6  TW2  61  TW3  10  TW3  25  TW4  29  TW5  32  TW6  78  TW6 The comparison deals with the problem, which as previously described is composed of 54 activations and 51 deactivations (reported in Table 2). Table 7 reports the results of the comparison. In particularly, there is an interesting decrease in running time (the computation time in minutes) as the total running time of the previous algorithm is 8.5 min while the proposed algorithm is completed in 3 min. Clearly, the difference between the two running times depends on the fact that the first algorithm takes into account more constraints with respect to the second solution, such as:

•
The load capacity (difficulty level of interventions) for the vehicle: the vehicle has a strict capacity, there is a maximum daily score per vehicle; • The workload for a worker: it takes into account the maximum daily score per operator, the maximum score per time window and additionally the maximum score per type of operation.
In light of the explanation above, it is important to point out that the solution provided by the proposed algorithm is clearly unfavorable. As Table 7 shows, the previous algorithm requires 9 operators and a total workload of 2673.2 while the current algorithm proposes 13 operators with a total workload of 2809.5. Furthermore, the imbalance increases from 182.9 min per operator per day in the previous case to 263.8 min in the case of the proposed algorithm. Table 7. Comparison between the previous and the current work. The analyzed operations (51 activations and 54 deactivations) represent the real data set analyzed and managed by the pilot company. In terms of time performance, by using a heuristic approach, a 'linear' increase in computational time would be expected in the case of a greater number of operations to be managed. Here, it is important to note that the proposed approach stems from the evident limitations of a previous math-heuristic approach that the authors developed to address the same problem. After almost a year of the pilot phase, the company showed difficulties in setting up and maintaining the set of constraints, as well as in adding or removing constraints and in managing frequent variations and unforeseen events occurring during day-by-day operations. Hence, the proposed approach is less 'optimal' (by using a set of heuristics) but allows the company to iteratively re-allocate operators to modified conditions in less time and in an easier way, thus, providing a more robust and effective solution to the real-time scenario. The benefits of the proposed approach are not just about pure performance but are related to its usability and applicability (less computational time, less constraints to be managed by adopting an appropriate simplification, better reactivity) and the ability to dynamically generate solution routes.

Number of Operators
To conclude, the previous algorithm performs well but it is less effective in a real-case scenario, as the case study pilot phase demonstrated. Clearly the performance of our approach is lesser but its ability to manage the real case scenario is greatly improved, showing that more interactive and easy-to-parametrize "optimization" algorithms could be more effective depending on the application context.

Conclusions
The proposed algorithm based on the Saving Algorithm heuristic approach which aims to solve a Capacitated Vehicle Routing Problem (CVRP) with time and distance constraints was developed. This paper starts with the results of a practical application to a real scenario as reported by [13], where an operational VRP characterized by the presence of both mandatory and not mandatory customers and partially overlapping time windows was presented. The authors adopted a VRP with time windows (VRP-TW) combined with Simulated Annealing (SA), which allows for real constraints optimization in reasonable computational time and represents an alternative approach with respect to existing contributions.
The implementation of the previous algorithm to the real scenario resulted in an overly complex system setting to manage all constraints and both the running time and usability were not effective. Furthermore, the highly-sophisticated optimization model was difficult to modify and maintain by the company, for example, to implement new additional paths and constraints. Because of these limitations, the presented paper proposes a reasonable reduction in complexity (accepted by the service company) in order to design and implement a multi-stage algorithm to solve the reduced problem. The heuristic approach proposed in this paper, based on the CVRP with time constraints, effectively tackles the complex problem regarding the field service operations of utility networks in a metropolitan area, by exploiting the well-known "cluster-first-route-second" problem decomposition technique to minimize the distance within the cluster by generalized assignment problem. In each cluster, a new route is generated by solving the 2-Opt, subdividing the original problem into sub-problems that are faced in a sequence of successive phases. This algorithm minimizes the transportation costs and maximizes the workload according to customer demand within the constraints of time window.
Advantages, documented by the results provided from the practical application of this algorithm, are:

•
The results demonstrated the efficiency of the approach in tackling a complex optimization task by providing highly pragmatic and realistic time-responsive solutions; • Flexibility since it can be applied to different routing problems not only to NGDNs; • Increased "adaptability" thanks to an implemented approach of optimizing routes (greater parametrization freedom), taking into account the different constraints, depending on the specific problem under study.
As a future study, this new approach can be tested and compared with some additional comparative numerical results from other relevant research works in order to show the benefit of adopting this proposed framework. Additionally, future developments could be the real-time connection of operators to this heuristic reactive approach in order to know the state and location of operators and to improve their allocation with respect to service visits inserted in the time window.

Conflicts of Interest:
The authors declare no conflict of interest.