Performance Comparison between Particle Swarm Optimization and Differential Evolution Algorithms for Postman Delivery Routing Problem

Generally, transportation costs account for approximately half of the total operation expenses of a logistics firm. Therefore, any effort to optimize the planning of vehicle routing would be substantially beneficial to the company. This study focuses on a postman delivery routing problem of the Chiang Rai post office, located in the Chiang Rai province of Thailand. In this study, two metaheuristic methods—particle swarm optimization (PSO) and differential evolution (DE)—were applied with particular solution representation to find delivery routings with minimum travel distances. The performances of PSO and DE were compared along with those from current practices. The results showed that PSO and DE clearly outperformed the actual routing of the current practices in all the operational days examined. Moreover, DE performances were notably superior to those


Introduction
The topic of global warming has been heavily addressed in recent decades due to increasing greenhouse gases [1]. Unsurprisingly, the transportation industry globally accounts for a large portion (mainly CO 2 emission), amounting to approximately 20% [2]. Since our economics and life wellness deeply depend on the industry, decent mitigation efforts are needed to enhance its operation efficiency. For logistics organizations, with transportation costs being their largest expenditure, an improved operation not only reduce the environmentally harmful gas, but also reduces the operational costs such as fuel, manpower and vehicle depreciation, significantly augmenting profits in the process [3]. As a result, a large logistics firm that handles thousands of transport vehicles per operational day should absolutely adopt an effective transportation management and scheduling system. The post office is one of the logistics establishments that has to deal with such management, involving area assignments for delivery vehicles and the type and number of vehicles, as well as delivery time windows. Some even engage in detailed routing of each vehicle, trying to determine the most optimized ways to deliver parcels to customers [4]. Such a challenge is essentially called the vehicle routing problem (VRP). However, most post offices do not feel the need to find solutions to the VRP and simply rely on the experience of their delivery drivers, which might be because of the high costs and complications involved in deploying such systems. Nevertheless, relying on skills and experience alone might not be suitable for a problem with exceptionally high complexities.
The VRP is one of the optimization problems that seeks optimized routing for vehicles that traverse either to deliver or pickup goods. An approach to finding a solution to the VRP might be to determine a set of routes having a minimized total distance traveled by all the vehicles involved [5]. Other VRP variations include enforcing the total delivery time to be within a certain time window [6], limiting vehicle capacity [7], having stochastic customer demands [8] or allowing customers to select their delivery options [9]. Generally, the problem sizes (such as number of vehicles and customer nodes) of the VRP that can be solved for exact optimized solutions by mathematical methods are limited [10], due to the inherent complexity of the VRP. However, actual VRP applications usually have large problem sizes, containing approximately hundreds or thousands of vehicles and customers (e.g., a beverage firm that has to deliver millions of bottles per day to a large number of depots and small businesses). Therefore, solving for optimized solutions by mathematical procedures alone is impractical for most applications. As a result, heuristic and metaheuristic optimization algorithms are usually selected as a pragmatic method to solve the VRP, since they can obtain feasible solutions relatively fast and the quality of those solutions is decent enough for real uses [11].
This study investigated the postman operations of the Chiang Rai post office, which is a logistics establishment that does not apply any VRP optimization techniques, and its delivery drivers have to rely on their own expertise and plans in order to traverse between customer locations. The routing of two delivery vehicles was tracked, in which they dispatched parcels to customers in their responsible areas for 50 operational days. The routing information consisted of traveled paths, distance covered and the actual geographic coordinates of customer locations which the vehicles visited each day. Then, we performed the routing optimization procedures based on those customer coordinates, using the metaheuristics algorithms for all the operational days of those delivery vehicles.
Since the delivery routes depended on the skills and experience in the area of their drivers, those routes might not be well optimized. Moreover, if new drivers were to replace the previous ones, they would have to learn anew and gain enough experience to better traverse the area. Thus, the efficiency of the delivery process was not consistent. As such, this research paper was focused on providing a smart routing decision, solving for VRP solutions such that the delivery performance could improve remarkably and be more consistent. Two metaheuristic algorithms were applied-particle swarm optimization (PSO) and differential evolution (DE)-and integrated with local search techniques. The algorithms sought to discover the most optimized routing those two delivery vehicles could take so that the total travel distance could be minimized. The solutions from the two algorithms, including those from current practices, had their performances compared in terms of the total distance covered.
The rest of the article is organized as follows. Section 2 provides literature reviews related to the VRP and the algorithms used in the study. The detailed dataset, mathematical model, algorithm frameworks and solution representation are given in Section 3. Finally, the experimental results and conclusion are provided in Sections 4 and 5, respectively.

Literature Review
The VRP, a nondeterministic polynomial time hardness (NP-hardness) problem [12], refers to an optimization problem related to routing management, first devised by Dantzig and Ramser [5] in 1959. It started as a simple problem of determining the optimal routes for a fleet of vehicles dispatching from a single terminal to multiple variants, including ones that limit vehicle capacity and delivery time, have multiple terminals or allow the customer to select a delivery option. Obtaining the exact optimal solution for the VRP is possible if the customer nodes are not too large for a mathematical solver. However, real VRP applications tend to have a substantial number of customer nodes. Therefore, heuristic methods are often chosen to solve the VRP [13].
There are several studies that focused on the VRP of post office deliveries with different strategies, objectives and constraints. Meira et al. [14] conducted an experiment to solve the multi-objective VRP with only one constraint-the route distance-using a large number of geographical delivery locations (~30,000 points) in Autur Nogueira, Brazil. Dumez et al. [9] tried to solve the VRP for a post office in which each customer could select a choice of vehicle routing problem delivery options (VRPDOs) with varying places and times, such as receiving parcels at home in the evening or at work during office hours. In addition, the VRP with a time window (VRPTW) was investigated by Niroomand et al. [15], who developed an optimization model for a network of post offices with aims to minimize the total operational cost within time constraints (i.e., office working hours).
Various algorithms were employed to solve the VRP and its variants. Most are metaheuristic, including differential evolution (DE) and particle swarm optimization (PSO). Other metaheuristic algorithms such as the large neighborhood search (LNS) [16] and hybrid intelligent algorithms [17] were also adopted for the VRP. DE is a population-based optimization method that is similar to other evolutionary algorithms in terms of conventional approaches (e.g., mutation and cross-over operations). However, it utilizes a less stochastic and greedier solving strategy [18]. DE has been successfully applied to many VRP studies. Kromer et al. [19] employed DE to find solutions to the stochastic VRP with simultaneous pickup and delivery by using actual bus service data in Anbessa, Ethiopia. Their objective was to minimize both the total travel distance and number of vehicles, and the algorithm yielded superior performance compared with that of a traditional savings algorithm. Sethanan [20] utilized DE (hybridized with a genetic operator) for a multi-trip VRP with backhauls and a heterogeneous fleet (MTVRPB). They devised an optimization method specifically for a company that needed to retrieve glass bottles (soft drinks) back from clients after delivery with a single objective (i.e., minimizing the total distance traveled). Erbao et al. [18] used both DE and genetic algorithms to find an optimal solution to the VRP with simultaneous delivery and pick-up and time windows (VRP-SDPTW). They built a mathematical model using mixed integer programming for their proposed VRP variant and yielded satisfying optimized results. Xing et al. [21] applied a hybrid discrete differential evolution (HDDE) to the split delivery VRP (SDVRP). The authors focused on vehicle routing management to deliver large supplies of hospital items split between multiple unmanned transport vehicles-a suitable solution in the pandemic era-with the objective of minimizing the total distance traveled, and the proposed HDDE algorithm produced superior performance compared with the forest-based tabu search (FBTS) and genetic algorithm (GA) in many problem instances. There are also studies that applied DE and local search to solve the VRP. For instance, Song and Don [22] proposed a DE algorithm integrated with local search to solve the capacitated VRP. They employed two strategies for the local search-two-random swap and two-opt testing of the combined method with the benchmark instances-and demonstrated that such a hybrid method performed better than other existing DE algorithms. Jun and Jian [23] modified DE further by using greedy subtour crossover in the mutation operation and modified ordered crossover in the crossover operation. In addition, local search with only the two-opt strategy was used in the decoding process.
In addition to DE, PSO was also used frequently in the VRP field. Moghaddam et al. [24] developed an advanced PSO to solve a VRP variant where the demands of customers were uncertain and demand distribution was undetermined, a special form of capacitated VRP (CVRP). They devised a new decoding method for the PSO and were able to obtain superior solutions compared with those from well-known studies using the same test instances. Another variant of the capacitated VRP, called the capacitated location-routing problem (CLRP), was also studied by Peng et al. [25], in which PSO was combined with clustering and local search techniques to solve the problem. Marinakis et al. [26] proposed the use of multi-adaptive PSO (MAPSO) with three different strategies to solve the VRPTW, where each strategy had a certain role in achieving optimal solutions faster, yielding promising results compared with those of other PSO variants. Chen and Shi [27] used hybrid PSO (HPSO) combined with simulated annealing (SA) to determine solutions for the multi-compartment VRP with a time window (MCVRPTW), where each vehicle possessed multiple compartments, each one containing goods with the same characteristics. They tested the algorithm with benchmark instances, and the results showed that HPSO had better performance than normal PSO (PSO without SA). Chen et al. [28] performed two types of PSO, namely discrete PSO (DPSO) and neural-like DPSO (NDPSO), on the periodic vehicle routing problem (PVRP). The PVRP deals with the tasks of finding an ideal customer service mode and determining optimal vehicle routing in accordance with the found mode. Sedighizadeh and Mazaheripour [29] combined PSO and the artificial bee colony (ABC) algorithm to solve the VRP, which converged faster than individual algorithms (PSO and ABC) and also yielded better solutions. Lastly, Yanwei et al. [30] employed PSO with selfadaptive inertia weight to solve the open VRP with time-dependent travel time (OVRPTD), a special VRP variant where vehicles do not have to return to the depot after finishing a delivery job and travel time is time-dependent, producing better solutions compared to those of PSO and the nearest-neighboring search (NNS) algorithm.

Materials and Methods
Details on the post office operation and how to obtain the necessary dataset and prepare the inputs for the study are presented in this section. Additionally, the mathematical model, algorithmic procedures of PSO and DE and process workflow for the postman delivery problem are included here.

Dataset
The dataset used in this study was obtained from the delivery operations belonging to the Chiang Rai post office in the Chiang Rai province of Thailand. The operations mainly involved postal delivery service to residential and commercial locations. The Chiang Rai post office holds delivery duties for customers in the area of the Viang and Rob-Viang districts, which are clustered into six delivery zones. Detailed data of the delivery operation were collected for 50 days, tracking two delivery vehicles equipped with Global Positioning System (GPS)-enabled android devices with location-tracking applications installed. The average number of customer locations to be delivered to per day by both vehicles was 97.22. Those two vehicles were motorcycles, each with large parcel pouches strapped on both sides, and both had duties to deliver parcels to customers in specific areas under the Chiang Rai post office's responsibility.
The GPS devices recorded routing data, including all customer locations and where the vehicles visited, in geographic coordinate form. Then, for each day of delivery operations, the total travel distance of both vehicles was extracted. Moreover, the stored latitude and longitude coordinates, which corresponded to the post office and customer locations, were used to construct a distance matrix of n × n dimensions, where n was the number of all location coordinates. Each cell of the distance matrix represented the smallest distance between a pair of location coordinates, and this distance value was determined by a python script using the googlemaps module from the Google Maps Application Programming Interface (API) [31]. The obtained matrix might not have been the symmetric type, since the shortest route from one coordinate to another may not have been the same the other way around, taking into account that one-way traffic could be enforced in some routes. In total, 50 matrices were constructed, each one representing one day of delivery operations. These matrices were essential inputs to the optimization algorithms, acting as a lookup distance table.

Mathematical Model
As mentioned earlier, this study addresses the delivery routing problem of the Chiang Rai post office. The problem is classified as a single-depot VRP, in which all vehicles that depart from the post office must return to the post office. For each operational day, a vehicle must decide on the delivery route for a set of different customers. All customers must be visited within an operational day with no specific time window. The traffic conditions were not considered in this case, due to the fact that there was little traffic in the investigated area. Thus, the delivery time and cost were primarily based on the total travel distance, and the objective of this problem was to find optimal routes with the lowest total travel distance. In addition, the capacity of the parcel pouches was assumed to be infinite, since the postman used two large pouches that could always contain all parcels, and there has been no case where the capacity was insufficient.
In this study, the postman delivery routing problem is represented as an integer programing (IP) model. The parameters and decision variables used in formulating the model are defined as follows: Indices i, j: Post office and customer locations (i, j = 1 for the post office location; i, j = 2, 3, . . . , n for customer locations). k: Vehicle (k = 1, 2, . . . , m).
Parameters d i, j : Distance from i to j U i , U j : Constant number used to avoid a subtour. N: Number of locations.

M:
Number of vehicles.
Decision Variable x i,j,k : If vehicle k travels from i to j or 0 otherwise.
The minimization of the total travel distance is expressed as The constraints are as follows: Appl. Sci. 2021, 11, x FOR PEER REVIEW 5 of 18 As mentioned earlier, this study addresses the delivery routing problem of the Chiang Rai post office. The problem is classified as a single-depot VRP, in which all vehicles that depart from the post office must return to the post office. For each operational day, a vehicle must decide on the delivery route for a set of different customers. All customers must be visited within an operational day with no specific time window. The traffic conditions were not considered in this case, due to the fact that there was little traffic in the investigated area. Thus, the delivery time and cost were primarily based on the total travel distance, and the objective of this problem was to find optimal routes with the lowest total travel distance. In addition, the capacity of the parcel pouches was assumed to be infinite, since the postman used two large pouches that could always contain all parcels, and there has been no case where the capacity was insufficient.
In this study, the postman delivery routing problem is represented as a mixed integer programing (MIP) model. The parameters and decision variables used in formulating the model are defined as follows: Indices The minimization of the total travel distance is expressed as The constraints are as follows: The objective of the model was to minimize the total travel distance as shown in Equation (1). Equations (2)

and (3) ensured that each customer was visited once and by
The objective of the model was to minimize the total travel distance as shown in Equation (1). Equations (2) and (3) ensured that each customer was visited once and by only one vehicle. The constraint in Equation (4) guaranteed the route continuity of all vehicles. Equations (5) and (6) ensured that each vehicle could take only one route, and it needed to depart from and return to the post office. Equation (7) was used to eliminate subtours. Equation (8) states that a decision variable is a binary number.

Application of PSO and DE to the Postman Delivery Routing Problem
As mentioned earlier, the postman delivery routing problem is NP-hard. Thus, solution methods based on metaheuristic approaches are more preferable than the traditional exact methods for finding solutions in real-world practices. This section presents the application of two effective metaheuristic approaches-particle swarm optimization (PSO) and differential evolution (DE)-to solve the postman delivery routing problem. A solution representation with an encoding and decoding scheme is proposed to transform a metaheuristic solution into a practical solution. The procedures of these algorithms are described below.

Particle Swarm Optimization (PSO)
Particle swarm optimization (PSO), proposed by Kennedy and Eberhart in 1995, is a population-based random search technique that imitates the behavior of fish schooling or birds flocking. In PSO, a solution is represented as a particle, and the population of solutions is called a swarm of NP particles. Each particle consists of two main attributes: The position and the velocity. The evolutionary concept of PSO is that each particle learns from the cognitive knowledge of its experiences (personal best, or pbest) and the social knowledge of the swarm (global best, or gbest) to guide itself to the better position. A particle moves to a new position using the updated velocity, which is adjusted according to the experiences of a particle. In particular, the updated velocity of a particle depends on three main control parameters: The inertia weight, personal best coefficient (cp) and global best coefficient (cg). Once a new position is obtained, the best position of each particle and the best position of the swarm are then updated as needed. These processes are repeated until a stopping criterion is met. Figure 1 shows the evolutionary procedure of the PSO.
The objective of the model was to minimize the total travel distance as shown in Equation (1). Equations (2) and (3) ensured that each customer was visited once and by only one vehicle. The constraint in Equation (4) guaranteed the route continuity of all vehicles. Equations (5) and (6) ensured that each vehicle could take only one route, and it needed to depart from and return to the post office. Equation (7) was used to eliminate subtours. Equation (8) states that a decision variable is a binary number.

Application of PSO and DE to the Postman Delivery Routing Problem
As mentioned earlier, the postman delivery routing problem is NP-hard. Thus, solution methods based on metaheuristic approaches are more preferable than the traditional exact methods for finding solutions in real-world practices. This section presents the application of two effective metaheuristic approaches-particle swarm optimization (PSO) and differential evolution (DE)-to solve the postman delivery routing problem. A solution representation with an encoding and decoding scheme is proposed to transform a metaheuristic solution into a practical solution. The procedures of these algorithms are described below.

Particle Swarm Optimization (PSO)
Particle swarm optimization (PSO), proposed by Kennedy and Eberhart in 1995, is a population-based random search technique that imitates the behavior of fish schooling or birds flocking. In PSO, a solution is represented as a particle, and the population of solutions is called a swarm of NP particles. Each particle consists of two main attributes: the position and the velocity. The evolutionary concept of PSO is that each particle learns from the cognitive knowledge of its experiences (personal best, or pbest) and the social knowledge of the swarm (global best, or gbest) to guide itself to the better position. A particle moves to a new position using the updated velocity, which is adjusted according to the experiences of a particle. In particular, the updated velocity of a particle depends on three main control parameters: the inertia weight, personal best coefficient (cp) and global best coefficient (cg). Once a new position is obtained, the best position of each particle and the best position of the swarm are then updated as needed. These processes are repeated until a stopping criterion is met. Figure 1 shows the evolutionary procedure of the PSO.

Differential Evolution (DE)
Differential evolution (DE), proposed by Storn and Price in 1995, is one of the evolutionary algorithms (EAs) used for global optimization over a continuous search space. Similar to other population-based random search techniques, DE starts by randomly generating an initial population of NV vectors. A solution in DE is represented by the D-dimensional vector, and each value in the D-dimensional space is represented as a real number. The crucial idea of DE is its mechanism for generating trial vectors by mutation and crossover operations with the use of a couple control variables: The scale factor (F) and the crossover rate (Cr). In this study, the classic but efficient mutation operation of DE/rand/1 was applied to generate a mutant vector by adding a weighted difference (F) between two randomly selected vectors to the third randomly selected vector. DE is able to generate better diverse solutions, since the best solution in the population does not exert any influence on the other solutions in the population. In addition, the mutant vector is always a solution that is not from the original population; therefore, the crossover operation in DE is always between a solution from the population and a newly generated one [32]. Then, the selection operation or replacement of an individual vector occurs only if the trial vector outperforms its corresponding current vector. As shown in Figure 2, the evolution procedure of a DE population continues through repeated cycles of three main operations-mutation, crossover and selection-until a stopping criterion is met. Differential evolution (DE), proposed by Storn and Price in 1995, is one of the evolutionary algorithms (EAs) used for global optimization over a continuous search space. Similar to other population-based random search techniques, DE starts by randomly generating an initial population of NV vectors. A solution in DE is represented by the D-dimensional vector, and each value in the D-dimensional space is represented as a real number. The crucial idea of DE is its mechanism for generating trial vectors by mutation and crossover operations with the use of a couple control variables: the scale factor (F) and the crossover rate (Cr). In this study, the classic but efficient mutation operation of DE/rand/1 was applied to generate a mutant vector by adding a weighted difference (F) between two randomly selected vectors to the third randomly selected vector. DE is able to generate better diverse solutions, since the best solution in the population does not exert any influence on the other solutions in the population. In addition, the mutant vector is always a solution that is not from the original population; therefore, the crossover operation in DE is always between a solution from the population and a newly generated one [32]. Then, the selection operation or replacement of an individual vector occurs only if the trial vector outperforms its corresponding current vector. As shown in Figure 2, the evolution procedure of a DE population continues through repeated cycles of three main operations-mutation, crossover and selection-until a stopping criterion is met.

Solution Representation
In order to apply PSO and DE for solving a combinatorial problem like the postman delivery routing problem, solution representation or solution mapping was required to transform the real number of the D-dimensional space into practical solutions for vehicle routing. Figure 3 presents the procedure of solution representation with the encoding and decoding scheme.

Solution Representation
In order to apply PSO and DE for solving a combinatorial problem like the postman delivery routing problem, solution representation or solution mapping was required to transform the real number of the D-dimensional space into practical solutions for vehicle routing. Figure 3 presents the procedure of solution representation with the encoding and decoding scheme.
The encoding process aimed to represent a solution as a string of dimensions. In this study, the postman delivery routing problem with n − 1 customers (noting that n is the number of locations, including the post office) and m vehicles was encoded as a string of D dimensions. The number of dimensions (D) was set to be equal to the total number of customers, and the dimension index stood for the customer index. Each value in a dimension was initially filled with a uniform random number in the range [0, 1]. Next, the process of decoding was divided into two stages: (1) the assignment of customers to vehicles and (2) the sequencing or routing for each vehicle. In the assignment stage, the range (max-min) of the dimension value was split by the number of vehicles as in Equation (9): where m is the number of vehicles, Min is the minimum value in the dimensions and Max is the maximum value in the dimensions. The encoding process aimed to represent a solution as a string of dimensions. In this study, the postman delivery routing problem with n − 1 customers (noting that n is the number of locations, including the post office) and m vehicles was encoded as a string of D dimensions. The number of dimensions (D) was set to be equal to the total number of customers, and the dimension index stood for the customer index. Each value in a dimension was initially filled with a uniform random number in the range [0, 1]. Next, the process of decoding was divided into two stages: (1) the assignment of customers to vehicles and (2) the sequencing or routing for each vehicle. In the assignment stage, the range (max-min) of the dimension value was split by the number of vehicles as in Equation (9): where m is the number of vehicles, Min is the minimum value in the dimensions and Max is the maximum value in the dimensions. According to Equation (9), the range of the dimension value was equally partitioned into m subsections corresponding to m vehicles. The span of each subsection is determined in Figure 4. It is noted that a subsection index denoted a vehicle index.  According to Equation (9), the range of the dimension value was equally partitioned into m subsections corresponding to m vehicles. The span of each subsection is determined in Figure 4. It is noted that a subsection index denoted a vehicle index. The encoding process aimed to represent a solution as a string of dimensions. In this study, the postman delivery routing problem with n − 1 customers (noting that n is the number of locations, including the post office) and m vehicles was encoded as a string of D dimensions. The number of dimensions (D) was set to be equal to the total number of customers, and the dimension index stood for the customer index. Each value in a dimension was initially filled with a uniform random number in the range [0, 1]. Next, the process of decoding was divided into two stages: (1) the assignment of customers to vehicles and (2) the sequencing or routing for each vehicle. In the assignment stage, the range (max-min) of the dimension value was split by the number of vehicles as in Equation (9): where m is the number of vehicles, Min is the minimum value in the dimensions and Max is the maximum value in the dimensions. According to Equation (9), the range of the dimension value was equally partitioned into m subsections corresponding to m vehicles. The span of each subsection is determined in Figure 4. It is noted that a subsection index denoted a vehicle index.  Then, a customer was assigned to a vehicle, considering its customer dimension value with the corresponding subsection. For example, if the dimension value of customer 1 fell into Section 2, that customer was assigned to vehicle 2 correspondingly. This process continued until all customers were considered. Once the assignment of customers to vehicles was decided, a heuristic called the nearest-neighbor search (NNS) was applied to determine the delivery routing of each vehicle. Thus, for each vehicle, the next visited customer was the one with the nearest distance to the current location. The completed route of a vehicle was obtained when all customers assigned to that vehicle were considered. It is noted that all vehicles needed to depart from and return to the post office.
To illustrate the process of solution representation, let us consider an example of the postman delivery routing problem with six customers and two vehicles. A distance matrix between the post office and customers is shown in Figure 5. The sample matrix in this example is symmetric; however, a matrix constructed from the actual operation might be asymmetric, depending on the locations where the delivery process takes place. Then, a customer was assigned to a vehicle, considering its customer dimension value with the corresponding subsection. For example, if the dimension value of customer 1 fell into Section 2, that customer was assigned to vehicle 2 correspondingly. This process continued until all customers were considered. Once the assignment of customers to vehicles was decided, a heuristic called the nearest-neighbor search (NNS) was applied to determine the delivery routing of each vehicle. Thus, for each vehicle, the next visited customer was the one with the nearest distance to the current location. The completed route of a vehicle was obtained when all customers assigned to that vehicle were considered. It is noted that all vehicles needed to depart from and return to the post office.
To illustrate the process of solution representation, let us consider an example of the postman delivery routing problem with six customers and two vehicles. A distance matrix between the post office and customers is shown in Figure 5. The sample matrix in this example is symmetric; however, a matrix constructed from the actual operation might be asymmetric, depending on the locations where the delivery process takes place. In the encoding process, the number of dimensions was set to six, which was equal to the number of customers. Each value in a dimension was initially generated with a uniform random number between 0 and 1, as shown in Figure 6.  Figure 7 illustrates the decoding procedure for obtaining the routing of a vehicle. First, according to Equation (9), the range of the dimension value is divided by the number of vehicles, which is (0.87 − 0.25)/2 = 0.31. In this example, two subsections are obtained according to the number of vehicles. Then, the condition of the customer assignment is set. The span of sections 1 and 2 are (0.25, 0.56) and (0.56, 0.87), respectively. When the customer dimension value falls into the first subsection, that customer is assigned to vehicle 1; otherwise, they are assigned to vehicle 2. Thus, customers 1, 2, 4 and 6 are assigned to vehicle 1, and customers 3 and 5 are assigned to vehicle 2. After that, the routing of each vehicle is determined using NNS method. According to the distance matrix in Figure 5, vehicle 1 starts from the post office, then goes to customer 6, customer 4, customer 2 and customer 1 before returning to the post office. Vehicle 2 starts from the post office, then goes to customer 5 and customer 3 before returning to the post office. The total distance of this particular example is 557 + 143 = 720 units. In the encoding process, the number of dimensions was set to six, which was equal to the number of customers. Each value in a dimension was initially generated with a uniform random number between 0 and 1, as shown in Figure 6. Then, a customer was assigned to a vehicle, considering its customer dimension value with the corresponding subsection. For example, if the dimension value of customer 1 fell into Section 2, that customer was assigned to vehicle 2 correspondingly. This process continued until all customers were considered. Once the assignment of customers to vehicles was decided, a heuristic called the nearest-neighbor search (NNS) was applied to determine the delivery routing of each vehicle. Thus, for each vehicle, the next visited customer was the one with the nearest distance to the current location. The completed route of a vehicle was obtained when all customers assigned to that vehicle were considered. It is noted that all vehicles needed to depart from and return to the post office.
To illustrate the process of solution representation, let us consider an example of the postman delivery routing problem with six customers and two vehicles. A distance matrix between the post office and customers is shown in Figure 5. The sample matrix in this example is symmetric; however, a matrix constructed from the actual operation might be asymmetric, depending on the locations where the delivery process takes place. In the encoding process, the number of dimensions was set to six, which was equal to the number of customers. Each value in a dimension was initially generated with a uniform random number between 0 and 1, as shown in Figure 6.  Figure 7 illustrates the decoding procedure for obtaining the routing of a vehicle. First, according to Equation (9), the range of the dimension value is divided by the number of vehicles, which is (0.87 − 0.25)/2 = 0.31. In this example, two subsections are obtained according to the number of vehicles. Then, the condition of the customer assignment is set. The span of sections 1 and 2 are (0.25, 0.56) and (0.56, 0.87), respectively. When the customer dimension value falls into the first subsection, that customer is assigned to vehicle 1; otherwise, they are assigned to vehicle 2. Thus, customers 1, 2, 4 and 6 are assigned to vehicle 1, and customers 3 and 5 are assigned to vehicle 2. After that, the routing of each vehicle is determined using NNS method. According to the distance matrix in Figure 5, vehicle 1 starts from the post office, then goes to customer 6, customer 4, customer 2 and customer 1 before returning to the post office. Vehicle 2 starts from the post office, then goes to customer 5 and customer 3 before returning to the post office. The total distance of this particular example is 557 + 143 = 720 units.  Figure 7 illustrates the decoding procedure for obtaining the routing of a vehicle. First, according to Equation (9), the range of the dimension value is divided by the number of vehicles, which is (0.87 − 0.25)/2 = 0.31. In this example, two subsections are obtained according to the number of vehicles. Then, the condition of the customer assignment is set. The span of sections 1 and 2 are (0.25, 0.56) and (0.56, 0.87), respectively. When the customer dimension value falls into the first subsection, that customer is assigned to vehicle 1; otherwise, they are assigned to vehicle 2. Thus, customers 1, 2, 4 and 6 are assigned to vehicle 1, and customers 3 and 5 are assigned to vehicle 2. After that, the routing of each vehicle is determined using NNS method. According to the distance matrix in Figure 5, vehicle 1 starts from the post office, then goes to customer 6, customer 4, customer 2 and customer 1 before returning to the post office. Vehicle 2 starts from the post office, then goes to customer 5 and customer 3 before returning to the post office. The total distance of this particular example is 557 + 143 = 720 units.

Experimental Results
Information for the experiment, including the algorithm parameters for PSO and DE and the computational results, are shown in this section. The results emphasize the performance comparison among PSO, DE and current practices. In addition, box plots comparing the results of DE and PSO in some instances are provided.

Parameter Setting of Metaheuristics
Typically, the performance of metaheuristics is influenced by the searching and solution mapping mechanisms, while the setting of parameters decides how well the solutions converge. For a fair comparison between PSO and DE, the number of function evaluations for both algorithms was equally set to 2500, with population sizes of 50 and 50 iterations. Tables 1 and 2 show the parameter settings of PSO and DE, respectively. It is noted that these settings were empirically determined from preliminary experiments.

Experimental Results
Information for the experiment, including the algorithm parameters for PSO and DE and the computational results, are shown in this section. The results emphasize the performance comparison among PSO, DE and current practices. In addition, box plots comparing the results of DE and PSO in some instances are provided.

Parameter Setting of Metaheuristics
Typically, the performance of metaheuristics is influenced by the searching and solution mapping mechanisms, while the setting of parameters decides how well the solutions converge. For a fair comparison between PSO and DE, the number of function evaluations for both algorithms was equally set to 2500, with population sizes of 50 and 50 iterations. Tables 1 and 2 show the parameter settings of PSO and DE, respectively. It is noted that these settings were empirically determined from preliminary experiments.

Computational Results
In this study, the metaheuristics were implemented using the Python programming language [33] under PyCharm IDE [34] version 2020.3.1 with Windows 10 as the operating system. The programs ran on an Intel ® Core TM i5-4690 CPU 3.50 GHz platform with 32 GB of RAM. The performances of the algorithms were evaluated using the dataset derived from a real case study of postman delivery in association with the Chiang Rai post office for a period of 50 days. Table 3 shows the VRP solutions (total travel distance) obtained from the proposed metaheuristics, which include the best, average and standard deviation of the total distance from 10 replicated runs for each instance. The average time duration to run one optimization cycle for each instance is also reported.  To further analyze the effectiveness of the metaheuristics for solving the postman delivery routing problem, the percentage deviation (PD) between solutions obtained from the metaheuristics and those from the current practices were calculated using Equation (10). The results for the PD are shown in Table 4: where Sol current_practice is the solution obtained from the current practice and Sol heu is the solution obtained from a metaheuristic.  According to the results in Table 3, it is obvious that both the PSO and DE algorithms performed very well in generating solutions with less distance, compared with those obtained from the current practices for all instances, with fast computing times. Both algorithms were able to provide robust solutions with small standard deviations. However, it was observed that DE clearly outperformed PSO, since DE was able to generate the better solution quality for the best and average travel distances in all instances. In Table 4, DE shows its superiority to PSO for all instances by providing higher PD values from the current practice solutions. The percentage deviation (PD) of the distance ranged from 33.34 to 45.91 with an average of 39.93 for the PSO algorithm and from 42.92 to 57.88 with an average of 51.62 for the DE algorithm.
Moreover, the performance differences between PSO and DE are visually demonstrated in Figure 8, showing box plots of the total travel distance from 10 replicated runs of the PSO and DE algorithms, which used solutions from instances 1, 20, 40 and 50. According to the results in Figure 8, it is clearly seen that the solutions obtained from DE were significantly better than those obtained by PSO.
Lastly, the behavior of the algorithms was investigated in terms of time and space complexity (i.e., how the algorithm runtime and space requirements scaled when the inputs grew). Specifically, the plots of the algorithm runtimes and memory usages versus the number of visited locations (customers plus the post office) were constructed as shown in Figures 9 and 10. The number of locations ranged from 10 to 100, with the distance matrix randomly generated specifically for the complexity examination of those algorithms. Moreover, the performance differences between PSO and DE are visually demonstrated in Figure 8, showing box plots of the total travel distance from 10 replicated runs of the PSO and DE algorithms, which used solutions from instances 1, 20, 40 and 50. According to the results in Figure 8, it is clearly seen that the solutions obtained from DE were significantly better than those obtained by PSO.  Lastly, the behavior of the algorithms was investigated in terms of time and space complexity (i.e., how the algorithm runtime and space requirements scaled when the inputs grew). Specifically, the plots of the algorithm runtimes and memory usages versus the number of visited locations (customers plus the post office) were constructed as shown in Figures 9 and 10. The number of locations ranged from 10 to 100, with the distance matrix randomly generated specifically for the complexity examination of those algorithms.   Lastly, the behavior of the algorithms was investigated in terms of time and space complexity (i.e., how the algorithm runtime and space requirements scaled when the inputs grew). Specifically, the plots of the algorithm runtimes and memory usages versus the number of visited locations (customers plus the post office) were constructed as shown in Figures 9 and 10. The number of locations ranged from 10 to 100, with the distance matrix randomly generated specifically for the complexity examination of those algorithms.   According to Figures 9 and 10, the quadratic regression was used to fit both the time and space complexity plots for DE and PSO. The coefficients of determination (R 2 ) of the selected regression model in the time complexity plot were 99.43% for DE and 99.62% for PSO, and those in the space complexity plot were 99.96% for DE and 99.97% for PSO. Therefore, the plots clearly demonstrate that the time and space complexity of both algorithms were classified as having an order of n 2 complexity, or O n 2 , since both the runtime and space requirements had quadratic relationships with the input sizes.

Discussion and Conclusions
This research employed metaheuristic algorithms, namely the PSO and DE algorithms, to solve the single-objective VRP for optimized routing solutions. The objective was to determine a solution with the lowest travel distance. The study used actual operational data from two delivery vehicles of the Chiang Rai post office over a period of 50 days. The data included geographic coordinates of the customer locations where those vehicles delivered parcels, which were then transformed into distance matrices for use by the optimization algorithms. For each delivery day, the algorithms tried to search and assign new routes to those two vehicles such that the total travel distance was minimized. The routing performances from the PSO and DE algorithms, as well as those from actual operations, were compared. The experimental results showed that both metaheuristic algorithms were far superior to the real delivery operations of the post office in all 50 days of the experiment. However, DE clearly outperformed PSO by providing outstanding solution quality for all instances. This result was mainly due to the different evolution mechanisms of DE and PSO. A PSO swarm has a high tendency to cluster rapidly and becomes stagnant, since the movement of a swarm is highly influenced by the global best solution. In addition, the diversification of DE is better than PSO, because the best solution in the population does not exert any influence on other solutions in the population. A new solution in DE would be selected only if it had better fitness, and thus the average fitness of the population would be equal or better from iteration to iteration. It is important to note that the results obtained from this study were based on the delivery constraints of the Chiang Rai post office. However, for the operations of other post offices, other constraints such as vehicle capacity, a customer's time window and the delivery time horizon could be taken into consideration. In such cases, the problem must be remodeled, and new solutions will be obtained under those circumstances.
Another essential factor was the duration for running one optimization cycle in order to obtain the optimized solution. For both the DE and PSO algorithms, the average duration was approximately six minutes. Combining the decent results and acceptable runtime duration, the methods used in this study have the potential to be employed in real post office operations. Additionally, even though there were two vehicles in this study, the proposed algorithms can also be applied to large cities in which the post office operates several delivery vehicles. In essence, the methods proposed in this study could drastically help reduce the delivery cost and time for postman delivery operations. Additionally, the vehicle routing optimization process (with about 100 customer locations) only takes minutes rather than hours, which is practical in the actual post office work environment where all parcels have to be processed and dispatched quickly.
There are some challenges to be considered if these methods are to be used in actual delivery operations. First, the post office operators need to determine a quick approach to convert customer addresses into geographic coordinates, which should be done systematically and automatically. Second, the duration to run an optimization cycle is proportional to the number of addresses and delivery vehicles. Therefore, the optimization parameters such as the population sizes and iterations should be well adjusted according to those numbers. Thus, the operators could utilize machine learning, artificial intelligence and image processing methods to quickly transform those customer addresses into definite GPS locations. They should also give enough lead time to account for the unexpected increase in processing time before all deliveries can start. In addition, the proposed PSO and DE use deterministic data, in which the customers and distance matrices are known in advance, and the routing is determined using these data. When some data are changed during delivery, the only way to adapt to such change is to rerun the algorithms. Thus, another important challenge is to develop an adaptive algorithm to respond to changes effectively.
In future studies, other metaheuristic algorithms, such as the adaptive large neighborhood search and genetic programming algorithms, are planned to be tested for their performance with datasets obtained from larger delivery operations, including not only postman delivery problems but also problems from other industries.