Mathematical models for the periodic vehicle routing problem with time windows and time spread constraints

Article history: Received: 10 December 2019 Accepted: 1 June 2020 Available Online: 10 September 2020 The periodic vehicle routing problem (PVRP) is an extension of the well-known vehicle routing problem. In this paper, the PVRP with time windows and time spread constraints (PVRP-TWTS) is addressed, which arises in the high-value shipment transportation area. In the PVRP-TWTS, period-specific demands of the customers must be delivered by a fleet of heterogeneous capacitated vehicles over the several planning periods. Additionally, the arrival times to a customer should be irregular within its time window over the planning periods, and the waiting time is not allowed for the vehicles due to the security concerns. This study, proposes novel mixed-integer linear programming (MILP) and constraint programming (CP) models for the PVRP-TWTS. Furthermore, we develop several valid inequalities to strengthen the proposed MILP and CP models as well as a lower bound. Even though CP has successful applications for various optimization problems, it is still not as well-known as MILP in the operations research field. This study aims to utilize the effectiveness of CP in solving the PVRP-TWTS. This study presents a CP model for PVRP-TWTS for the first time in the literature to the best of our knowledge. Having a comparison of the CP and MILP models can help in providing a baseline for the problem. We evaluate the performance of the proposed MILP and CP models by modifying the well-known benchmark set from the literature. The extensive computational results show that the CP model performs much better than the MILP model in terms of the solution quality.


Introduction
The periodic vehicle routing problem (PVRP) is an extension of the standard vehicle routing problem (VRP), which plans the vehicle routes over the several planning periods. In each period, routes are constructed for a fleet of capacitated vehicles such that: each vehicle starts and ends its tour at the depot; demands of the customers are delivered, and vehicle capacities are regarded. Then, the goal of the PVRP is to find a set of vehicle routes that minimizes the total travel time (or distance) over the planning horizon. This paper addresses the periodic vehicle routing problem with time windows and time spread constraints (PVRP-TWTS), which arises in high-value shipment transportation. Since the high-value shipment companies generally plan the vehicle routes over several days and most of their customers are visited each day, their transportation processes can be formulated as a PVRP. Besides the aforementioned assumptions and limitations, in the PVRP-TWTS, there are also time windows for the customers, where each customer can only be visited within its time window. Furthermore, due to security concerns, changing arrival times are preferred in high-value shipment transportation. Hence, there are also time spread constraints for the customers. Namely, the arrival times of any two visits to the same customer must vary by a given time constant over the planning periods. Additionally, due to the security reasons, the waiting times are not allowed for the vehicles. The PVRP-TWTS can be seen in real-life during the planning of the daily trips for the loading/unloading operations of the Automated Teller Machines (ATM). Generally, the banks make agreements with the corporations, which are specialized in the transportation of the expensive products, for the loading/unloading operations of these ATMs. The Mathematical models for the periodic vehicle routing problem with time windows… 11 purpose of these Cash in Transit (CIT) companies is to move valuable things like coins, banknotes, or jewels from one location to a different one safely. The customers of these CIT corporations also involve supermarkets and various shopping stores in addition to the banks. Since these CIT companies transport valuable products, they must arrange their daily visits using different routes due to security reasons [1,2]. In this study, we consider the PVRP-TWTS with heterogeneous vehicles and customers with periodspecific demands, which is a novel extension of the problem. We propose a novel mixed-integer linear programming (MILP) model and a novel constraint programming (CP) model for the problem, as exact solution approaches. Furthermore, we develop several valid inequalities to strengthen the proposed MILP and CP models as well as a simple lower bound. Then, we evaluate the performance of the proposed models on the well-known benchmark set of [3] by modifying these instances. In the literature, the CP technique has been successfully applied to various combinatorial optimization problems, and it has become a strong competitor to state-of-the-art mathematical programming techniques such as MILP [4]. However, it is still not as well-known as the MILP in the operations research field. The MILP models have mainly been employed to solve various real-life problems, including the VRP and its variants. Even though some studies have applied CP to solve the VRP [5][6][7][8], MILP models have been widely used for those types of problems in the literature. This study aims to utilize the effectiveness of CP in solving such a complex PVRP-TWTS. To the best of our knowledge, this study presents a CP model for the PVRP with time spread constraints for the first time in the literature. In this paper, we compare the CP and MILP techniques for the PVRP-TWTS. Having a comparison of CP with a MILP model can help to provide a baseline for the problem and increase the usage of CP for similar problems in the literature by showing its effectiveness. The rest of the paper is structured as follows. A literature review is provided in Section 2. Section 3 formally defines the problem. In Section 4, the MILP and CP formulations are presented as well as the valid inequalities and the lower bound. Section 5 presents the computational results to assess the performance of the proposed models. In Section 6, the conclusions are mentioned as well as future research directions.

Literature review
Many exact and heuristic methods have been proposed in the literature for vehicle routing problems such as [9,10]. Three recent comprehensive reviews are provided for the VRP in [11][12][13]. A recent survey on PVRP is also presented by [14]. There are two main approaches in the literature for constructing unpredictable routes: order dissimilarity and arrival time dissimilarity. The order dissimilarity is generally addressed under the peripatetic vehicle routing problems. The m-peripatetic vehicle routing problem (m-PVRP) aims to obtain a set of edge-disjoint routes with a minimum cost over periods such that: each customer is visited only once in a period, and the edge among a pair of customers can be utilized at most once through the periods. The m-PVRP was initially studied by [15], in which several upper and lower bounds were proposed for the problem. Later, the same authors also proposed another lower bounding procedure and an exact method based on polyhedral and column generation approaches for the m-PVRP in [16]. A hard constraint of the m-PVRP is that: no edge is traversed twice during the periods. Later, in [17], this limitation was loosened by proposing a related problem named the k-dissimilar VRP, where the aim is to obtain various VRP solutions such that the similarity among a couple of solutions is under a given level. Note that the k-dissimilar VRP does not forbid the multiple usages of an edge explicitly. Akgün et al. [18] stated that edge dissimilarity does not always result in geographically dispersed routes. In [19], spatially different routes were developed by evaluating the intersection zones between the routes, where Pareto-optimal paths were found by employing a multi-criteria shortest path algorithm. In [20], a biobjective function was employed to minimize the shipping cost and the security risk of transferring highvalue goods. The authors developed a combined risk criterion, which is the weighted sum of the route similarity and the risk of a route. Talarico et al. [21,22] also developed a risk measure related to the quantity of loaded money and the distance traversed by the vehicles. Another strategy to construct unforeseeable routes is to differ the visiting (arrival) time of each customer over the planning periods. In [23], diversification of the visiting times was studied for the case of overnight safety guarding, where the aim is to obtain a set of dissimilar routes with regard to visiting times. Note that they presumed that all customers have the same time window. Constantino et al. [24] studied the problem of gathering the money from parking meters, where the aim is to minimize the total traveling time and keep a similarity measure under a given threshold. They developed a MILP model and a metaheuristic for the problem, where the similarity measure of two tours is expressed as the fraction of tasks that are served in the same period. Yan et al. [25] also studied the arrival time diversification by assuming that the visiting time of a customer must vary by at least minutes from the ℎ previous visiting time. They formulated the problem as an integer multi-commodity network flow problem. To the best of our knowledge, there are two studies [1,2] that consider similar PVRP to our PVRP-TWTS. In the first study of Michallet et al. [1], the PVRP with arrival time dissimilarity was modeled by differing the arrival time at each customer by at least a given time spread constant over the planning periods. In the study of [1], there are also time windows for the customers, and waiting time is not permitted for the vehicles. They developed an iterated local search for the problem, where the infeasible solutions that violate the arrival time dissimilarity and the time window limitations are penalized with a time-dependent penalty function. Later, Hoogeboom and Dullaert [2] also studied a similar problem to the one of [1], where a specific time spread constant is defined for each customer. In [2], a rolling horizon approach is employed instead of using the periodic setting of [1]. According to their rolling horizon approach, arrival time slots from the previous periods are deleted from the solution region in each period, then, a VRP with multiple time windows is solved for each period. They proposed an iterated granular tabu search algorithm for the problem, as well as four penalization approaches to permit infeasible solutions through the heuristic procedure. Note that these two studies [1,2] mentioned above assume that all vehicles are identical, i.e., have the same capacity, and the demand of a customer is the same in all periods. However, in real-life, the CIT companies generally have a heterogeneous fleet of vehicles, and the demand of a customer varies from one period to another due to the seasonality. Consequently, in this paper, we extend the PVRP-TWTS of [1] by considering heterogeneous vehicles with different capacities and customers with period-specific demands to reflect the real practice. Lack of exact solution procedures for this novel extension of the PVRP-TWTS is a significant gap in the literature that needs to be filled. This study aims to fill this gap by proposing novel MILP and CP models for this more realistic extension of the PVRP-TWTS, as exact solution methods. Additionally, several valid inequalities are developed to strengthen the proposed MILP and CP models, as well as a simple lower bound. Even though some studies have applied CP to solve the VRP [5][6][7][8], to the best of our knowledge, the CP model is proposed for the PVRP with time spread constraints in this study, for the first time in literature. This study also aims to emphasize the impact of CP approaches and their increasing attractiveness by showing its effectiveness on such a complex PVRP-TWTS. Finally, we intend this work to be useful to different researchers and experts in the field as a beginning point for their CP and/or PVRP research attempts.

Problem definition
In this paper, we consider the PVRP-TWTS for the ATM loading process of the banks. Generally, banks plan the visits of the vehicles to the customers (ATMs) over several days, where the customers must be visited within their time windows each day. Namely, in each day (period), they assign customers to the vehicles, and, then, they plan the routes of the vehicles regarding these assignments. In the ATM loading process, a team of workers departs from the depot, loads the ATMs, which are on their path, and turns back to the depot again. Note that the workers should load each ATM within specified periods, i.e., time windows, preferably outside the working hours of the customers. Furthermore, since each vehicle has an individual capacity, the total demand of the customers on the vehicle's route should not exceed the capacity of the vehicle. As the demand of the customers (ATMs) generally varies from day to day, each customer has a different demand in each period. Since the vehicles are carrying valuable goods, i.e., money for the ATMs, arrivals earlier than the start of the customer's time window are forbidden to ensure the security of the vehicles. Furthermore, the arrival times of the vehicles to a customer should be irregular over the planning periods. Namely, the arrival times of any two visits to the same customer must vary by at least a certain amount of time over the periods, without violating the time window of the customer. In this way, the arrival times to a customer (ATM) vary between the periods. Since it is going to be hard to track the vehicles, the exact time of the visits to the ATMs would not be known. Consequently, the goal of the banks is to provide the security of transfers over several planning periods with minimum travel time. Note that this problem with additional limitations, i.e., the PVRP-TWTS, can also be seen in other high-value shipment transportation processes of the CIT companies. The PVRP-TWTS can be formally defined on a directed graph with a planning horizon P of several periods, i.e., days. The set of nodes = {0, 1, … , , + 1} contains the set of customers, where nodes 0 and + 1 denote the depot. There is a set of heterogeneous vehicles with different capacities, where each vehicle has an individual capacity . The service at each customer includes providing an amount , of goods in a period ∈ , with a service time , without arriving earlier than and not after the . Note that waiting at a customer until the beginning of its time window is forbidden. Furthermore, for any two periods, and ′, a minimum time constant should split the arrival times at a customer, corresponding to time spread constraints. The PVRP-TWTS aims to construct a set of vehicle routes in each period that minimizes total travel time while distinguishing the arrival times. In a route, a vehicle departs from the depot with a load that complies with its capacity; visits a subgroup of customers, which are assigned to it, within their time windows; delivers to these customers their demands; and returns to the depot. Note that the PVRP-TWTS considers both the assignment of the vehicles to customers and the routing of the vehicles. The PVRP-TWTS is NP-hard, since the basic single-period version of the problem, i.e., the vehicle routing problem, is known to be NP-hard [26]. A variant of the PVRP-TWTS was initially studied in the literature by [1], and a MILP formulation was provided for the problem. However, the authors [1] assume that all vehicles are identical, and the demand of a customer is the same in all periods. In this paper, we consider heterogeneous vehicles with different capacities and customers with period-specific demands to reflect the real practice. Consequently, we extend the MILP model of [1] by including heterogeneous vehicles and customers with period-specific demands. Furthermore, we propose a novel CP model for the studied PVRP-TWTS as well as valid inequalities and a lower bound.

Mathematical models
As mentioned in Section 3, this paper formulates a MILP model and a CP model for the PVRP-TWTS using the OPL modeling language [27]. All parameters and sets of the PVRP-TWTS used in the proposed MILP and CP models are presented as follows.

Sets & parameters
: Set of nodes = {0, 1,…, n, n+1} : Set of customers = {1,…, n} : Set of vehicles : Set of periods : Capacity of vehicle ∈ , : Travel time from node ∈ to node ∈ , : Demand of node ∈ at period ∈ . Note that the demand of the depot nodes equal to 0 in all periods. : Start time of the time window of node ∈ : End time of the time window of node ∈ : Service time at node ∈ . Note that the service time of the depot nodes equal to 0.

Objective function
, , + + , − ( , , + + , + (1 − , , , ) ≥ , , ∀ , ∈ | ≠ , ∈ , ∈ , , ≥ ∀ ∈ , ∈ , ∈ (10) , , ≤ ∀ ∈ , ∈ , ∈ , , − , ′, ′ + (1 − , , ′, , ′ ) ≥ ∀ ∈ , , ′ ∈ , , ′ ∈ | ≠ ′ (12) , ′, ′ − , , + * , , ′, , ′ ≥ ∀ ∈ , , ′ ∈ , , ′ ∈ | ≠ ′ (13) The objective function of the PVRP-TWTS (Eq.1) minimizes the total travel time of the vehicles. The model aims to reduce the total cost of the routes within the planning horizon. The constraint (Eq.2) ensures that, in each period, each customer is visited precisely once by a single vehicle, while the constraint (Eq.3) guarantees that the number of departures from the depot is not greater than the total number of vehicles in each period. Constraints (Eq.4-Eq.5) state that there must be arrival to each node except the depot node 0, and there must be a departure from each node except the depot node + 1, respectively. All vehicles should return to the depot node + 1, specified by the constraint (Eq.6). Since vehicles have different capacities, the customervehicle assignments should respect those vehicle capacities. The constraint (Eq.7) guarantees that the total demand assigned to a vehicle does not exceed its capacity. Constraints (Eq.8-Eq.9) ensure that arrival times are consistent on a route, waiting times are not allowed, and sub tours are avoided. The following constraints (Eq.10-Eq.11) are written to respect the time window restrictions of the customers by providing that arrival times at each node must be within its time window. There should be a minimum time interval between the visits of the same customer in different periods due to security reasons, which is provided by the constraints (Eq.12-Eq.13).

Valid inequalities for the MILP model
In order to strengthen the MILP model in Section 4.1, two valid inequalities are developed as follows: As mentioned before, waiting times are not allowed for the vehicles while traveling from a customer to another. The first valid inequality (Eq.14) is written concerning this assumption. Note that the latest time for a vehicle to arrive at a customer is the upper limit of its time window. If a vehicle visits a customer at the upper limit of its time window ( ) and serves it with a service time ; after that, travels to another customer with travel time ( , ), and arrives at the customer earlier than the lower limit of its time window ( ), then the vehicle should wait to serve the customer until its time window is opened. Since the waiting times are not allowed, the vehicle cannot visit the customer immediately after visiting the customer if the aforementioned situation occurs. The second inequality (Eq.15) is written concerning the time window restrictions of the nodes. Note that the earliest time for a vehicle to serve a customer is the lower limit of its time window. If a vehicle visits a customer at the lower limit of its time window ( ) and serves it with a service time ; after that, travels to another customer with travel time ( , ), and arrives at the customer later than the upper limit of its time window ( ), then, the vehicle cannot visit the customer immediately after visiting customer , as the upper limit of the time window of customer is exceeded.

Constraint programming (CP) model
The following parts present the decision variables, the objective function, and the constraints for the CP model. In the CP model, optional interval variables are defined to represent serving node with a service time regarding its time window by vehicle at period . Furthermore, the sequence variables are defined to represent the route of a vehicle at period . Each sequence constraint collects all the optional interval variables related to a certain vehicle.

Decision variables
, , : Optional interval variable for visiting the node by vehicle in period with a service time , which is defined in a domain of [ , + ]. ℎ , : Sequence variable of vehicle for period over { , , | ∈ }.

Objective function
In the objective (Eq.16), the CP model minimizes the total travel time of the vehicles over all periods. In the equation (Eq.16), the next operator displays the next node in the route of the vehicle , immediately after visiting node in period . The constraint (Eq.17) ensures that each customer is visited by only a single vehicle at each period, and the constraint (Eq.18) guarantees that each customer is visited exactly | | times over the planning horizon. Constraints (Eq.19-Eq.20) ensure that the vehicles start and end their routes at the depot in each period. Note that nodes 0 and + 1 represent the depot. Since a vehicle visits a subset of customers at each period, the total demand of those customers should not exceed the capacity of the vehicle, which is assured by the constraint (Eq.21). According to the constraints (Eq.22, Eq.23), the depot should be the first and the last node in the route of a vehicle for each period. A vehicle cannot visit more than one node at the same time. Thus, the constraint (Eq.24) states that visits of a vehicle in a route should not be overlapped, and there must be at least the traveling time of those visited customers between their visits. To provide the security of the vehicles, and differentiate the routes between periods, the model inserts at least a certain interval between the visits of the same customer in different periods using the constraint (Eq.25). The constraint (Eq.26) prohibits the waiting time for a vehicle by limiting its tour as the sum of the service time and travel time to the customers on the route.

Valid inequalities for the CP model
In order to strengthen the CP model presented in Section 4.3, the following two valid inequalities are added: ( ℎ , , , , , , ) ≠ ∀ , ∈ | ≠ , + + , < , ∈ , ∈ (27) ( ℎ , , , , , , , ) The constraint (Eq.27) ensures that if there is a waiting time between the subsequent visits of customers and , then customer should not be visited immediately after customer since the waiting times are not allowed between consecutive visits of the customers. Note that, the valid inequality (Eq.27) for the CP model is equivalent to the valid inequality (Eq.14) of the MILP model. The constraint (Eq.28) forces a customer to be visited before another customer due to the limits of their time windows if both of these customers are assigned to the same vehicle. Namely, suppose the upper limit of the time window of customer is less than the lower limit of the time window of customer . In that case, the customer should be visited before customer , if both customers and are assigned to the same vehicle.

Lower bound
To strengthen the aforementioned MILP and CP models, we also present a simple lower bound for the studied PVRP-TWTS as follows: As stated in the aforementioned models, in each period , there must be arrival to each customer from another node except the node + 1. Hence, taking the minimum traveling time for visiting customer provides a lower bound for the traveling time to the customer . Furthermore, in each period , there must be at least one return to the depot ( + 1) from one of the customers. Thus, taking the minimum traveling time for returning to the depot + 1 provides a lower bound for the traveling time to the depot + 1.
In the presented lower bound (LB) in (Eq.29), the sum of minimum traveling time of each customer and the depot over all the periods provides a lower bound for the total traveling time of the vehicles over all periods, i.e., the objective function of the PVRP-TWTS. In order to strengthen the proposed MILP and CP models, we also included this lower bound (LB) to the MILP and CP models by including a constraint ≥ .

Experimental results
Experimental results are divided into two subsections.
In the first part, the data generation procedure is explained in detail. Then, the computational results of the MILP and the CP models are compared in the second section.

Data generation
In order to evaluate the performance of the proposed MILP and CP models, we generate 280 instances by modifying the well-known benchmark set of Solomon [3] for the VRP with time windows. Originally, there are 56 instances in this set [3], where there are six subsets for the instances. The given data for these instances is the number of customers, the vehicle capacity, Euclidean coordinates of the nodes (customers and the depot), demand and service time of the nodes, and the time windows of the nodes. Instance sets R1 (12 instances) and R2 (11 instances) include randomly scattered nodes, sets C1 (9 instances), and C2 (8 instances) include clustered customers, while sets RC1 (8 instances) and RC2 (8 instances) include a combination of random and clustered settings. Instance sets of type 1 have a short scheduling horizon, allowing only a few customers to route. On the other hand, the instance sets of type 2 have a long scheduling horizon, permitting more customers to be served by the same vehicle. Within a set, the instances vary with respect to the width of the time windows.
In order to use these benchmark sets, we include additional parameters to these instances as follows. We assume three periods for the planning horizon and = ⌊ * 0.5⌋ as proposed by [1]. The travel time between two nodes is equal to the corresponding Euclidean distance, which is rounded down to the nearest integer. Since the demand of each customer varies from period to period in our PVRP-TWTS, we generate the demand of each customer in each period according to a discrete uniform distribution in the range [1,50]. We also assume that there are two types of vehicles with different capacities: the first type of vehicle (large vehicle) has a capacity of 100 units, whereas the other (small vehicle) has a capacity of 75 units. Originally, the benchmark set of Solomon [3] has 100 customers and 25 identical vehicles. Since exact solution approaches have been presented in this study, these instances are very large to test the performance of the mathematical models. Therefore, we cropped the original instances to have small-sized instance sets with varying customers from 5 to 15. Namely, we take the first customers in the original instances, and we generate 56 instances with customers, where = {5, 8, 10, 12, 15}. Note that there are 56 instances for each customer size, resulting in 280 instances in total. For each instance, we assume that the number of vehicles is equal to the number of customers. Note that there are ⌊ /2⌋ large vehicles in each instance, whereas the rest are small vehicles.

Computational results
As mentioned in Section 5.1, the computational tests have been organized on five different sets of instances, where the number of customers ranges from 5 to 15. The MILP and CP models are coded on the IBM ILOG CPLEX Optimization Studio and solved with CPLEX 12.8. Both MILP and CP results are obtained for all instances with a 30-and 60-minutes time limit on a Core i5, 3.20 GHz, 8 GB RAM computer. Note that only instances with five customers can be optimally solved within the 1-hour time limit. Since the studied PVRP-TWTS is NP-hard, optimal results cannot be obtained for the instances with a higher number of customers. The relative percentage deviation (RPD) is computed for each instance as follows: where is the objective function value obtained by a given mathematical model (MILP or CP), and is the best objective function value that is found among the proposed models under the same time limit. The average relative percentage deviation (ARPD) of the MILP and CP models under 30-and 60-minutes time limits are reported for each instance set in Table 1 and Table 2, respectively, as well as the number of best solutions found. As mentioned in Section 5.1, there are six subsets of instances (R1, R2, RC1, RC2, C1, C2) for each customer size. The first column of the table represents each subset of instances grouped according to the scheduling horizon, the distribution of the nodes, and the number of customers. In other words, the average RPD of the results that belong to the same cluster type (R, C, or RC) have the same scheduling horizon (type 1 or type 2), and the same number of customers (5, 8, 10, 12, or 15) are computed. For example, instance set R101-R112-5 represents the instances with randomly scattered customers (R) having a short scheduling horizon (type 1) and 5 customers. Also, this example has a total of 12 instances ranges between R101 and R112. Then, for a given model, "# of Best" and ARPD columns of the table report the number of best solutions found and the ARPD value for each subset of instances. Finally, for a given model, the last row of the table reports the total number of best solutions found and the average ARPD value overall. In both tables, the MILP model cannot reach any feasible solution within the given time limit for some of the sets with 12 and 15 customers. For these instance sets (with 12/15 customers), the number of instances in which the MILP cannot find any feasible solution is specified by a superscript in the "# of Best" column. As seen in Table 1, the CP model clearly outperforms the MILP model under a 30 minutes time limit in terms of the overall ARPD value. Specifically, the overall ARPD value of the MILP is 9.22%, whereas the overall ARPD value of the CP model is 0.53%. Additionally, the effectiveness of the CP model is obviously seen when the number of customers increases. Note that the CP has smaller ARPD and finds a higher number of best solutions than the MILP for the instances with a higher number of customers. Moreover, the MILP cannot obtain a feasible solution within the 30 minutes time limit for some of the instances with 12 and 15 customers and one instance with 10 customers, whereas the CP finds a feasible solution for all instances within the 30 minutes time limit. Furthermore, the MILP model provides 121 out of the 280 best solutions, whereas the number of best solutions found by the CP model is nearly twice as high as the MILP model, i.e., the CP finds 235 out of the 280 best solutions. Consequently, it can be said that the proposed CP model is very effective in solving the PVRP-TWTS in terms of both ARPD value and the number of best solutions found, as compared to the MILP model. The performances of the models are also analyzed under a longer time limit to see their changes on the solution quality. Table 2 summarizes the results of the MILP and CP models under 60 minutes time limit. As seen in Table 2, again, the CP model clearly outperforms the MILP model in terms of the overall ARPD value. Specifically, the overall ARPD value of the MILP is 9.18%, whereas the overall ARPD value of the CP model is 0.50%. Moreover, the MILP model can find 130 out of the 280 best solutions, whereas the CP model can find 236 out of the 280 best solutions. Similar to the results under a 30 minutes time limit, the effectiveness of the CP model is clearly seen when the number of customers increases. Additionally, the MILP cannot obtain a feasible solution within the 60 minutes time limit for some of the instances with 12 and 15 customers, whereas the CP finds a feasible solution for all instances within the 60 minutes time limit. As seen in Table 2, ARPD values and the number of best solutions found do not have a dramatic change compared to the results under a 30 minutes time limit. Namely, the proportion of the best solutions found and the ARPD values of the MILP and CP models are nearly the same with those obtained under a 30 minutes time limit. However, according to the "# of Best" column's superscript numbers, as the solution time limit increases, the MILP model starts to find more feasible solutions; specifically, it finds nine more feasible results.
In terms of CPU times, the CP model has always operated up to the time limit, so the CPU time of the CP model is equal to the given time limit. This also applies to the MILP model, except for the instances with five customers. The MILP model has yielded results in less than the given time limit in some of the instances with five customers. Namely, the average CPU time for the instances with five customers is 524.67 seconds for the MILP model when the limit is 30 minutes, and it is 984.43 seconds when the limit is 60 minutes. Table 3 presents the effect of the time limit on the number of best solutions found and the ARPD values. Instead of comparing the MILP and CP models under the same time limit, Table 3 compares the MILP and CP results under 30 minutes time limit with the best results obtained under the 60 minutes limit. As the time limit increases, the objective function decreases, and the solution quality increases, the best results are assumed as the best results obtained among MILP and CP models under 60 minutes time limit. Table 3 indicates the number of best solutions found and ARPD values of the MILP and CP models under 30 minutes time limit with respect to the best results. Note that these results and deviations are calculated with respect to the best solutions obtained under 60 minutes time limit.  Table 3 indicates that the solution quality of both models gets worse when the time limit decreases to 30 minutes. The number of best solutions found by the MILP model decreases from 130 to 101, while its ARPD increases from 9.18% to 11.91%. The difference between the results is not too much for the MILP model because the CP model mostly finds the best solutions. In spite of obtaining the MILP results under different time limits, the performance of the MILP model is worse than the CP model in both configurations, resulting in close numbers in the solutions. Therefore, it seems that increasing the solution time limit does not change the performance of the MILP model significantly.
As seen from Table 3, when the time limit decreases, the number of best solutions provided by the CP model decreases from 236 to 137, and the ARPD increases from 0.50% to 3.22%. We can infer from those results that the level of a time limit has a vital effect on the performance of the CP model. In order to compare all the results of both models under two different time limits, Figure 1 is depicted. In Figure  1, we compare the ARPD values of MILP and CP models under 30-and 60-minutes time limits with respect to the best results obtained among MILP and CP models under 60 minutes time limit. The detailed RPD values for all model and time limit configurations (MILP-30 min., MILP-60 min., CP-30 min., CP-60 min.) are also provided in the Appendix (Table A1-  As seen from Figure 1, all ARPD results are nearly the same for the instances with 5 and 8 customers. The ARPD results start to be differentiated as the customer size increases, and the differences are obvious when the customer size equals 15. The CP model generates the best solutions under 60 minutes time limit, and the worst results belong to the MILP model under a 30 minutes time limit. The ARPD results of the CP model under 30 minutes time limit are much better than the ARPD results of the MILP model under 60 minutes time limit, highlighting the superior performance of the CP model over the MILP model. As mentioned before, and it can also be seen from Figure 1, the level of a time limit has a vital effect on the performance of the CP model. To see the limitation of this effect, we increased the time limit to 120 minutes just for the CP model and obtained the results for some of the instance configurations. Note that the CP model dominates the MILP model for the instances with the larger sizes of the customers. Additionally, data sets having only five customers are not realistic when the real-life situation is considered. Therefore, we did not obtain the results under 120 minutes for the instances with five customers. Consequently, only a single representative instance from each combination is taken for the remaining data sets, resulting in 24 instances in total. For these instances, Figure 2 is generated to show the performance of the CP model under different time limits, i.e., 30, 60, and 120 minutes. In this figure, we compare the RPD values of the CP model under different time limits with respect to the best results obtained by the CP model under 120 minutes time limit. As seen in Figure 2, for the instances with 8 and 10 customers, the CP model found the same result under all the time limitations. However, for the instances with 12 and 15 customers, Figure 2 indicates that as the time limit increases, the RPD value of the CP model decreases. The difference is not obvious for the smaller instances, but apparent when the customer size equals 15. There is a substantial difference between the solutions under the 30-and 60-minutes time limit. Nevertheless, there is no significant difference between the solutions under the 60-and 120-minutes time limit. Finally, the detailed analysis of the experiments designates the superior performance of the CP model over the MILP model. As the given time limit increases, both models' solution quality improves, but more improvement is achieved for the CP model. However, there is not a noteworthy difference between the results under the 60-and 120-minutes time limit of the CP model. Consequently, it can be said that the CP model is very effective in solving the PVRP-TWTS in terms of the solution quality, as compared to the MILP model.

Conclusion
This paper introduces novel MILP and CP models for the PVRP-TWTS to provide safe routes for a fleet of heterogeneous vehicles while distributing valuable products, i.e., distributing the money between the ATMs throughout the city. In the PVRP-TWTS, to avoid unsecured situations and protect the transferred goods, waiting times are not allowed. The routes of the vehicles are differentiated over the planning periods by forcing the time spread constraints. In addition, the availability of the customers, i.e., working and off hours of ATMs, is regarded by employing time windows. Therefore, under those realistic constraints, this paper aims to construct safe routes for the vehicles by minimizing the total travel time, which reduces the fuel consumption of the vehicles. When the managerial implications of the study are concerned, the proposed solution approaches can be used to construct efficient and safe routes in the area of high-value shipment transportation, regarding many realistic constraints such as arrival time diversification, time windows, and waiting time restrictions. The proposed models can help the manager in solving such a complex problem effectively. The contributions of this study to the literature can be summarized as follows. Two different novel exact solution approaches, namely the MILP and CP models, were proposed to solve the PVRP-TWTS with heterogeneous vehicles and customers with periodspecific demands, which was a novel extension of the problem. In order to improve the solution quality, valid inequalities were added to both mathematical models, and a lower bound was developed for the studied problem. The instances were derived from the literature and modified according to the specifications of the considered problem. Then, the solutions of the models were compared using the modified benchmark set. The extensive computational results showed that the CP model outperforms the MILP model in terms of the solution quality. We hope that our paper can serve as a reference to researchers and experts studying in this field. The results demonstrate that, for a small number of customers up to 15, it is now possible to obtain highquality solutions using the proposed exact solution methods for such a complex problem, which is a significant development. It is useful to emphasize that the CP model is highly compact and solved by a blackbox, commercial solver. Although the potential future contribution that can be obtained through the CP model is evident, the possible drawback, i.e., the loss of solution quality and the high CPU time caused by the unavoidable increase of customer nodes, can be substantial. Therefore, as future research, new techniques for the CP solution approach, such as employing constructive heuristics with CP, will be experimented extensively to handle this drawback before the practical implementation. Since exact solution methods were proposed in this study, the limitation of the paper is that the proposed models can solve the problem for a small number of customers, up to 15 customers. Additionally, in the deeper analysis of the CP model, an effective lower bound is also necessary to find the proven optimal solution. Therefore, we are planning to formulate a hybrid model that employs both MILP and CP models to reduce the CPU time and provide more qualified solutions. Furthermore, heuristic, meta-heuristic, or matheuristic algorithms can be proposed for the PVRP-TWTS to solve larger instances.    This work is licensed under a Creative Commons Attribution 4.0 International License. The authors retain ownership of the copyright for their article, but they allow anyone to download, reuse, reprint, modify, distribute, and/or copy articles in IJOCTA, so long as the original authors and source are credited. To see the complete license contents, please visit http://creativecommons.org/licenses/by/4.0/.