A Multiobjective Large Neighborhood Search Metaheuristic for the Vehicle Routing Problem with Time Windows

: The Vehicle Routing Problem with Time Windows (VRPTW) is an NP-Hard optimization problem which has been intensively studied by researchers due to its applications in real-life cases in the distribution and logistics sector. In this problem, customers deﬁne a time slot, within which they must be served by vehicles of a standard capacity. The aim is to deﬁne cost-e ﬀ ective routes, minimizing both the number of vehicles and the total traveled distance. When we seek to minimize both attributes at the same time, the problem is considered as multiobjective. Although numerous exact, heuristic and metaheuristic algorithms have been developed to solve the various vehicle routing problems, including the VRPTW, only a few of them face these problems as multiobjective. In the present paper, a Multiobjective Large Neighborhood Search (MOLNS) algorithm is developed to solve the VRPTW. The algorithm is implemented using the Python programming language, and it is evaluated in Solomon’s 56 benchmark instances with 100 customers, as well as in Gehring and Homberger’s benchmark instances with 1000 customers. The results obtained from the algorithm are compared to the best-published, in order to validate the algorithm’s e ﬃ ciency and performance. The algorithm is proven to be e ﬃ cient both in the quality of results, as it o ﬀ ers three new optimal solutions in Solomon’s dataset and produces near optimal results in most instances, and in terms of computational time, as, even in cases with up to 1000 customers, good quality results are obtained in less than 15 min. Having the potential to e ﬀ ectively solve real life distribution problems, the present paper also discusses a practical real-life application of this algorithm.


Introduction
Vehicle Routing Problems (VRP) have been intensively studied and researched for the last 60 years. The Capacitated VRP was first introduced in the seminal work of Dantzig and Ramser [1], under the name of "Truck Dispatching Problem". They generalized the Traveling Salesman Problem (TSP), adding the parameter of multiple vehicles. Since then, many changes in the initial problem have been made, and many different variants of the VRP have been created, in an attempt to correlate VRP variants with real-life distribution problems [2,3]. The Vehicle Routing Problem with Time Windows (VRPTW) is a well-known variant of the VRP which has received considerable attention in recent years and reflects many real-world scenarios, as, in distribution operations, the time window of the delivery is a crucial parameter of the problem. In conjunction with the time windows, multiple other to another promising one in the neighborhood. Tabu Search (TS) and Simulated Annealing (SA) are among the most common metaheuristics for solving the VRPTW. In the present paper, a Large Neighborhood Search (LNS) is proposed. The LNS algorithm was proposed by Shaw [19], and, from then on, multiple researchers have studied and proposed LNS algorithms, not only for the VRPTW but for other variants of the problem too. Generally, LNS algorithms work by partially destroying and then repairing the solution through insertion operators [20], providing the opportunity to explore a larger neighborhood of the current solution.
The VRPTW in many cases is faced as a single objective as the individual objective functions are composed into a single one. This approach may significantly simplify the problem; however, it is not very reliable as a small variation in the weights may lead to different solutions [21]. Tan et al. [22] and Ombuki et al. [23] are some of the first that considered the VRPTW as multiobjective, while also applying the Pareto optimality concept. In both articles, evolutionary algorithms were developed and proposed, and then tested in Solomon's benchmark instances, offering new Pareto optimal solutions. Ever since, multiple researchers have studied and considered the multiobjective nature of the problem, as advances in technology and increasing computing capabilities have emerged. In many cases, multiobjective evolutionary algorithms (MOEA) are proposed or combined with other methods due to their search capabilities and their efficiency. Some very efficient multiobjective algorithms, according to their solutions in Solomon's database, were proposed by Chiang and Hsu [24], Garcia-Najera and Bullinaria [25] and Baños et al. [26].
The contribution of this paper is noted in encountering the problem as multiobjective, as well as in the development of a Multiobjective Large Neighborhood Search (MOLNS) algorithm which exploits the solution obtained from the construction heuristic algorithm. The proposed MOLNS algorithm is sufficient in minimizing both the number of vehicles and the total traveled distance by applying removal and insertion operators while also maintaining the concept of Pareto optimality [27]. To the best of our knowledge, this is the first time a Large Neighborhood Search algorithm is applied in the multiobjective VRPTW. The algorithm is tested in Solomon's 56 benchmark instances with 100 customers as well as in Gehring and Homberger's benchmark instances with 1000 customers, and the results of the algorithm are compared to the best-published results of the literature, indicating a very efficient algorithm. The proposed algorithm considers the problem as multiobjective and can offer more competitive solutions with a viable trade-off between the quality of the solution and the computational time.
In the remainder of the paper, the VRPTW is defined and formulated in Section 2, while, in Section 3, the proposed MOLNS algorithm aiming to solve real-life VRPTW cases is thoroughly presented. The computational results of the proposed algorithm in Solomon's and in Gehring and Homberger's benchmark instances are presented in Section 4. In the same section, the expected application of the algorithm in routing and scheduling software is discussed. Finally, Section 5 contains the conclusions of the paper and highlights the need for further research.

Problem Definition and Formulation
The VRPTW is defined by a set of identical vehicles denoted by K and by a direct network G (N, A), where N is the set of nodes and A = (i, j) : i j, i, j ∈ N is the set of arcs. Node 0 represents the central depot, while N * = {N/0} represents customers. Every arc (i, j), which is a path from node i to node j, is characterized by a distance, indicated as d ij . t ij stands for the travel time from node i to node j and has the same value with d ij , as the assumption that each distance unit corresponds to one specific time unit is made. Each customer i is characterized by the demand q i , the service time needed s i and the time window (e i , l i ), where e i is the opening time and l i the closing time. Customers must be served by exactly one vehicle that may arrive any time within the time window. The arrival time is given by t i . Additionally, if the vehicle arrives before the beginning of the time window (e i ), it must wait for w k i time, until service is possible, while no vehicle may arrive after the end of the time interval (l i ).
The decision variable x k ij is equal to 1 if vehicle k drives from node i to node j, and 0 otherwise. The objective of the VRPTW is to service all customers, minimizing the number of vehicles and the total traveled distance while simultaneously ensuring that all the constraints are satisfied. The mathematical model of the VRPTW is presented below: Ob jective 1 : Ob jective 2 : k∈K j∈N The objective function in Equation (1) states that the total number of vehicles and the total traveled distance should both be minimized based on the objective function (2). The constraint (3) states that a vehicle's capacity cannot be exceeded and (4) that each customer must be served exactly once and by one vehicle. Constraints (5)-(7) make sure that each vehicle starts from the depot {0}, visits and serves a certain number of customers and finally returns to the depot {0}. Constraints (8)- (10) indicate that, for a trip from node i to node j, no vehicle may arrive at customer j after the end of the time window (l i ). Simultaneously, the time of arrival to a customer depends on the time of arrival, the waiting time and the service time of the previous served customer.

Problem Solution
The proposed MOLNS algorithm, as a first step, obtains the initial solution by applying the Time-oriented nearest neighbor algorithm, which is a construction heuristic one. Consequently, in each iteration, the removal and insertion operators are applied for improving the solution. If the newly obtained solution is better than the current best one, then the new one replaces it and is used as input in the next iteration. To evaluate whether a solution is better than another one, the Pareto optimality concept is applied, which is thoroughly analyzed in Section 3.1. The construction heuristic algorithm is described in Section 3.2, while, in Section 3.3, the removal and insertion operators of the MOLNS algorithm, as well as the general framework, are presented.

Pareto Optimal
Deliveries to end-customers and consumers constitute a key part of supply chain and logistics operations. Multiple stakeholders are involved, such as third-party logistics (3PL) companies, transportation companies or even production companies that have their own fleet of vehicles. Depending on the case, either the number of vehicles or the total traveled distance is more important or contributes more to the cost. Therefore, it is important to obtain every possible non-dominated solution, and each stakeholder makes the decision and selects the solution that is most profitable [28]. More specifically, as shown in Figure 1a, when moving from one Pareto optimal solution to another, there is always a decrease in one objective and a simultaneous increase in the second objective. More specifically, when comparing solutions X 2 and X 3 of a minimization problem, it is observed that the solution X 2 has a lower value in Objective 2 and a higher value in Objective 1. Therefore, with no further information, we cannot decide which solution is more profitable as it depends on the case.

Pareto Optimal
Deliveries to end-customers and consumers constitute a key part of supply chain and logistics operations. Multiple stakeholders are involved, such as third-party logistics (3PL) companies, transportation companies or even production companies that have their own fleet of vehicles. Depending on the case, either the number of vehicles or the total traveled distance is more important or contributes more to the cost. Therefore, it is important to obtain every possible non-dominated solution, and each stakeholder makes the decision and selects the solution that is most profitable [28]. More specifically, as shown in Figure 1a, when moving from one Pareto optimal solution to another, there is always a decrease in one objective and a simultaneous increase in the second objective. More specifically, when comparing solutions X2 and X3 of a minimization problem, it is observed that the solution X2 has a lower value in Objective 2 and a higher value in Objective 1. Therefore, with no further information, we cannot decide which solution is more profitable as it depends on the case. In the present paper, both in the phase of construction and in the phase of the LNS algorithm, every obtained solution is checked for belonging either to the dominated or the non-dominated solutions, as shown in Figure 1a. If a solution obtained from the algorithm belongs to the nondominated solutions, then the Pareto optimal front is updated. More specifically, in each iteration of the MOLNS, the set of optimal solutions (X1, X2, X3, X4, X5), that belong in the initial Pareto optimal front are selected, as presented in Figure 1a. In case solution X3' is obtained after the implementation of the removal and insertion operators of the proposed MOLNS algorithm, which improves the value of Objective 2 compared to solution X3, and the values of both Objective 1 and Objective 2 compared to solution X2, the Pareto optimal front must be updated, as shown in Figure 1b. Therefore, every solution obtained from the proposed algorithm is compared to the non-dominated solutions in order to ensure that optimal solutions are maintained and continuously updated.

Time-Oriented Nearest Neighbor Algorithm
Several heuristic algorithms have been developed for obtaining solutions for the VRP. In the present paper, the Time-Oriented Nearest Neighbor (TONN) algorithm is selected for obtaining the initial solutions, due to its simplicity and speed. The algorithm works by building routes which start from the depot (Node 0) and, consequently, the "closest" to the last visited node unrouted customer is added. In the process of searching the "closest" customer, time window as well as capacity constraints must be respected; otherwise, no customer can be added. A new route is started if no customer is found, unless no more customers are left to add, meaning that all routes have been constructed and the process is completed. Solomon [18] proposed this method in an attempt to optimally solve the VRPTW.

Objective 1
Objective 2 In the present paper, both in the phase of construction and in the phase of the LNS algorithm, every obtained solution is checked for belonging either to the dominated or the non-dominated solutions, as shown in Figure 1a. If a solution obtained from the algorithm belongs to the non-dominated solutions, then the Pareto optimal front is updated. More specifically, in each iteration of the MOLNS, the set of optimal solutions (X 1 , X 2 , X 3 , X 4 , X 5 ), that belong in the initial Pareto optimal front are selected, as presented in Figure 1a. In case solution X 3 ' is obtained after the implementation of the removal and insertion operators of the proposed MOLNS algorithm, which improves the value of Objective 2 compared to solution X 3 , and the values of both Objective 1 and Objective 2 compared to solution X 2 , the Pareto optimal front must be updated, as shown in Figure 1b. Therefore, every solution obtained from the proposed algorithm is compared to the non-dominated solutions in order to ensure that optimal solutions are maintained and continuously updated.

Time-Oriented Nearest Neighbor Algorithm
Several heuristic algorithms have been developed for obtaining solutions for the VRP. In the present paper, the Time-Oriented Nearest Neighbor (TONN) algorithm is selected for obtaining the initial solutions, due to its simplicity and speed. The algorithm works by building routes which start from the depot (Node 0) and, consequently, the "closest" to the last visited node unrouted customer is added. In the process of searching the "closest" customer, time window as well as capacity constraints must be respected; otherwise, no customer can be added. A new route is started if no customer is found, unless no more customers are left to add, meaning that all routes have been constructed and the process is completed. Solomon [18] proposed this method in an attempt to optimally solve the VRPTW.
The "closeness" of each customer is measured through metric c ij . This metric measures the direct distance between the two customers d ij , the time difference between the completion of service at i and the beginning of service at j, T ij , and the urgency of delivery to customer j, u ij .

Multiobjective Large Neighborhood Search
As described in Section 3.2, customers are added in each route based on their "closeness" to the last added customer. However, time windows affect to a great extent the initial solution, which in many cases deviates from the optimum. Therefore, after the initial solutions are constructed, the MOLNS algorithm, which partially destroys and then repairs the solution, is applied for simultaneously minimizing both the number of vehicles and the total distance traveled. The general framework of the LNS algorithm was designed by Ropke and Pisinger [20]. The removal and insertion moves are thoroughly analyzed in Sections 3.3.1 and 3.3.2.

Removal Operators
In this section, five different removal operators are described for selecting the customers which will be removed from their routes. Most operators are adapted from Ropke and Pisinger [20] and Demir et al. [29], while a new one (distant and waiting-time removal) is presented. From all removal operators, we manage to select s customers and remove them from their routes.

Route Removal
The specific operator removes an entire route from the solution. In each iteration, the route is randomly selected and "destroyed" from the solution. In a later phase, and since all customers must be served, the customers of the "destroyed" route are reinserted according to the insertion operator.

Random Customers Removal
In this operator, a specific number of customers (s) that will be removed from the existing solution is initially determined. Consequently, an empty list of size s is created and filled with customers who are randomly selected. The random nature of this operator increases the diversity, as well as the searching space. As shown in Figure 2, three customers are randomly selected (s = 3) and removed from the route that each one belongs.

Distant Customers Removal
In this operator, as with the case of random customers removal, the number of customers (s) that will be removed from the solution is determined. Consequently, for each customer k, we calculate the metric = + , where i is the previous delivered customer and j is the next delivered customer, according to k, while d indicates the distance between two customers. We then sort all customers according to the metric ck and select the s "worst" customers (highest ).

Waiting-Time Removal
For each customer to be served, the metric calculates the difference between the start of the time window and the time of arrival. More specifically, the metric for customer i is calculated by = max {0, − }, where is the time of arrival at customer i. Consequently, customers are sorted according to metric and s number of customers that have the highest waiting times ( ) are selected.

Distant and Waiting-Time Removal
This operator combines in a single metric both the distant customers removal and the waitingtime removal operators, by assigning weights to both metrics. Similar to the previous operators, a specific number of customers (s) with the worst metrics is selected.
Regardless of the removal operator, the initial solution is presented as an I list containing all customers. The sequence of the elements of I determines the order of customer visitations. When a constraint is violated (capacity or time window), then a new route begins, as shown in Figure 2a. In addition, after the removal operator is applied, the customers that no longer belong to the route are saved in another list L, as shown in Figure 2b

Insertion Operators
After a removal operator has been applied, an insertion operator must be implemented so that all customers are reinserted in the solution. In the VRPTW, there is a necessary precondition to serve all customers' needs. In the present paper, three insertion operators are applied and analyzed for improving the initial constructed solution.

Distant Customers Removal
In this operator, as with the case of random customers removal, the number of customers (s) that will be removed from the solution is determined. Consequently, for each customer k, we calculate the metric c k = d ik + d k j , where i is the previous delivered customer and j is the next delivered customer, according to k, while d indicates the distance between two customers. We then sort all customers according to the metric ck and select the s "worst" customers (highest c k ).

Waiting-Time Removal
For each customer to be served, the metric wt calculates the difference between the start of the time window and the time of arrival. More specifically, the metric wt for customer i is calculated by wt i = max{0, e i − t i }, where t i is the time of arrival at customer i. Consequently, customers are sorted according to metric wt k and s number of customers that have the highest waiting times (wt i ) are selected.

Distant and Waiting-Time Removal
This operator combines in a single metric both the distant customers removal and the waiting-time removal operators, by assigning weights to both metrics. Similar to the previous operators, a specific number of customers (s) with the worst metrics is selected.
Regardless of the removal operator, the initial solution is presented as an I list containing all customers. The sequence of the elements of I determines the order of customer visitations. When a constraint is violated (capacity or time window), then a new route begins, as shown in Figure 2a. In addition, after the removal operator is applied, the customers that no longer belong to the route are saved in another list L, as shown in Figure 2b.

Insertion Operators
After a removal operator has been applied, an insertion operator must be implemented so that all customers are reinserted in the solution. In the VRPTW, there is a necessary precondition to serve all customers' needs. In the present paper, three insertion operators are applied and analyzed for improving the initial constructed solution.

Closer Customer Insertion
This operator cooperates with all removal operators. The selected customers of list L are reinserted in the routes and more specifically in the sequence of deliveries where the total distance traveled is minimized. The specific insertion procedure is separated into two different operators. The first one searches for the minimum distance while the number of vehicles must remain the same as before the removal operator was applied. The second operator searches for the minimum distance regardless of the number of vehicles.

Insertion for Minimizing the Number of Vehicles
The specific insertion operator is combined only with the route removal operator in order to decrease the number of vehicles. If all customers of the selected route can be inserted in the rest of the routes, then the number of vehicles is reduced by one, as shown in Figure 3. Simultaneously, we seek to find the best sequence and schedule of deliveries in order to manage the distance traveled as well. However, the fact that VRPTW is a multi-objective optimization problem does not guarantee us that both the number of vehicles and the total distance traveled are simultaneously minimized. Instead, decreasing the number of vehicles by reinserting "destroyed" customers in the best-fitted sequence of deliveries may cause increased total distance traveled. Figure 3 describes how the route removal operator is combined with the specific insertion operator. More specifically, among Routes , the second one (Route 2) is randomly selected and "destroyed". In the search for inserting customers n and m in the rest of the routes, we form Algorithms 2020, 13, x FOR PEER REVIEW 8 of 17

Closer Customer Insertion
This operator cooperates with all removal operators. The selected customers of list L are reinserted in the routes and more specifically in the sequence of deliveries where the total distance traveled is minimized. The specific insertion procedure is separated into two different operators. The first one searches for the minimum distance while the number of vehicles must remain the same as before the removal operator was applied. The second operator searches for the minimum distance regardless of the number of vehicles.

Insertion for Minimizing the Number of Vehicles
The specific insertion operator is combined only with the route removal operator in order to decrease the number of vehicles. If all customers of the selected route can be inserted in the rest of the routes, then the number of vehicles is reduced by one, as shown in Figure 3. Simultaneously, we seek to find the best sequence and schedule of deliveries in order to manage the distance traveled as well. However, the fact that VRPTW is a multi-objective optimization problem does not guarantee us that both the number of vehicles and the total distance traveled are simultaneously minimized. Instead, decreasing the number of vehicles by reinserting "destroyed" customers in the best-fitted sequence of deliveries may cause increased total distance traveled. Figure 3 describes how the route removal operator is combined with the specific insertion operator. More specifically, among Routes , the second one (Route 2) is randomly selected and "destroyed". In the search for inserting customers n and m in the rest of the routes, we form Routes When implementing the removal and insertion operators, it is essential to ensure the diversity of searching and the global search capabilities. Therefore, after saving the removed customers in list L, as described in Section 3.3.1, the sequence of insertion must be considered. Repeating the same sequence of customers insertion may lead to getting stuck in a local optimum. In the presented insertion operators, a random order for inserting customers is produced, for increasing the search space of the algorithm.

General Framework
Initially, the Pareto optimal solutions constructed from the TONN algorithm are utilized as an input to the MOLNS algorithm, as described in Figure 4. In each iteration that the MOLNS algorithm is implemented, the "removal" and "insertion" operators are applied in each one of the Pareto optimal solutions. Consequently, for each solution, a new one (Xnew) is obtained. If Xnew solution constitutes a Pareto optimal solution, as it belongs to the non-dominated solutions, then the Pareto optimal front is renewed; otherwise, no changes are made. When implementing the removal and insertion operators, it is essential to ensure the diversity of searching and the global search capabilities. Therefore, after saving the removed customers in list L, as described in Section 3.3.1, the sequence of insertion must be considered. Repeating the same sequence of customers insertion may lead to getting stuck in a local optimum. In the presented insertion operators, a random order for inserting customers is produced, for increasing the search space of the algorithm.

General Framework
Initially, the Pareto optimal solutions constructed from the TONN algorithm are utilized as an input to the MOLNS algorithm, as described in Figure 4. In each iteration that the MOLNS algorithm is implemented, the "removal" and "insertion" operators are applied in each one of the Pareto optimal solutions. Consequently, for each solution, a new one (X new ) is obtained. If X new solution constitutes a Pareto optimal solution, as it belongs to the non-dominated solutions, then the Pareto optimal front is renewed; otherwise, no changes are made.

Comparisons with Best Known Results
This section describes computational experiments carried out to evaluate the performance of the proposed algorithm. The algorithm was coded in Python and ran on a personal computer with the following specifications: Intel Core i7 (Intel Corporation, California, USA)-8550U 1.8 GHz with 16 GB memory. The proposed algorithm was tested in Solomon's VRPTW benchmark instances with 100 customers (see [18]), as well as in Gehring and Homberger's instances with 1000 customers (see [30]). In the specific datasets, C1 and C2 problems have geographically clustered customers; in R1 and R2, geographical data are randomly generated; and in RC1 and RC2, a mix of random and clustered structures is included. R1, C1 and RC1 have a short scheduling horizon, narrow time windows and low vehicle capacities, respectively, while R2, C2 and RC2 have a long scheduling horizon, wide time windows and high vehicle capacities, defining the number of customers serviced by the same vehicle.
Tables 1 and 2 present a summary of the results obtained from our proposed MOLNS algorithm in Solomon's instances, as well as the best-published solutions in the literature. The column labeled BP gives the best-published results, while column MOLNS the results of our Large Neighborhood Search algorithm. Each solution is characterized by the number of vehicles (NV) needed and the total distance traveled (TD). The published work in which the optimum solution is proposed is presented in column "Reference", the computational time needed to obtain each solution is given in column "Time" in seconds and the percentage deviation in the total distance traveled between the BP and the LNS is presented in column "Deviation".
Each instance was run three times, while the computational time limit was set to be 5 min, which is a reasonable time for a logistics company to obtain the routing of their vehicles and the scheduling of their deliveries. In most cases, the best-obtained results needed less time than the time limit. Since we face the multiobjective VRPTW in some instances, more than one solution was obtained. In order for the BP solution and the solution obtained from the MOLNS to be comparable, the number of vehicles must be the same. Therefore, the deviation is measured in the total distance traveled between the BP solution and the solution of the MOLNS. It can be easily observed that instances with geographically clustered customers (C1 and C2) have more efficient results than problems with random geographical data, as the mean deviation in that case is 0.30%. In randomly generated geographical data (R), the mean deviation is higher and reaches 3.37%, while in mixed geographical customers' data (RC) it is 3.47%. Finally, the mean deviation for all instances is less than 2.85%, which can be considered very efficient, since the optimal results have emerged from multiple researchers over time. In total, 32.47% of the comparable solutions have a deviation of less than 1%, most of them

1.
Construct solutions with the Time-oriented Nearest Neighbor heuristic 2.
Save Pareto optimal solutions in P 3.
While stopping criteria not met do 4.
For each Pareto optimal solution -X p do 5.
Apply Removal and Insertion operators in X p 6.
New solution is produced -X new 7.
If X new is a Pareto optimal solution then 8.
Add X new in Pareto optimal solutions (P) 9.
Update the Pareto optimal front 10.
Pareto optimal solutions (P) do not change 12.
Next 14. End while

Comparisons with Best Known Results
This section describes computational experiments carried out to evaluate the performance of the proposed algorithm. The algorithm was coded in Python and ran on a personal computer with the following specifications: Intel Core i7 (Intel Corporation, California, USA)-8550U 1.8 GHz with 16 GB memory. The proposed algorithm was tested in Solomon's VRPTW benchmark instances with 100 customers (see [18]), as well as in Gehring and Homberger's instances with 1000 customers (see [30]). In the specific datasets, C1 and C2 problems have geographically clustered customers; in R1 and R2, geographical data are randomly generated; and in RC1 and RC2, a mix of random and clustered structures is included. R1, C1 and RC1 have a short scheduling horizon, narrow time windows and low vehicle capacities, respectively, while R2, C2 and RC2 have a long scheduling horizon, wide time windows and high vehicle capacities, defining the number of customers serviced by the same vehicle.
Tables 1 and 2 present a summary of the results obtained from our proposed MOLNS algorithm in Solomon's instances, as well as the best-published solutions in the literature. The column labeled BP gives the best-published results, while column MOLNS the results of our Large Neighborhood Search algorithm. Each solution is characterized by the number of vehicles (NV) needed and the total distance traveled (TD). The published work in which the optimum solution is proposed is presented in column "Reference", the computational time needed to obtain each solution is given in column "Time" in seconds and the percentage deviation in the total distance traveled between the BP and the LNS is presented in column "Deviation".
Each instance was run three times, while the computational time limit was set to be 5 min, which is a reasonable time for a logistics company to obtain the routing of their vehicles and the scheduling of their deliveries. In most cases, the best-obtained results needed less time than the time limit. Since we face the multiobjective VRPTW in some instances, more than one solution was obtained. In order for the BP solution and the solution obtained from the MOLNS to be comparable, the number of vehicles must be the same. Therefore, the deviation is measured in the total distance traveled between the BP solution and the solution of the MOLNS. It can be easily observed that instances with geographically clustered customers (C1 and C2) have more efficient results than problems with random geographical data, as the mean deviation in that case is 0.30%. In randomly generated geographical data (R), the mean deviation is higher and reaches 3.37%, while in mixed geographical customers' data (RC) it is 3.47%. Finally, the mean deviation for all instances is less than 2.85%, which can be considered very efficient, since the optimal results have emerged from multiple researchers over time. In total, 32.47% of the comparable solutions have a deviation of less than 1%, most of them in C instances. Simultaneously, 98.70% of the results deviate less than 10% from the optimum. Finally, three new Pareto-optimal solutions are proposed in Solomon's benchmark instances, in RC202, RC207 and R201. In RC202 problem, a decrease in the total distance traveled is accomplished, but at the expense of an increased number of vehicles. In R201 and RC202 instances, we managed to improve the solutions by decreasing the total distance traveled, while the number of vehicles remains steady. The exact routes and the schedule of deliveries in each of the three new optimal solutions that are proposed in the current research are given in Appendix A.
Additionally, when comparing results obtained from algorithms, it is important to take into consideration the computational time needed. Chiang and Hsu [24] managed to produce high-quality solutions within 0.5 min and thoroughly presented data related to the computational time for other algorithms addressing this problem. Ghoseiri and Ghannadpour [31] and Ombuki et al. [23], who also considered the VRPTW as multiobjective, did not present the computational time needed, while the proposed algorithm manages to produce high efficiency results in less than a minute (52 s on average). Of course, the CPU as well as the programming language are also important factors affecting the computational time needed. Finally, in Gehring and Homberger's benchmark instances, the number of customers vary from 200 to 1000. In the current research, the algorithm was tested in instances with 1000 customers in order to identify its efficiency in terms of quality of results and computational time, when the number of customers is extremely high. These instances, as in the case of Solomon's, are categorized into six cases: C1, C2, R1, R2, RC1 and RC2. In Table 3, the average number of vehicle as well as the average traveled distance, of both the best-published (BP) solutions and the solutions produced by the MOLNS algorithm, are presented for each case. In addition, the percentage deviation between the solutions was calculated for each objective (deviation in the NV and deviation in the TD). More specifically, the mean deviation for all cases concerning the number of vehicles is 6.63%, and 10.60% for the traveled distance. On the contrary, the corresponding percentage in the TD in Solomon's instances is 2.85%, indicating that the quality of results is affected to a small extent when the number of customers increases. However, the mean computational time needed (Time column) for each case to be solved remains reasonable and low, given the volume of customers.

Discussion on Real-Life Implementation of the Algorithm
Every algorithm presents specific advantages and disadvantages. However, regardless of the algorithm, when multiple VRP variants are involved in the optimization process, the computation time is increased. Heuristic algorithms used to be applied and implemented in routing and scheduling software, due to their ability to perform in a relatively limited search space and simultaneously produce efficient solutions in limited computing time [51]. On the other hand, metaheuristic algorithms operate a wider search space, but this does not necessarily mean that the computational time will significantly increase. The advances of technology and computer science give researchers the ability to produce near-optimal solutions, and in most cases better than those produced by heuristics. Therefore, metaheuristics should be considered for implementation in routing software solutions due to offering reliable and on-time deliveries to transportation companies.
Gayialis et al. [52], Rincon-Garcia et al. [53] and Veres et al. [54] described the design and requirements of various routing and scheduling software that need fast and efficient solutions both for the initial planning and for the rescheduling of deliveries. It is clear that practitioners and companies are looking for efficient algorithms to integrate into their software solutions for logistics operations. According to Tables 1-3, the proposed MOLNS algorithm offers near-optimal solutions in a reasonable time, while also having the ability to solve cases with a higher number of customers and rendering the algorithm ideal for the routing of vehicles and scheduling of deliveries in real-life cases. Therefore, the proposed algorithm can be implemented in a software solution, to solve the problem of routing and scheduling of deliveries effectively. That was the scope of the developed algorithm from the beginning and the current research proves that it is in line with this scope.
As noted in the Introduction, information available to the planner commonly changes in real-life cases after the completion of the initial routing and scheduling. Consequently, a rapid rescheduling procedure that is fed with dynamic data becomes a necessity to create a viable real-life routing and rescheduling system. This becomes even more clear when considering that, in real-life cases, changes may even occur while vehicles are on the road to deliver the orders, and, therefore, rapid rescheduling becomes mandatory. The algorithm presented in this paper can also support the dynamic aspect of the VRPTW, as the removal and insertion operators of the MOLNS algorithm can handle changes of the input data, such as travel times between nodes. Hence, the algorithm can be applied for both the initial plan and the rescheduling procedure by applying some minor modifications.

Conclusions and Further Research
The results obtained from the MOLNS algorithm were compared to both the best-published solutions in Solomon's dataset, as well as to those of other multiobjective algorithms, providing us significant insights about this method. Based on the conducted research, as it becomes clear, no algorithm can obtain all Pareto optimal solutions in each instance in Solomon's dataset. Instead, the results have been conjunctly obtained from multiple algorithms since the VRPTW first appeared. The results of the proposed MOLNS algorithm are near-optimal, with only small deviations from the optimal, and even include three new non-dominated solutions. The computational experiments indicate that the results are efficient and are obtained in a very realistic time. As both efficiency and speed are critical factors when implementing an algorithm in software solutions, the proposed MOLNS algorithm can be examined for implementation in such software, and this will in fact be the next step of our research.
In real-life scenarios, performing an initial plan of routes and deliveries does not necessarily fulfill the purpose of distribution services. Unforeseen events and traffic congestion, in many cases and especially in the urban environment, cause delays and generate the need for rescheduling, in order to ensure that all deliveries are executed and served on-time. The presented algorithm can be easily adapted and applied in dynamic routing cases to calculate new schedules of deliveries efficiently and quickly. The implementation of the proposed algorithm in a software solution for logistics operations, therefore, needs to be further researched. Additionally, the potential reinforcement of the algorithm to solve cases that require the dynamic rescheduling of deliveries should also be studied. Recent advances in technology, such as traffic management sensors, real-time traffic data and the Internet of things technology, can provide data to enhance such rescheduling procedures and deal with the dynamic nature of the real-life problems. In addition, the algorithm could be further enhanced by applying multicriteria decision making approaches such as the one presented by Kechagias et al. [55] to optimize its function in accordance with different use-case scenarios. Finally, provided the aforementioned technologies are successfully combined with the proposed algorithm, in the next step of our research, an advanced software solution that can support the distribution processes of logistics companies, in the urban environment, while, at the same time, focusing on the fulfillment of time-window restrictions, will be created.

Funding:
The present work was co-funded by the European Union and Greek national funds through the Operational Program "Competitiveness, Entrepreneurship and Innovation" (EPAnEK), under the call "RESEARCH-CREATEINNOVATE" (Project Code: T1EDK-00527 and Acronym: SMARTRANS).

Acknowledgments:
In this section you can acknowledge any support given which is not covered by the author contribution or funding sections. This may include administrative and technical support, or donations in kind (e.g., materials used for experiments).

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