Simulated Annealing with Restart Strategy for the Path Cover Problem with Time Windows

This research presents a variant of the vehicle routing problem known as the path cover problem with time windows (PCPTW), in which each vehicle starts with a particular customer and finishes its route at another customer. The vehicles serve each customer within the customer’s time windows. PCPTW is motivated by a practical strategy for companies to reduce operational cost by hiring freelance workers, thus allowing workers to directly service customers without reporting to the office. A mathematical programming model is formulated for the problem. This research also proposes a simulated annealing heuristic with restart strategy (SARS) to solve PCPTW and test it on several benchmark datasets. Computational results indicate that the proposed SARS effectively solves PCPTW.


Introduction
The vehicle routing problem (VRP) is widely used to address logistics and transportation problems, and on the basis of the literature over the last few decades many variants of it have been extensively studied. Examples include VRP with time windows [1][2][3][4], VRP with stochastic demand [5][6][7], time-dependent VRP [8][9][10], and VRP with split delivery [11][12][13]. Most of these studies are devoted to dealing concerns about vehicle routes in a logistics and transportation system, the resources of which belong permanently to the companies. However, in real-life situations, resources are commonly not owned permanently by companies, and this situation must be considered by decision makers.
The present work is motivated by the practical strategy currently adopted by many several small-and medium-sized enterprises (SMEs) that hire freelance workers to reduce operational costs. This practice allows workers to start driving from either their homes or certain geographical locations and return directly to them without reporting to the main office. The freelance workers can be private car owners or occasional drivers provided by delivery service companies such as Uber and Grab. Hence, hiring freelance workers provides vehicle fleet operators or companies with flexibility in planning service routes. By using this strategy, companies can assign efficient paths that disregard the distance between workers' homes and the depot. Companies will only be concerned with the path and thus this research extends PCP by considering time windows to propose the path cover problem with time windows (PCPTW). Figure 1 depicts an illustration of PCPTW, whose structure is similar to PCP, but also covers the service time windows of customers. As shown in Figure 1, a dummy depot and a dummy route disregard the distance from the depot to the customer.
The vehicle routing problems that consider the hiring of outsourced vehicles are commonly modeled as open vehicle routing problems (OVRPs), as presented by Salari et al. [15]. Outsourced vehicles provided by third-party logistics operators are hired to fulfill customer demands. When each customer is restricted to an earliest and latest visit time, the problem is extended to the 'open vehicle routing problem with time windows' (OV-RPTW), as discussed in Repoussis et al. [16]. After finishing the tasks, the outsourced vehicles do not return to the depot, since they belong to third-party logistics operators. PCPTW is different from OVRP and OVRPTW in that the drivers are allowed to visit their first customer directly from their homes instead of from the main office. In practice, each route is assigned to a driver who can visit the customers within their respective time windows. This is similar to the underlying assumption of VRPTW that each driver is available to start the route at a time that allows her/him to meet the customers' time window restrictions, including that of the first customer. VRP is an NP-hard combinatorial optimization problem, and hence PCPTW is also an NP-hard problem. To address it, this work proposes a simulated annealing with restart strategy (SARS) heuristic, which is a diversified version of simulated annealing (SA) that utilizes a restart strategy to escape from local optima. This study conducts extensive computational experiments on a PCP dataset, well-known VRPTW benchmark problems, and a PCPTW dataset to investigate the performance of the proposed method on the three problems. The vehicle routing problems that consider the hiring of outsourced vehicles are commonly modeled as open vehicle routing problems (OVRPs), as presented by Salari et al. [15]. Outsourced vehicles provided by third-party logistics operators are hired to fulfill customer demands. When each customer is restricted to an earliest and latest visit time, the problem is extended to the 'open vehicle routing problem with time windows' (OVRPTW), as discussed in Repoussis et al. [16]. After finishing the tasks, the outsourced vehicles do not return to the depot, since they belong to third-party logistics operators. PCPTW is different from OVRP and OVRPTW in that the drivers are allowed to visit their first customer directly from their homes instead of from the main office. In practice, each route is assigned to a driver who can visit the customers within their respective time windows. This is similar to the underlying assumption of VRPTW that each driver is available to start the route at a time that allows her/him to meet the customers' time window restrictions, including that of the first customer.
VRP is an NP-hard combinatorial optimization problem, and hence PCPTW is also an NP-hard problem. To address it, this work proposes a simulated annealing with restart strategy (SARS) heuristic, which is a diversified version of simulated annealing (SA) that utilizes a restart strategy to escape from local optima. This study conducts extensive computational experiments on a PCP dataset, well-known VRPTW benchmark problems, and a PCPTW dataset to investigate the performance of the proposed method on the three problems.
The four main contributions of this study are as follows. First, it considers time windows for PCP, which is closer to real-life situations. In addition, PCPTW is different from VRPTW [17], because it considers the open-ended paths of vehicles instead of the closed vehicle loops in VRPTW. Second, the proposed SARS is compared with an exact solver on the PCPTW dataset, which is adapted from Solomon's VRPTW dataset. SARS obtains optimal solutions for all the small instances and can improve the upper bounds for medium and large instances. This research also investigates the advantage of the restart strategy in improving SA. Third, this work evaluates the proposed SARS's performance against other heuristics on Solomon's VRPTW and PCP datasets, finding that the algorithm performs at par with other heuristics. This research further compares it with another single-based solution heuristic and finds that SARS performs more efficiently. Fourth, the proposed SARS obtains new best-known solutions for six instances in Solomon's VRPTW dataset.
The remainder of this paper is structured as follows. Section 2 presents the literature review. Section 3 presents a mathematical model for PCPTW. Section 4 discusses the proposed SARS for PCPTW. Section 5 discusses the computational study. Section 6 concludes the paper.

Literature Review
To the best of our knowledge, the literature has no articles devoted to PCPTW. While OVRP, OVRPTW and VRPTW are variants of VRP that are close to PCPTW, PCPTW extends PCP toward real-life aspects by considering a range of service times at the customer location. In PCP, a vehicle starts at a customer and ends at another customer without visiting the depot. This problem is introduced by Yu et al. [14] and motivated by a practical strategy of employing private vehicles to undertake the transportation function of a company or organization. In OVRP, vehicles start from the depot and terminate at a customer location, and thus there is no transportation cost associated with the distance between the last customer and the depot. This variant was first introduced by Sariklis and Powell [18] as a practical technique for companies that require delivery operations from distribution companies, called freelance vehicles, which deliver goods to customers, but are not owned by the distribution companies. In contrast to OVRP, in PCPTW there is both no distance cost between the last customer and the depot, and no distance cost between the depot and the first customer. In addition, PCPTW also ensures that vehicles serve customers in a certain time range, which is not considered in OVRP. Due to the more complex characteristics of PCPTW, solving PCPTW needs more effort compared to solving OVRP. Classical heuristic methods such as the Clarke and Wright algorithm [19], petal algorithm [20], and sweep algorithm [21] are the early solution approaches for VRP and it variants [22]. In the last few decades, metaheuristics has gained much interest from VRP researchers due to their ability to solve VRP related problems compared to the classical heuristics [23].
Many OVRP studies consider additional characteristics in real-world applications, including OVRPTW. Xia and Fu [24] proposed a variant of OVRP by considering split deliveries and solved the problem using an adaptive tabu search. The algorithm was tested on instances with up to 100 customers. Xia and Fu [25] considered time window and satisfaction rate in OVRPTW. They developed an improved tabu search algorithm for the problem. The algorithm was tested on Solomon benchmark instances with 100 customers. Lalla-Ruiz and Mes [26] presented a two-index-based mathematical formulation to solve the multi-depot OVRP. The model was tested on instances with 48 to 288 customers.
In VRPTW, vehicles start from the depot, serve all customers without violating the vehicle capacity and constraints of the customer time range, and then return to the depot. Solomon [3] solved VRPTW using sequential and parallel methods. Sequential procedures construct one route at a time until all customers are scheduled. Parallel procedures are characterized by the simultaneous construction of a number of routes. The author generated six sets of problems, which are well-known as Solomon's VRPTW dataset. This dataset serves as a reference to subsequent researchers for testing their approaches to solving VRPTW efficiently. Thereafter, several researchers have presented metaheuristic methods as strong tools for solving VRPTW. This section underlines the methods from recently published papers on VRPTW, and the average performance of these methods is close to the best-known solution (BKS). In recently published papers, several algorithms have been tested on Solomon's VRPTW dataset with 100 customers. The methods include hybrid multi-objective evolutionary algorithm [27], column generation heuristic (CGH) [28], a hybrid of tabu search and ant colony optimization [17], parallel iterated tabu search heuristic (PITSH) [29], set-based particle swarm optimization [30], cooperative population learning algorithm (CPLA) [31], hybrid shuffled frog leaping algorithm (HSFLA) [32], meta-harmony search algorithm [33], and the most recent tabu-artificial bee colony (Tabu-ABC) [34].
Yang [35] stated that metaheuristic algorithms could be classified as single-based and population-based solution heuristics. The methods mentioned in recently published papers on VRPTW, such as the PITSH, are included in the single-based solution heuristics, whereas the other methods are included in the population-based solution heuristic. Our SARS is included in the single-based solution heuristic.
In addition to the VRPTW studies that focused on improving the solution method as described in the previous paragraph, researchers have also proposed many variants of VRPTW that accommodate characteristics in real problems. For example, Keskin and Çatay [36] proposed a VRPTW variant in which freight is distributed by fast-charging electric vehicles. They solved the problem by an adaptive large neighborhood search. Goel et al. [37] considered stochastic customers' demands and stochastic service times in VRPTW. The authors developed a mathematical model and a modified ant colony system to solve the problem. Energy consumption in cold chain logistics also employed VRPTW as reported in Song et al. [38]. An improved artificial fish swarm algorithm was proposed to solve realistic instances.

Mathematical Model
PCPTW is characterized as follows. Let G = (V, A) be a complete directed graph, where V = {0, 1, . . . n} is a vertex set, and A = {(i, j) : i, j ∈ V, i = j} is an arc set. Vertex 0 represents a dummy depot, and V 0 = V\{0} denotes a customer set in which each customer (i ∈ V 0 ) has demand q i and service duration s. Each (i, j) is associated with travel cost c ij , which is also interpreted in this research as travel time. Each customer i can be visited within time window (e i , l i ), where e i is opening time, and l i is closing time; m denotes vehicle availability to serve the customers, and each vehicle has capacity C, activation cost h, and maximum tour length T; and u i and y i represent vehicle load and arrival time at customer i, respectively. The objective function of PCPTW is to determine the service paths that minimize the total transportation costs by employing freelance vehicles to serve the demand of customers under the following constraints: the customers must be serviced exactly once, the total load of a vehicle does not exceed C, the duration of any route does not exceed T, and the vehicle serves the customer within their time windows. The assumptions of this model are realized in the constraints, which are reasonable in reality. For example, the total load of the vehicle cannot exceed its capacity, the route representing the driver has a standard work time, and customers have a certain time range to receive services from other parties.
The mathematical model of PCPTW is a slightly modified version of that shown in Yu et al. [17] by disregarding distance costs between customers and the depot and adding the activation cost of the vehicles to calculate the objection function. A binary variable, x ij , is used, and x ij = 1 if a vehicle moves from customer i to customer j. The mathematical model is as follows.
The objective function (1) minimizes the total transportation cost, which consists of total travel distance and total vehicle activation cost. Constraints (2) and (3) ensure that all customers are served exactly once. Constraints (4) and (5) ensure that vehicle availability is not exceeded. Constraints (6) and (7) impose the capacity and connectivity of a route. In constraints (8) and (9), L is a large number and maintains that each vehicle ends the service no later than T. Constraint (10) denotes that the customers must be served within the time windows. Note that the first customer must be served in the relevant time window, although this model disregards the distance between the depot and its location, since only the possible freelance worker as described in Section 1 can handle the task. Binary integrality is guaranteed through constraint (11).

Solution Method
SARS targets the solution of PCPTW. This algorithm is a diversified version of SA to produce improved results. SA itself was introduced by Metropolis et al. [39] and then by Kirkpatrick et al. [40]. Eglese [41] popularized the heuristic by conducting an SA overview in the combinatorial optimization problem. According to the literature, SA has been successfully applied to many NP-hard combinatorial problems [42][43][44][45][46].
SA is a local search-based algorithm consisting of a mechanism to escape from local optima. It belongs to a single-solution-based algorithm that iteratively explores a better objective value. This algorithm is a simple local search, which starts with an initial solution. It then selects a neighborhood solution of the current solution at each iteration. If the objective value of the neighborhood solution is better than that of the current solution, then the neighborhood solution replaces the current solution; otherwise, the current solution is retained. SA allows a worse-neighborhood solution to replace the current solution with a small probability; thus, the procedure can escape being trapped at a local optimum. This feature advantage of SA motivates this research to employee the SA-based heuristic in order to solve PCPTW.
In many cases, SA-based search procedures may require more diversification mechanisms to avoid being trapped at the local optima. Without these diversification mechanisms, the search space of these SA-based search procedures may be confined to a small region of the solution space. Therefore, the probability of achieving a global optimum may be reduced [47]. A popular diversification mechanism is the restart strategy, which can be used to guide the search escapes from local optima. Thus, this work proposes SARS, which combines the advantages of an SA-based algorithm in fast search convergence and the restart strategy in escaping the local optima, herein to obtain better results compared with merely using SA. The algorithm has been successful in solving combinatorial optimization problems [48]. The remaining subsections explain the solution representation, neighborhood mechanism, parameters used, and details of the SARS procedure.

Solution Representation
The solution is represented by a string of numbers consisting of a permutation of n customers denoted by the set {1, 2, . . . , n} and N dummy zeros. Each of the N dummy zeros represents the dummy depot and is a sign of constructing a new route, although the current route does not violate any restrictions. N dummy is calculated as x represents the smallest integer that is larger than or equal to x. Furthermore, the jth non-zero string in the solution representation denotes the jth customer to be visited. Thus, the first non-zero number is the first customer of the first route. Subsequent customers are added one at a time to complete the current route, provided that the time window constraints of the customers and the maximum length of the route of the vehicle are not violated. If adding a customer to the current route violates either the time window constraint of the customer, the time window constraint of the depot, or the maximum length of route constraint, then this customer is cancelled and the current route is terminated; hereafter a new route is started whenever feasible.
The vehicle must start within the range of the service time window of the customer. If the vehicle arrives at the customer earlier than the opening time, then its service must wait until the opening time. This solution representation always maintains a feasible PCPTW solution.

Neighborhood
This study uses a standard neighborhood search mechanism that includes swap, insertion, and reverse movements, similar to Yu et al. [49] and Yu and Lin [50]. The swap movement randomly selects the ith and jth customers of X and then exchanges their positions. The insertion movement is performed by randomly selecting the ith customer of X and then inserting it into the position immediately before another randomly selected jth customer of X. The reverse movement is conducted by reversing the sequence from the ith until the jth customers of X. The probabilities of performing each movement is set to be 1/3. After a movement is executed, the new solution must be encoded again to ensure its feasibility. Yu and Lin [50] illustrated the implementation of the movements.

Parameters Used
The proposed SARS consists of six parameters: T 0 , T F , α, I iter , N non-improving , and M restart . T 0 represents the initial temperature, T F denotes the final temperature, α denotes controlling the temperature reductions, I iter denotes the number of iterations performed at a certain temperature, and N non-improving declares the number of allowable maximum consecutive non-improving temperature reductions that the value of the objective function has not improved. Finally, M restart is the maximum number of restarts.

SARS Heuristic
The SA starts with an initial solution X generated randomly. The initial temperature is set to be T 0 ; and the objective function value of X is denoted as Obj(X). The current best solution X best is then set to be X. F best denotes the value of the best objective function obtained in the current SA run and is initialized as Obj(X). F Gbest is the value of the best objective function value obtained in all SA runs and is set as zero.
At each iteration, a new solution Y is generated from N(X), the neighborhood of the current solution X, and its feasibility is checked to ensure that no constraints are violated. If Y is feasible, then its objective function value is evaluated; otherwise, it is disregarded. If Y is better than X, then it replaces X as the new current solution; otherwise, Y is accepted as the new current solution with a low probability. The probability of accepting the solution is calculated based on the Boltzmann function exp(−∆/T), with ∆ = Obj(Y) − Obj(X). This process is performed by randomly generating a number 0 < r < 1 and replacing X with Y if r < exp(−∆/T).
The current temperature is gradually reduced after performing I iter iterations at the current temperature. The cooling process occurs in this stage. The next temperature is obtained by multiplying the current temperature with cooling rate α (T = T* α). N non-improving represents the maximum allowable number of reductions in temperature if the value of the objective function has not improved.
The algorithm restarts if N non-improving is reached. When the algorithm restarts, the current temperature is reset to the initial temperature, and a new initial solution is generated randomly to initiate a new SA run. The algorithm is terminated once it reaches the maximum number of restart (M restart ). The best PCPTW solution can be derived from X Gbest after the SARS algorithm is terminated. Figure 4 presents the pseudocode of the proposed SARS heuristic.
Case p ≤ 1/3 Generate a solution Y based on X by reverse move 9.
Case 1/3 < p ≤ 2/3 Generate a solution Y based on X by insertion move 10.
Case 2/3 < p ≤ 1 Generate a solution Y based on X by swap move 11.
if Y is feasible then 12.

Computational Study
SARS was coded in C++ and performed with Intel ® Core™ i7 at 3.4 GHz CPU, 8 GB of RAM, and a 64-bit Windows 10 Pro operating system. CPLEX 12.2 was used to obtain exact solutions. The following subsections describe test instance, parameter setting, and computational results.
Three experiments are conducted to investigate the performance of the proposed SARS. The first experiment compares the proposed SARS's performance with another algorithm on PCP dataset. In the second experiment, the proposed SARS is compared with other algorithms on VRPTW benchmark problems. In the last experiment, the upper

Computational Study
SARS was coded in C++ and performed with Intel ® Core™ i7 at 3.4 GHz CPU, 8 GB of RAM, and a 64-bit Windows 10 Pro operating system. CPLEX 12.2 was used to obtain exact solutions. The following subsections describe test instance, parameter setting, and computational results.
Three experiments are conducted to investigate the performance of the proposed SARS. The first experiment compares the proposed SARS's performance with another algorithm on PCP dataset. In the second experiment, the proposed SARS is compared with other algorithms on VRPTW benchmark problems. In the last experiment, the upper bounds for each PCPTW instance are generated by CPLEX, and the SA algorithm is employed to solve the instances, in order to investigate the performance of the proposed SARS in solving PCPTW.

Test Instances
To evaluate the proposed algorithm, three datasets are used. The dataset used for the first experiment consists of 10 small and 14 large instances of PCP problems, and is adopted from Yu et al. [14]. The feature of this dataset includes total demand (TD), vehicle activation cost (h), vehicle capacity (C), service time (s), maximum route length (T), and the number of vehicle routes (m).
The dataset used for the second experiment is Solomon's VRPTW instances adopted from Sintef website (https://www.sintef.no/projectweb/top/vrptw/solomon-benchmark/, accessed on 5 July 2020). In particular, two types of datasets are used, i.e., small and large instances, each consisting of 25 and 100 customers, respectively. There is a total of 102 instances including 56 small instances and 56 large instances. These datasets vary in vehicle capacity, travel time, spatial distribution of customers, time window density, and width and are thus classified into three types of problem as follows: C-type (clustered customers), R-type (uniformly distributed customers), and RC-type (a mix of R and C types). Each of these three types of problem consists of two categories. Problem sets C1, R1, and RC1 exhibit a narrow scheduling horizon, whereas problem sets C2, R2, and RC2 show a large scheduling horizon. Vehicles with small capacities and short route times are considered setbacks to the narrow scheduling horizon. Thus, a vehicle can only serve a few customers. A large scheduling horizon includes the use of a vehicle with large capacity and long travel times; thus, more customers can be served using the same vehicle. Problem sets C1, C2, R1, R2, RC1, and RC2 consist of 9, 8, 12, 11, 8, and 8 instances, respectively.
The datasets used for the third experiment contains 168 instances including 56 small, 56 medium, and 56 large instances. Because no prior PCPTW dataset exists, the aforementioned VRPTW datasets are converted into the PCPTW datasets in the following manner. The calculation of the objective value is similar to that of VRPTW. However, the objective function value of PCPTW excludes the distances between the depot and the customers and includes the vehicle activation cost of 100. Note that double precision is utilized for calculating the distances.

Parameter Settings
The parameter values may influence the quality of the computational results. A parameter setting experiment is conducted to determine the appropriate parameter setting for the proposed algorithm. This study employs a two-level (2 k ) factorial design to set the parameters. This type of factorial design has been successfully used as a procedure to determine an effective parameter setting in heuristics [51]. For the preliminary setting, some instances are randomly selected from each of the six problem sets of Solomon's benchmark problems [3].
Six parameters are required for the 2 k factorial design: T 0 , T F , α, I iter , N non-improving , and First, the one-factor-at-a-time (OFAT) approach is used to generate low and high levels for the 2 k factorial design prior to the actual parameter setting. Low and high levels are generated for each parameter from the OFAT approach to conduct the factorial design of the experiment, as depicted in Table 2. Parameter values may influence the solution quality and computational time, and therefore this study conducts sensitivity analyses to observe the effects.  Figure 5 illustrates that I iter and T F tend to influence solution quality and computational time. As I iter or T F increases, the space for exploring solutions becomes larger, and so the chances of obtaining better solutions will increase. However, these efforts need more computational time. In contrast to I iter and T F , α does not affect solution quality or computational time, whereas T 0 does slightly affect the performance of the proposed algorithm. Increasing the value of this parameter does not always improve solution quality or reduce computational time.   Figure 5 illustrates that Iiter and TF tend to influence solution quality and computational time. As Iiter or TF increases, the space for exploring solutions becomes larger, and so the chances of obtaining better solutions will increase. However, these efforts need more computational time. In contrast to Iiter and TF, α does not affect solution quality or computational time, whereas T0 does slightly affect the performance of the proposed algorithm. Increasing the value of this parameter does not always improve solution quality or reduce computational time.   When N non-improving or M restart increases, one needs significantly more time to solve the problem, but it is not always followed by an increase in the quality of the solution. Based on the design of the experiment and the sensitivity analyses, it seems that the setting I iter = 15,000, α = 0.99, T 0 = 10, T F = 0.001, N non-improving = 100, and M restart = 7 provide the best results. Therefore, these parameter values are used to run the proposed SARS.

Verification of SARS for Solving the PCP Datasets
The first experiment's results are presented in Tables 3 and 4. The first column denotes the name of the instances. The subsequent five columns denote the feature of the instances described above. The following columns present the performance of HVTNS and the proposed SARS, respectively. Based on executing SARS in ten replications for each instance, the performance of the two algorithms is represented by both the best objective value (Obj) and the computational time (CT) of the best objective value. The last column describes the percentage gap between the best objective values of the two algorithms. Since there is no report of the computational time for small instances of PCP dataset in Yu et al. [14], the performance of HVTNS in Table 3 is only represented by the best objective value.   Table 3 shows that the proposed SARS can give a solution quality as good as HVNTS with an average gap of 0.14% in small PCP instances. Table 4 shows the computational results in large PCP instances. It implies that SARS implementation is unable to provide a better solution than HVNTS for instances p2, p3, p4, p5, and p7. However, for instances p8, p9, p10, p11, p13, and p14, the result obtained by SARS is superior compared to HVNTS in terms of objective value. Overall, the proposed SARS outperforms HVTNS in terms of average objective value of all large instance. Particularly, the average percentage gap of the objective value between SARS and HVNTS is −1.42%. Although the average computational time required by SARS is higher than that of HVTNS, it is still acceptable. The result indicates that the performance of SARS is as good as HVNTS, if not better, while the computational time is still reasonable.

Comparison of SARS with Other Algorithms on VRPTW Benchmark Instances
Tables 5-7 summarize the results of the second experiment, which show the performance of SARS in solving both small and large VRPTW benchmark instances. In each table, column "TD" represents the total travel distance, "NV" denotes the number of vehicles, and "CT" indicates the computational time (s). Table 5 presents the results of the experiment of VRPTW with 25 customers. The SARS heuristic is compared with the optimal solution obtained by previous heuristics [52]. TD is the best result of 10 runs. On the basis of the experimental results, SARS obtains the optimal solution for all instances.   [53]. b = no mention of specific computing time [54]. Gap = (TD of SARS − TD of optimal solution)/TD of optimal solution × 100%.  The italic font indicates that the proposed SARS outperforms PITSH on the corresponding instance. The bold font implies that the proposed SARS obtains a new BKS. Table 6 presents the comparison of SARS to other heuristics from recently published papers on VRPTW. BKSs are summarized from Zhang et al. [34]. The table implies that in terms of solution quality, the results of SARS are equal to other algorithms and BKS on C1 and C2 instances and still comparable on other instances. Moreover, in terms of computational time, the proposed SARS outperforms CGH, PITSH, and tabu-ABC and is comparable with HSFLA. Note that Barbucha [31] did not report computational times. From the table, the other heuristics consist of two categories, namely single-based solution algorithms for PITSH by Cordeau and Maischberger [29] and population-based solution algorithms for CGH by Alvarenga et al. [28], CPLA by Barbucha [31], HSFLA by Luo et al. [32], and Tabu-ABC by Zhang et al. [34]. The table implies that the results of SARS are equal to other heuristics in terms of BKS on C1 and C2 instances and also still comparable in the other instances. In particular, SARS performs better in terms of the average gap with BKS when compared with another single-solution based algorithm from [29]. According to the table, the average gaps are 0.75% and 3.78% for SARS and PITSH [29], respectively. Table 7 shows the comparison between SARS and PITSH [29] for each instance. The table shows that, in terms of solution quality, the proposed SARS outperforms PITSH for 34 of the 56 instances (italic font). Moreover, SARS generates six new BKSs (in bold font) in instances R204, R207, R208, RC202, RC204, and RC208. Furthermore, in terms of computational time, the proposed SARS also outperforms PITSH. Note that Cordeau and Maischberger [29] reported only average computational times of PITSH as presented in Table 6.
The third experiment compares the results of SARS and CPLEX. In addition, since SARS is modified from SA, this experiment also compares results of SARS and SA. Tables 8-10 show the comparison of the PCPTW datasets with 25, 50, and 100 customers, respectively. Each table contains the results of the three methods. Each method consists of columns NV, TD, and CT representing the number of vehicles, objective value, and computational time (s), respectively. Columns Gap1 and Gap2 show the gap between CPLEX with SA and SARS, respectively. Gap1 is calculated by ((TD of SA − TD of CPLEX)/TD of CPLEX) × 100%, whereas Gap2 is calculated by ((TD of SARS − TD of CPLEX)/TD of CPLEX) × 100%.
In Table 8, CPLEX generated an optimal solution for 37 instances and ran for 5 h to generate the upper bound for 19 instances. Both SA and SARS could also achieve all optimal solutions provided by CPLEX in significantly shorter computational time. CPLEX could not solve the remaining instances to optimality, while our SARS yielded a better solution for three instances and the same solutions for 16 other instances. The average gap of SARS is −0.01%. Table 9 shows that CPLEX achieved optimal solution for 12 instances and ran for 5 h to generate the upper bound for 44 instances. The proposed SARS achieved an optimal solution for eight instances and produced equal or better results compared with those of the upper bounds for 33 instances. The average gap is −0.48%. Table 10 shows that CPLEX also achieved an optimal solution for 8 instances, and ran for 5 h to generate upper bounds for the remaining 48 instances. The proposed SARS obtained optimal solutions for 6 instances, and the average gap is −3.33%.
The three tables report that the proposed SARS is advantageous in improving SA. The proposed SARS does not appear superior to solve PCPTW with 25 customers. As the size of the problem increases, the algorithm shows its superiority compared with SA heuristic. This phenomenon is shown in Table 8 for PCPTW with 25 customers, where averages of Gap1 and Gap2 are equal at −0.01%. The values are different in Tables 9 and 10 for PCPTW with 50 customers and 100 customers, respectively. In these two tables, average Gap2 is larger than average Gap1, implying that SARS is better than SA in terms of solution quality. In terms of computational time, SARS needs more time to implement the restart strategy. The higher computational time of SARS is still acceptable since the results of SARS are better than those obtained by SA.    The values of the table are summarized from the results presented in  Tables 8-10. The percentage gap in the last column shows that for small size instances, both methods result in the same number of instances, yielding the optimal solution. The last row of the Gap column presents the percentage gap between SARS and CPLEX in terms of solution quality. It shows that the proposed SARS has a gap of less than 1% compared to that of CPLEX, yet SARS still has faster computational time than CPLEX.

Analysis of the Restart Strategy
This subsection analyzes the restart strategy. Figure 6 shows the convergence history of SA and SARS, which appears in blue and black lines, respectively. Both methods are executed with one run on instance RC207 of the PCPTW dataset. Based on the figure, both methods are able to converge during the searching process. However, SARS is able to give a better solution compared to SA.
OR PEER REVIEW 22 of 25

Analysis of the Restart Strategy
This subsection analyzes the restart strategy. Figure 6 shows the convergence history of SA and SARS, which appears in blue and black lines, respectively. Both methods are executed with one run on instance RC207 of the PCPTW dataset. Based on the figure, both methods are able to converge during the searching process. However, SARS is able to give a better solution compared to SA.

Conclusions and Future Research
The scientific novelty of this study is the development of a mathematical programming model and an effective SARS algorithm for solving PCPTW, which can be implemented to extend PCP by considering customer time windows in small-and mediumsized enterprises that hire freelance workers to reduce operational cost. The mathematical model of PCPTW is developed and solved by using CPLEX. SARS, which is an extended SA that utilizes a restart strategy, is proposed to tackle PCPTW.
To investigate the effectiveness of the proposed algorithm, we test it on PCP and well-known VRPTW datasets. The experimental results show that SARS is comparable with HVTNS in solving PCP instances and with several algorithms on VRPTW datasets. Moreover, on VRPTW datasets, SARS outperforms another single-solution based algorithm in terms of the average gap with BKS, and generates new BKSs for six large instances.
In solving the PCPTW dataset, the experimental results show that SARS outperforms CPLEX in terms of the solution quality and the computational time for all sizes of the 1,300 1,600 1,900 2,200 2,500 0 500,000 1,000,000 1,500,000 2,000,000 2,500,000 3,000,000 Objective Iteration SA SARS Figure 6. Convergence history of SA and SARS.

Conclusions and Future Research
The scientific novelty of this study is the development of a mathematical programming model and an effective SARS algorithm for solving PCPTW, which can be implemented to extend PCP by considering customer time windows in small-and medium-sized enterprises that hire freelance workers to reduce operational cost. The mathematical model of PCPTW is developed and solved by using CPLEX. SARS, which is an extended SA that utilizes a restart strategy, is proposed to tackle PCPTW.
To investigate the effectiveness of the proposed algorithm, we test it on PCP and well-known VRPTW datasets. The experimental results show that SARS is comparable with HVTNS in solving PCP instances and with several algorithms on VRPTW datasets. Moreover, on VRPTW datasets, SARS outperforms another single-solution based algorithm in terms of the average gap with BKS, and generates new BKSs for six large instances.
In solving the PCPTW dataset, the experimental results show that SARS outperforms CPLEX in terms of the solution quality and the computational time for all sizes of the problem scale. Furthermore, utilizing the restart strategy enables SARS to outperform SA in terms of solution quality. Although employing the restart strategy needs more computational time, the results indicate that the longer computational time of SARS is still reasonable.
The developed formulation and solution techniques can be extended to a practical strategy for companies in order to reduce operational cost by hiring freelance workers. However, current studies are limited to considering cost minimization. In practice, the proposed model might not yet cover several realistic situations that could arise.
Future research can consider other realistic situations and use a more practical dataset. Factors beyond cost optimization such as customer satisfaction and workers' preferences are becoming a social topic that cannot be ignored and should be addressed. Moreover, the involvement of a heterogeneous fleet, split delivery, and multiple time windows can also extend the problem. Optimization methods such as particle swarm optimization, genetic algorithm, and other advanced algorithms can be used further to increase solution quality.