Optimal multi-depot location decision using particle swarm optimization

Abstract

The depot locations have a significant effect on the transportation cost in the multi-depot vehicle routing problem. A two-tier particle swarm optimization framework is proposed, in which an external particle swarm optimization and an internal particle swarm optimization are used to determine the optimal depot locations and the optimal multi-depot vehicle routing problem solution, respectively. In the internal particle swarm optimization, a novel particle encoding scheme is used to minimize the computational cost by concurrently allocating the customers to depots, assigning the customers to vehicles, and determining the optimal routing path for each vehicle. The quality of the solutions is enhanced through a designed mutation local search with savings scheme. To verify the effectiveness of the proposed scheme, six standard multi-depot vehicle routing problem instances are tested and compared. It is shown that the use of the external particle swarm optimization scheme to optimize the multi-depot locations reduces the average routing distance obtained by the internal particle swarm optimization by around 13.16% on average. Furthermore, for a real-world case, the proposed two-tier particle swarm optimization scheme reduces the total routing cost by around 18%. Restated, the proposed particle swarm optimization algorithm provides an effective and efficient tool for solving practical multi-depot vehicle routing problems. Notably, the proposed scheme can be used as a reference model for obtaining the optimal locations in a variety of scheduling problems.

Introduction

Scheduling problems are one of the most common forms of problem faced by managers in the business, production, and management fields. Among such problems, the vehicle routing problem (VRP) encountered in industry is one of the most important and complex problems. For example, in Industry 4.0 (also known as smart manufacturing), it is necessary not only to optimize the production process in terms of the machines used, the manufacturing schedule applied, the resources consumed, and so on but also to integrate the production process with an optimized logistics process which ensures the timely and cost-efficient distribution of the finished products to the customer. VRPs are found not only in industry but also throughout daily life, including the delivery of letters, supermarket food, pharmaceuticals, online-purchased products, and so on. VRP applications typically vary in terms of their constraints, for example, the number of service depots, the capacity and running time of the vehicles, the time required for receiving and shipping, and the delivery time. Consequently, several variants of the VRP have been proposed, including the capacity vehicle routing problem (CVRP),^1–3 the periodic vehicle routing problem (PVRP), ⁴ the multi-depot vehicle routing problem (MDVRP), ⁵ and the vehicle routing problem with time windows (VRPTW). ⁶

Previous studies have shown the VRP to be a non-deterministic polynomial-time (NP)-hard problem. ⁷ For small-scale VRPs, optimal solutions can be obtained using an exact algorithm. For example, Baldacci et al. ³ solved the CVRP using an exact algorithm based on a two-commodity flow formulation. Contardo and Martinelli ⁵ proposed an exact algorithm for the MDVRP under capacity and route length constraints, in which the problem was formulated using vehicle-flow and set-partitioning formulations. Baldacci et al. ⁶ presented a comprehensive review of the most effective exact methods proposed in the literature for the CVRP and VRPTW over the preceding ten years. However, solving large-scale VRPs using an exact approach is cost and time prohibitive. Accordingly, various metaheuristic methods have been proposed for obtaining near-optimal solutions of large-scale VRPs within an acceptable time, including genetic algorithms (GAs), ⁸ ant colony optimization (ACO), ⁹ tabu search (TS), ¹⁰ iterated local search (ILS),^11,12 and particle swarm optimization (PSO).^13,14

Among these methods, PSO involves only a few parameters and has an easy implementation, a low complexity, and a rapid convergence performance. Furthermore, PSO achieves a good balance between an exploitation-type search behavior and an exploration-type search behavior through the use of a non-linear adjustment mechanism for the inertia weight parameter. Moreover, PSO involves both cognitive and social learning processes during evolution. As a result, it has been applied to many optimization problems, including the traveling salesman problem (TSP), ¹⁵ the multi-mode resource-constrained project scheduling problem (MRCPSP),^16,17 and the permutation flowshop scheduling problem (PFSP).^18,19Table 1 summarizes some of the more notable metaheuristics proposed in the literature for solving the MDVRP.

Table 1.

Metaheuristics proposed in literature for solving MDVRP.

Principal method employed	Author(s)	Description
Particle swarm optimization (PSO)	Geetha et al. 20	The authors used an improved k-means algorithm to reduce the MDVRP to multiple VRPs, which were then solved individually using a nested PSO scheme.
Genetic algorithm (GA)	Ho et al. 21	Two hybrid GAs were proposed for solving the MDVRP, in which the initial solutions were generated either randomly (HGA1) or using a savings method and nearest neighbor heuristic (HGA2).
Tabu search (TS) algorithm	He et al. 22	A variable cluster grouping method was used to transform the MDVRP to multiple single depot VRPs, which were then solved individually using TS.
Ant colony optimization (ACO) algorithm	Yu et al. 23	The authors presented an improved ACO method based on an enhanced pheromone updating rule and an improved mutation operation.
Hybrid heuristics algorithm	Peng et al. 24	A hybrid heuristics scheme involving GA, ACO, and simulated annealing was proposed for solving large-scale MDVRPs.

MDVRP: multi-depot vehicle routing problem; VRP: vehicle routing problem.

In previous MDVRP studies, the multi-depot locations are known in advance, and the aim is thus only to find the routing paths which minimize the total delivery cost without exceeding the vehicle capacity constraints. However, in practice, the location of the depots has a significant effect on both the choice of routing path and the total transportation cost. According to previous studies, ²⁵ the transportation cost should account for no more than 12%–15% of the total operation cost of an enterprise. Hence, in solving the MDVRP, it is desirable to consider not only the optimal vehicle routing paths but also the optimal locations of the multi-depots. Accordingly, this study proposes a two-tier PSO framework which first searches for the optimal depot locations using an external PSO and then uses an internal PSO to solve the optimal vehicle routing sequence given a knowledge of these multi-depot locations. In general, solving the MDVRP involves answering three questions, namely, (1) “which depot should offer service to which customers?” (2) “which vehicle should deliver goods to which customers?” and (3) “in what order should each vehicle visit its assigned customers?” Therefore, algorithms for solving the MDVRP must be not only effective but also computationally efficient. Exact algorithms are effective but not efficient. By contrast, approximation methods are efficient but not effective. Therefore, metaheuristic algorithms such as those presented in Table 1 are generally preferred, in which either the same scheme or different schemes are used to solve the three MDVRP sub-problems. In this study, the computational complexity of the internal PSO used to solve the MDVRP is reduced by means of a novel extended particle encoding scheme, which facilitates customer-to-depot allocation, customer-to-vehicle assignment, and optimal vehicle routing determination in a single process. The proposed encoding scheme basically treats MDVRP as a clustering problem with one vehicle one cluster. In addition, a clustering and sweeping heuristic is employed to improve the search performance by identifying a set of diverse initial solutions from which to commence the search of the solution space. As with all metaheuristic algorithms, PSO does not guarantee that the optimal solution will be found. Previous studies^17,21 have addressed this problem by utilizing different PSO topologies, or by integrating the PSO process with other optimization schemes such as local search or GAs. In this study, the quality of the solutions is improved by means of a so-called local search with savings (LSS) enhancement scheme which integrates a mutation local search procedure and a savings algorithm designed to minimize the route length for each vehicle.

The method proposed in this study offers several advantages over existing MDVRP methods. First, in contrast to existing methods, in which the locations of the multi-depots are taken as a given, the two-tier PSO scheme proposed in this study first optimizes the locations of the depots and then solves the MDVRP accordingly. Furthermore, to reduce the computational complexity, the proposed method utilizes a novel particle encoding scheme to solve all three sub-problems of the MDVRP in a single pass. Finally, the method improves the accuracy of the final solutions using a clustering and sweeping heuristic to set the initial search direction of the PSO algorithm and a LSS enhancement mechanism to minimize the routing cost.

The remainder of this article is organized as follows. Section “MDVRP” introduces the MDVRP considered in this study. Section “Two-tier PSO framework” describes the two-tier PSO framework used to solve the MDVRP and presents the proposed particle encoding scheme, the clustering and sweeping heuristic, and the LSS mechanism. Section “Experimental results” presents and discusses the experimental results. Finally, section “Conclusion” provides some brief concluding remarks and indicates the intended direction of future research.

MDVRP

The truck dispatching problem was first posed by Dantzig and Ramser ²⁶ in 1959 and subsequently led to the development of the MDVRP, in which a company has multi-depots from which to service its customers. The MDVRP involves three problems, namely, (1) allocating the customers to the depots; (2) assigning the customers to individual vehicles; and (3) determining the optimal routing sequence for each vehicle. Consider the illustrative example shown in Figure 1, in which a company has two depots (D₁ and D₂) and is required to service 26 customers (c₁, c₂, …, c₂₆). Assume that Depot 1 has four vehicles with which to service 12 customers, while Depot 2 has three vehicles to service 14 customers. The goal of the MDVRP is then to minimize the total traveling cost of the seven vehicles.

Figure 1.

Illustration of MDVRP.

Consider a general MDVRP involving n customers and m depots. Let the node set C  = {c₁, c₂, …, c_n, c_n + 1, …, c_n + m} be the combined set of customers and depots and let D  = {D₁, D₂, …, D_m} be the set of depots. In addition, let each depot have k vehicles, where the vehicle set is denoted as V  = {v₁, v₂, …, v_k}. Finally, let each customer i have a particular transportation demand represented by r_i. In solving the MDVRP, each depot dispatches its vehicles to deliver goods to the customers in such a way that each customer is visited only once and each vehicle starts from the depot and returns to the same depot on completion of its delivery run. The objective and constraints of the MDVRP can thus be formulated as follows

$$\mathrm{Fitness}=min\sum _{p=1}^m\sum _{v=1}^k\sum _{i=1}^{n+m}\sum _{j=1}^{n+m}{d}_{\mathrm{ij}}{x}_{\mathrm{ij}}^{\mathrm{pv}}i\ne j$$

(1)

$$x_{\mathrm{ij}}^{\mathrm{pv}}=\{\begin{array}{c}1,\mathrm{if}\mathrm{vehicle}v\mathrm{of}\mathrm{depot}p\mathrm{from}i\mathrm{to}j\hfill \\ \multicolumn{1}{c}{0,\mathrm{otherwise}}\end{array}$$

(2)

$$\sum _{i=1}^{n+m}\sum _{j=1}^{n+m}{x}_{\mathrm{ij}}^v{r}_{i\le Q_v}i\ne j,v\notin V$$

(3)

As shown in equation (1), in solving the MDVRP, the aim is to determine the shortest total routing distance of the vehicles. In other words, equation (1) defines the fitness function as the total traveling cost of all the vehicles of all the depots. Specifically, d_ij indicates the distance between node c_i and node c_j, and $x_{\mathrm{ij}}^{\mathrm{pv}}$ indicates whether or not vehicle v of depot p visits node c_i and then node c_j (see equation (2)). For instance, $x_{53}^{13}=1$ indicates that a certain vehicle (v₃) of depot D₁ services customers 3 and 5, with the vehicle first visiting Customer 5 then Customer 3, as shown in Figure 1. Equation (3) represents the capacity constraint of the MDVRP, where $x_{\mathrm{ij}}^{\mathrm{pv}}$ indicates that a certain depot p dispatches vehicle v to service customer i and then customer j, and Q_v denotes the maximum carry capacity of vehicle v. Notably, the total carry demands of the customers serviced by vehicle v must be less than (or equal to) the vehicle carry capacity. For example, for a vehicle route v₂: D₁ → c₁ → c₂ → c₄ → D₁, the sum of the carry demands r₁ + r₂ + r₄ must be no greater than Q_v2.

Two-tier PSO framework

This study proposes a two-tier PSO framework for solving the MDVRP in which an external PSO is first used to identify the optimal depot locations and an internal PSO is then applied to determine the optimal vehicle routing path between these depots and the customers.

PSO

PSO was first proposed by Kennedy and Eberhart ²⁷ in 1995. In the PSO procedure, the solutions are modeled as a particle swarm scattered in the solution space. The location of each particle represents a candidate solution to the problem, and, in each iteration, the particles move with a certain velocity to a new region of the solution space in order to search for new and improved solutions. In performing particle movement, three vectors are referenced, namely, the velocity vector of the original movement, the best location the particle has visited thus far (i.e. the individual optimal experience), and the best location visited by any of the particles in the swarm thus far (i.e. the global optimal experience). Suppose that n particles are scattered in an L-dimensional solution space. The particle swarm can thus be represented as X_i = {X₁, X₂, …, X_n}, where the location of each particle is given by X_ij = {X_i1, X_i2, …, X_iL}, in which X_ij refers to the jth location component of the ith particle. In addition, the velocity vector of each particle is represented by V_ij = {V_i1, …, V_iL}, the local experience of the particle by p_ij = {p_i1, …, p_iL}, and the global experience of the swarm by g = {g₁, …, g_L}. The PSO solution process is thus formulated as follows

$$\begin{array}{c}{v}_{\mathrm{ij}}^{\mathrm{new}}=v_{\mathrm{ij}}\times w+c_1\times r_1\times (p_{\mathrm{ij}}-x_{\mathrm{ij}})+c_2\times r_2\times (g_j-x_{\mathrm{ij}})\hfill \\ x_{\mathrm{ij}}^{\mathrm{new}}=x_{\mathrm{ij}}+v_{\mathrm{ij}}^{\mathrm{new}}\hfill \end{array}$$

(4)

Here, v^new refers to the updated velocity of the particle and x^new indicates the new location of the particle following movement. In addition, c₁ and c₂ are learning factors, and r₁ and r₂ are random numbers with values ranging from 0 to 1 and are used to generate diverse movement directions of the particles such that the particles explore new regions of the solution space.

Particle encoding schemes

The two-tier PSO framework proposed in this study utilizes two particle encoding schemes, namely, one scheme for the outer PSO and another scheme for the inner PSO. In the outer PSO, the particles represent candidate depot locations. The depot locations are expressed in the form of rectangular coordinates (x, y). Hence, the particles in the PSO swarm are encoded as a D × 2 matrix, $\mathbf{X}=\left[\begin{array}{c}\begin{array}{cccc}{x}_1\hfill & x_2\hfill & \cdots \hfill & x_D\hfill \end{array}\hfill \\ \multicolumn{1}{c}{\begin{array}{cccc}{y}_1\hfill & y_2\hfill & \cdots \hfill & y_D\hfill \end{array}}\end{array}\right]$ , in which x₁ refers to the x coordinate of the first depot, y₁ refers to the y coordinate of the first depot, and so forth.

Having obtained the optimal depot locations from the outer PSO, the inner PSO is applied to solve the corresponding MDVRP. In implementing the inner PSO, the particle encoding scheme proposed in this study is designed in such a way as to solve all three sub-problems of the MDVRP at the same time, that is, “which depot should offer service to which customers?”“which vehicle should deliver goods to which customers?” and “in what order should each vehicle visit its assigned customers?” Distributing the services of D depots over n customers involves dividing the customers into D groups. Furthermore, given m vehicles at each depot, the customers serviced by the depot are divided into m groups. In other words, the MDVRP involves dividing the customers into D × m groups. In this study, the MDVRP is thus solved using an extended particle encoding scheme ¹ which considers (D × m − 1) segmentation points in addition to n components corresponding to the n customers. In other words, each encoded particle string X consists of (n + D × m − 1) components, that is, X  = {x₁, …, x_n, x_n + 1, …, x_{n + D × m − 1}}.

In solving the MDVRP, the particles are decoded using a random key method. Therefore, in the encoding scheme, particle X corresponds to a set of key values, that is, KEY = {c₁, c₂, …, c_n, S₁₁, …, S_1m, …, S_D1, …, S_{D(m − 1)}}, in which S_ij denotes the key value of the segmentation point of the jth vehicle in depot i. The particle encoding scheme and key values are shown in Figure 2. In accordance with the random key decoding method, X is sorted in ascending order, and the corresponding decoding key values change accordingly. Figure 3 shows the decoding process for an illustrative MDVRP involving six customers (n = 6), two depots (D = 2), and two vehicles per depot (m = 2). The coded particle string thus comprises nine vectors (i.e. n + D × m − 1 = 6 + 2 × 2 − 1 = 9). Having sorted X in ascending order, the segmentation points define not only which vehicle in each depot visits which customers but also the sequence in which these customers are visited. For example, the jth vehicle in the ith depot, represented by segmentation point S_ij, offers service to all the customers before the segmentation point and visits the customers in the order indicated in the sorted sequence. Notably, the customers at the end of the sorted sequence (i.e. the customers after the final segmentation point) are simply serviced by the last vehicle of the last depot. Figure 4 shows the scheduling results obtained for the illustrative MDVRP. As shown, Customer 5 and Customer 1 are serviced by the first vehicle of the first depot in accordance with the following delivery route: Start from Depot → Service Customer 5 → Service Customer 1 → Return to Depot (i.e. Depot 1 → Customer 5 → Customer 1 → Depot 1). Similarly, Customer 3 and Customer 6 are serviced by the first vehicle of the second depot in accordance with the route: Depot 2 → Customer 3 → Customer 6 → Depot 2. Customer 2 is serviced by the second vehicle of Depot 1 in accordance with the route: Depot 1 → Customer 2 → Depot 1. Finally, Customer 4 is serviced by the second vehicle of the second depot in accordance with the route: Depot 2 → Customer 4 → Depot 2.

Figure 2.

PSO encoding and key values.

Figure 3.

Random key decoding.

Figure 4.

Scheduling results for illustrative MDVRP.

Initial solution heuristic mechanism for internal PSO

Random initial solutions cause the PSO algorithm to search the solution space in an aimless manner and hence the solution procedure converges very slowly. By contrast, good initial solutions limit the solution space which must be searched and therefore enable the solution procedure to converge more rapidly to the optimal solution. Accordingly, this study develops a heuristic mechanism for obtaining a certain proportion of effective initial solutions for the inner PSO. As described above, the MDVRP involves solving three sub-problems, namely, grouping the customers in accordance with the depots by which they are serviced, determining which vehicles at each depot should service which customers, and establishing the order in which each vehicle should visit its assigned customers. Therefore, in searching for suitable initial solutions for the MDVRP, the heuristic mechanism commences by grouping the customers using a method similar to that of the k-means clustering algorithm. Specifically, the customers are grouped using a set of weight values determined on the basis of their respective distances from each depot. The weight values are determined using the formulation shown in equation (5), in which w_ij is the weight of the distance between the ith customer and the jth depot, dist_ij is the distance between customer i and depot j, and dist_ik is the distance between customer i and depot k. In every case, a higher weight value indicates a shorter distance between the customer and the depot. For illustration purposes, consider the simple case of a single customer and two depots. Let the distances between the customer and the two depots be notionally assigned as 92.6 and 16.1, respectively. The corresponding weights are thus computed as w₁₁ = 1/(92.6 × (1/92.6 + 1/16.1)) = 0.148, w₁₂ = 1/(16.1 × (1/16.1 + 1/92.6)) = 0.852, w₁₁ + w₁₂ = 1

$$w_{\mathrm{ij}}=\frac{1}{\mathrm{dis}{t}_{\mathrm{ij}}\times \sum _{k=1}^D\frac{1}{\mathrm{dis}{t}_{\mathrm{ik}}}}$$

(5)

To achieve a more diverse initial grouping of the customers, the customers are assigned to particular depots using either a greedy method or a roulette method, where the choice of method is determined using a user-defined grouping probability (gp) parameter. More specifically, a random value rand is generated in the range of 0–1, and the greedy method is adopted if rand < gp; or the roulette method is adopted otherwise. In this study, gp is assigned a value of 0.8, that is, the greedy method is adopted with a probability of 80%, while the roulette method is adopted with a probability of 20%.

Having grouped the customers, a sweep algorithm is applied at each depot to screen the matching between the customers and the vehicles and the vehicle routing sequence, as shown in Figure 5. In executing the sweep algorithm, the polar coordinate angle of each customer relative to the depot is calculated as θ = tan⁻¹(y/x), where x and y are the location coordinates of the customer. The polar coordinate angles of all the customers are arranged in ascending order in order to generate a customer list. A random value is then used to choose one of the customers in the list as the initial customer to be serviced by the first vehicle. After assigning this customer to the vehicle, the next customer on the list is chosen and added to the same vehicle. The process is repeated in this way until the vehicle lacks sufficient capacity to service any more customers, at which point the next customer is assigned to the second vehicle at the depot. The process continues in this way until all of the customers have been serviced. For each vehicle, the customers are serviced in the order in which they were added to the vehicle.

Figure 5.

(a) Customer sweep and (b) customer-to-vehicle matching and corresponding vehicle routing.

As described above, if all of the initial solutions are generated randomly, the PSO algorithm will converge only very slowly to the optimal solution. Conversely, if all of the initial solutions are generated by the proposed heuristic method, the algorithm will converge rapidly, but may become trapped at a local sub-optimal solution. Accordingly, the initial heuristic rate (IHR) parameter shown in equation (6) is employed to determine the probability that the heuristic method proposed above is actually adopted when executing the inner PSO. Note that rand in equation (6) refers to a randomly generated value in the range of 0–1

$$\mathrm{initial}\mathrm{solution}=\{\begin{array}{c}\mathrm{heuristic}\mathrm{generation},\mathrm{if}\mathrm{rand}<\mathrm{IHR}\hfill \\ \multicolumn{1}{c}{\mathrm{heuristic}\mathrm{generation},\mathrm{otherwise}}\end{array}$$

(6)

Inertia weight adjustment in internal PSO

Ideally, the solution space should be searched in a wider manner at the beginning of the solution procedure, but more narrowly as the preferred region of the solution space becomes clear. In other words, the search process should gradually evolve from a global exploration-type search behavior to a local exploitation-type search behavior as the iteration process proceeds. Many linear and non-linear methods have been proposed for adjusting the inertia weight of the particles in PSO schemes. ²⁸ This study considers the linear method shown in equation (7), in which i is the current iteration number and i_max is the maximum number of iterations (i.e. the number of iterations at which the algorithm terminates)

$$w_i=(i_{\max }-i)/i_{\max }$$

(7)

Four non-linear methods are also considered, namely, the top chord decline curve (DC) method, ²⁹ the decline curve down (DCD) method, ²⁹ the S-decreasing curve (S-curve) method, ³⁰ and the sigmoid curve method. ³¹ In the initial period of the DC methods, the range of particle movement is narrowed quickly; later, it is narrowed more slowly. Throughout the iteration process, the method aims to provide a more small-range movement search. The DC has the form shown in equation (8), where u controls the curvature of the top chord of the curve

$$w_i=(1-\frac{i}{i_{\max }})/(1+u\times \left(\frac{i}{i_{\max }}\right))$$

(8)

Enhancement mechanism for internal PSO

To enhance the quality of their solutions, movement-based algorithms generally undertake a local search in the neighborhood of each newly-obtained solution to see whether or not a better solution can be found. Accordingly, the internal PSO proposed in this study also performs a local search each time it obtains a new solution. If the new solution is found to be better than the current one, the current solution is replaced with the new one; else, the current solution is retained. The quality of the solutions is further improved by applying the savings algorithm proposed in Lysgaard ³² to the customer visit sequence obtained in the first stage of the solution process with the hope of finding a new customer visit sequence which further reduces the length of the vehicle routing path.

In implementing the enhancement mechanism, seven local search methods are considered, namely, a developed mutation method, and six methods presented in the literature, that is, swap, ³³ section swap, ³³ reverse section swap, ³³ reverse, ³³ insertion, ³³ and double swap. ³⁴ In the mutation method, random values (rand) ranging from 0 to 1 are generated for all of the components in the particle code, and these values are then compared with a mutation rate (mr) parameter. If rand < mr, another random value (δ) ranging from −r to r is generated and added to the component; else, the component remains unchanged, as shown in equation (9). Note that mr and r are determined through testing

$$x_i=\{\begin{array}{c}{x}_i+\delta ,\mathrm{if}\mathrm{rand}<\mathrm{mr}\hfill \\ \multicolumn{1}{c}{x_i,\mathrm{otherwise}}\end{array}$$

(9)

The savings algorithm commences by calculating the cost of visiting each customer separately and then calculates the cost incurred in visiting all of the customers in accordance with the determined visit sequence. The reduction in cost is referred to as the savings value. Consider the simple illustrative example shown in Figure 6, in which 0 represents the depot, and a and b represent two customers (Customer a and Customer b, respectively). The main steps in the savings algorithm can be described as follows.

Figure 6.

Illustrative example of savings algorithm.

1. Calculate the savings value S_ab between Customer a and Customer b using equation (10)

$$\begin{array}{c}\multicolumn{1}{c}{S_{\mathrm{ab}}=(c_{0a}+c_{a0}+c_{0b}+c_{b0})}\\ -(c_{0a}+c_{b0}+c_{\mathrm{ab}})=c_{a0}+c_{0b}-c_{\mathrm{ab}}\hfill \end{array}$$

(10)

Note that for n customers, the number of savings values obtained from equation (10) is equal to n × (n − 1)/2 (see Table 2).

Table 2.

Savings values S_ab between Customer a and Customer b.

a	b
	1	2	3	4
1	–	14	9	7
2		–	18	11
3			–	15
4				–

• 2. Sort the savings values in descending order and identify the corresponding customer visit sequence (see Table 3).

• 3. If two customers in the list are neighbors, connect them; else, ignore the connection (see Figure 7).

Table 3.

Sorted savings values.

Savings value	18	15	14	11	9	7
Visit sequence	2–3	3–4	1–2	2–4	1–3	1–4

Figure 7.

Result of savings algorithm for illustrative case.

In this study, the combined local search and savings algorithm is denoted simply as the LSS enhancement mechanism. The overall operation of the proposed two-tier PSO scheme for solving the MDVRP is shown schematically in Figure 8.

Figure 8.

Flowchart of proposed two-tier PSO scheme for solving MDVRP.

Experimental results

The experiments commenced by identifying the optimal values of the inertia weight and IHR parameters in the internal PSO. Further experiments were then performed to evaluate the effectiveness of the LSS enhancement mechanism. The experiments were performed using MDVRPs developed by Cordeau (http://neo.lcc.uma.es/vrp/vrp-instances/multiple-depot-vrp-instances/) as test instances. Six cases were chosen for test purposes (two large, two average, and two small), denoted as p01, p02, p03, p04, p05, and p06, respectively. The effectiveness of the external PSO in identifying the optimal locations of the multiple depots was investigated by comparing the total vehicle routing distance obtained by the internal PSO on the basis of these locations with that obtained by the internal PSO only (i.e. the depot locations were not optimized in advance). The internal PSO was implemented using 40 particles, learning factors (c₁ and c₂) equal to 2, and a maximum number of iterations equal to 2500. The performance of the internal PSO scheme was evaluated using the deviation (Dev) and average deviation (Dev.Avg) metrics shown in equation (11), where the best known solution is the optimal known solution in the corresponding test set, the instance number is 6, and N is the number of trials and has a value of 10 in this study

$$\begin{array}{l}Dev(\%)={\displaystyle \sum _{i\in \mathrm{instances}}\left(\frac{obtainedsolutio{n}_i-bestknownsolutio{n}_i}{bestknownsolutio{n}_i}\times 100\%\right)}/\left|\mathrm{instances}\right|\\ Dev.Avg(\%)=\frac{{\displaystyle \sum _{t=1~T}Dev(\%)}}{N}\end{array}$$

(11)

The experiments were conducted in a Windows 7 SP1 environment on a PC with an Intel Core i7 CPU 4770 (3.40 GHz), 8 GB of memory, and Visual Studio 2010 C# software. As described above, the experiments commenced by examining the effects of the inertia weight and IHR parameter settings on the optimization performance of the internal PSO scheme. Table 4 shows the experimental results (Dev.Avg) obtained when using various fixed inertia weights (ω) in the range of 0.1–1.0 and different values of the IHR parameter. As shown, the optimal performance (i.e. the lowest Dev.Avg) is obtained using an inertia weight of ω = 0.3 and an IHR parameter setting of IHR = 0.3. In other words, 30% of the initial solutions for the internal PSO are generated using the heuristic method, while the remaining solutions are generated randomly. In searching a solution space, it is desirable to commence by performing large-range movements (i.e. exploitation) in the initial stage of the search process and then to shift gradually toward small-range movements (i.e. exploration) as the search process proceeds. Accordingly, as described in section “Inertia weight adjustment in internal PSO,” various linear and non-linear adjustment policies were applied to update the inertia weight parameter during the internal PSO search procedure. Table 5 shows the Dev.Avg results obtained when using the different adjustment policies. It is seen that the optimal performance is obtained using the DC method with a reduction value range of 0.6–0.2.

Table 4.

Dev.Avg results obtained using different inertia weights and IHR parameter settings in internal PSO.

IHR	Inertia weight (ω)
	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
0	142.31	133.93	134.01	139.35	144.63	155.90	N/A	N/A	N/A	N/A
0.1	17.38	19.81	23.75	16.22	19.81	17.71	49.54	47.11	49.80	53.76
0.2	17.78	15.55	16.09	16.39	16.44	19.08	46.63	44.33	49.54	49.48
0.3	17.42	16.34	14.95	16.23	16.65	18.53	46.21	45.27	49.05	49.47
0.4	18.11	16.45	17.04	16.58	16.68	18.75	45.69	44.62	49.28	49.42
0.5	18.28	16.14	17.30	16.98	17.92	19.10	45.47	44.15	49.16	48.25
0.6	19.43	17.04	18.09	17.92	18.55	21.25	44.62	44.32	48.68	48.49
0.7	18.80	17.63	20.13	20.04	17.77	21.43	45.37	44.75	48.47	48.33
0.8	20.62	21.13	21.84	20.73	20.61	23.76	45.44	44.33	48.55	48.54
0.9	25.98	26.77	31.92	30.09	28.11	26.15	45.61	44.26	47.65	48.34
1	46.96	47.64	46.43	46.63	46.75	46.33	45.68	43.30	48.45	47.81

IHR: initial heuristic rate; PSO: particle swarm optimization.

Table 5.

Dev.Avg results obtained using different inertia weight adjustment policies in internal PSO.

Inertia weight	Decreasing curves
	DC	DCD	Linear	S-curve	Sigmoid
0.5–0.1	16.86	16.56	15.57	16.45	16.36
0.5–0.2	16.67	16.85	16.68	15.56	15.85
0.6–0.1	16.05	16.92	16.21	15.62	16.78
0.6–0.2	14.84	16.73	16.6	16.64	17.07
0.9–0.4	17.1	32.25	21.08	18.92	19.34

PSO: particle swarm optimization; DC: decline curve; DCD: decline curve down.

Overall, the results presented in Table 4 show that a larger inertia weight increases the deviation of the obtained solution from the optimal solution reported in the literature. This result is reasonable since a larger inertia weight results in a larger particle search step size and hence increases the risk of the solution procedure converging toward a local minimum. Moreover, parameter settings IHR = 0 and IHR = 1.0 both lead to a poor optimization performance. Again this finding is to be expected since if all of the initial solutions are generated randomly (IHR = 0) or by the heuristic (IHR = 1.0), that is, the initial solutions are too diverse or lack diversity, respectively, the PSO procedure either converges too slowly or converges rapidly to a local sub-optimal solution. As commented above, the non-linear top chord DC method with a reduction range of 0.6–0.2 yields the best PSO outcome. In other words, the optimization performance is improved when the range of particle movement is narrowed quickly in the initial stages of the search process and is then narrowed more slowly toward the later stages.

Adopting an IHR parameter value setting of 0.3, and either a fixed inertia weight setting of ω = 0.3 or the DC inertia weight adjustment policy with a reduction value ranging from 0.6–0.2, a further series of experiments was performed to identify the optimal settings of the mutation rate (mr) and random value parameter (δ) in the LSS enhancement mechanism of the internal PSO. Table 6 shows the Dev.Avg results obtained when using different mr and δ adjustment policies in the internal PSO and a constant inertia weight of ω = 0.3. Table 7 shows the Dev.Avg results obtained when using decreasing inertia weights (ω) in the range of 0.6–0.2 in the internal PSO. As shown, for both weighing scenarios, the optimal routing performance is obtained using a mutation rate of mr = 0.01 and a random value parameter of r = 3 (δ, −3 to 3). In other words, the PSO optimization performance is improved using a low mutation rate of mr = 0.01 (which prevents a random walk in the solution process) and a higher perturbation strength of r = 3 (which provides the iteration process with sufficient momentum to escape local sub-optimal solutions).

Table 6.

Dev.Avg results obtained using different mr and δ (IW = 0.3, IHR = 0.3).

r	mr
	0.01	0.02	0.03	0.04	0.05	0.06	0.07	0.08	0.09	0.1
0.25	14.66	14.34	15.18	15.74	14.86	14.50	16.61	15.48	15.67	15.59
0.5	13.48	14.77	15.19	15.92	14.25	14.74	16.04	15.18	17.68	17.18
0.75	14.18	13.64	14.59	15.04	14.29	15.83	16.59	16.01	16.55	18.06
1	13.90	14.14	14.92	14.80	16.13	15.28	15.67	17.11	17.20	17.67
1.25	13.97	13.83	14.24	14.38	14.92	16.78	15.94	16.81	18.94	18.77
1.5	13.51	13.23	14.49	14.74	14.67	15.40	18.27	20.00	18.44	18.47
1.75	12.20	12.81	12.76	13.65	16.02	16.31	17.17	17.92	19.12	19.59
2	13.20	13.92	14.22	14.14	16.35	18.23	16.92	18.90	20.68	19.87
3	12.08	12.28	13.91	15.85	16.35	17.01	18.27	19.36	21.02	20.86
4	12.11	12.95	15.30	14.60	16.04	17.06	18.85	19.34	19.84	22.22

Table 7.

Dev.Avg results obtained using different mr and δ (IW = 0.6–0.2, IHR = 0.3).

r	mr
	0.01	0.02	0.03	0.04	0.05	0.06	0.07	0.08	0.09	0.1
0.25	15.47	15.76	15.78	14.90	16.08	16.98	17.35	17.23	16.60	17.80
0.5	14.93	14.58	17.08	17.00	14.70	16.02	17.49	16.54	17.81	17.29
0.75	14.44	13.42	15.02	15.54	15.17	15.35	16.41	16.75	16.75	16.84
1	15.06	14.81	15.85	14.58	14.53	17.21	16.91	18.16	17.30	19.89
1.25	14.20	14.58	15.44	14.78	15.31	16.81	16.78	17.21	18.25	19.32
1.5	13.79	13.90	14.28	14.76	16.21	17.43	16.67	18.37	18.58	22.53
1.75	12.71	13.40	14.71	15.74	16.34	16.74	18.82	17.65	19.72	20.12
2	14.69	14.02	14.27	15.42	16.16	17.68	18.56	18.92	21.07	20.10
3	12.13	14.41	14.18	14.19	15.68	19.22	19.10	21.91	19.50	20.27
4	13.15	15.34	15.53	16.84	17.23	19.09	20.03	19.78	22.17	23.37

The performance of the LSS enhancement scheme was evaluated for each of the different local search methods (i.e. swap, section swap, reverse section swap, reverse, mutation, insertion, and double swap). Among them, the LSS enhancement scheme involving mutation yields the best performance. As in the previous set of experiments, the performance of the LSS scheme was evaluated using both a constant inertia weight setting (IW = 0.3) and a reducing weight setting (DC, 0.6–0.2). In addition, the experiments were performed with and without the savings algorithm, respectively. The IHR parameter was set as IHR = 0.3 in every case. The Dev.Avg results obtained using the two different weight adjustment policies are presented in Tables 8 and 9, respectively. It is seen that for both adjustment policies, the best routing performance is obtained using the savings algorithm in conjunction with the mutation local search method.

Table 8.

Dev.Avg results obtained using different local search schemes (IW = 0.3, IHR = 0.3).

	Ls1	Ls2	Ls3	Ls4	Ls5	Ls6	Ls7
Without savings	15.11	16.91	17.32	15.77	12.08	15.61	16.94
With savings	12.55	15.17	16.24	16.26	11.08	16.93	15.46

Table 9.

Dev.Avg results obtained using different local search schemes (DC curve 0.6–0.2, IHR = 0.3).

	Ls1	Ls2	Ls3	Ls4	Ls5	Ls6	Ls7
Without savings	15.24	16.04	16.12	17.13	12.13	16.52	15.23
With savings	12.31	18.03	17.49	16.52	11.78	18.04	15.02

DC: decline curve.

The effectiveness of the proposed two-tier PSO scheme was evaluated by comparing the routing performance obtained by the internal PSO using the optimal depot locations identified by the external PSO with that obtained using the internal PSO alone based on non-optimized depot locations. In implementing the external PSO, the number of particles was set as 10, the particle inertia weight was set equal to ω = 0.8, and the maximum number of iterations was given as 300. In the internal PSO, the particle inertia weight and IHR were set as ω = 0.3 and IHR = 0.3, respectively, and the full LSS enhancement mechanism was deployed. The initial solutions of the depot locations were generated randomly in accordance with the customer distribution range. Table 10 presents the experimental results (total routing distance) obtained with and without the external PSO for the six MDVRP test cases. Note that Best1 and Avg1 refer to the best (i.e. shortest) vehicle routing distance obtained from the internal PSO only over 10 repeated tests and the corresponding average vehicle routing distance over the 10 tests, respectively. Best2 and Avg2 give the equivalent results for the full two-tier PSO scheme (i.e. the external PSO and the internal PSO). Finally, the “Improvement” values indicate the percentage reductions in the best and average routing distances, respectively, when the external PSO is employed in addition to the internal PSO (i.e. Best (%) = ((Best2 − Best1)/Best1) × 100% and Avg (%) = ((Avg2 − Avg1)/Avg1) × 100%). It is seen that for each of the considered MDVRPs, the use of the external PSO to optimize the depot locations improves both the BEST and the AVG routing distance. Overall, the results show that the external PSO reduces the routing distance by 13.16% on average compared to the case where the depot locations are not optimized in advance.

Table 10.

BEST and AVG routing distances obtained using internal PSO only and full two-tier PSO framework.

Case	Without outer PSO		With outer PSO		Improvement
	Best1	Avg1	Best2	Avg2	Best (%)	Avg (%)
p01	604.23	625.24	539.11	544.39	−10.78	−12.93
p02	507.96	536.46	458.39	463.07	−9.76	−13.68
p03	705.30	738.57	632.91	640.87	−10.26	−13.23
p04	1108.90	1147.25	1033.32	1040.41	−6.82	−9.31
p05	865.46	977.96	757.56	768.19	−12.47	−21.32
p06	962.40	992.05	902.40	907.90	−6.23	−8.48
Avg					−9.39	−13.16

PSO: particle swarm optimization.

To better visualize the dynamic convergence performance of the proposed two-tier PSO scheme, Figures 9 and 10 show the evolution of the PSO solutions for the vehicle routings and depot locations, respectively, over the course of the respective solution procedures. In particular, Figure 9(a)–(d) shows the vehicle routings obtained by the internal PSO in the initial state and after 500, 1500, and 3000 iterations, respectively, given the p01 test case with four depots. Figure 10(a)–(d) shows the evolution of the multi-depot locations during the external PSO solution process. Note that in each sub-figure in Figure 10, the solid circular symbols labeled “g” indicate the best depot locations identified by the swarm so far. Figure 11(a) and (b) shows the fitness evolutions of the external PSO and internal PSO, respectively. The validity of the proposed two-tier PSO scheme was further evaluated by considering a real-world example involving a wholesaler in Taiwan having two depots from which to service 20 retailers using two vehicles per depot. Based on an inspection of the wholesaler’s records, the capacity of each vehicle was set as 150, and the demand of each customer was set in the range of 10–30. The company reported that the original locations of the two depots were not well planned and therefore resulted in a significant delivery time and cost. Accordingly, the proposed PSO framework was employed to replan the depot locations and reschedule the delivery routes accordingly. The two new depots have to be determined to reduce the traveling costs. The two new depots need to send out two trucks to carry all demand needed from all customers and back to the depot. Figure 12(a) and (b) shows the original depot locations and vehicle routes and those determined using the proposed PSO framework, respectively. A detailed inspection of the results shows that the proposed scheme reduces the total vehicle routing distance from 227 to 186 km, representing a cost saving of 18%.

Figure 9.

Vehicle routings obtained from internal PSO for p01 test case: (a) initial, (b) 500 iterations, (c) 1500 iterations, and (d) 3000 iterations.

Figure 10.

Depot locations obtained from external PSO for p01 test case: (a) initial, (b) 50 iterations, (c) 150 iterations, and (d) 300 iterations.

Figure 11.

Fitness evolution of external and internal PSOs: (a) internal PSO and (b) external PSO.

Figure 12.

Real-world case example: (a) original depot locations and vehicle routes and (b) optimized depot locations and vehicle routes.

Conclusion

This study has proposed a two-tier PSO framework consisting of an external PSO and an internal PSO for solving the MDVRP. In the proposed framework, the external PSO is used to determine the optimal depot locations based on the geographic distribution of the customers. These depot locations are then used by the internal PSO to establish the shortest vehicle routing path. The internal PSO is implemented using a novel particle encoding scheme which enables the simultaneous solution of all three MDVRP sub-problems, namely, “which depot should offer service to which customers?”“which vehicle should deliver goods to which customers?” and “in what order should each vehicle visit its assigned customers?” Consequently, the computational complexity of the MDVRP solution process is significantly reduced compared to existing schemes in which each sub-problem is solved individually. The performance of the proposed internal PSO is further improved through the use of a heuristic mechanism to generate suitable initial solutions for the search process and an LSS enhancement mechanism consisting of a local search process and a savings algorithm. The performance of the proposed two-tier PSO scheme has been evaluated using six test instances and a real-world case in Taiwan involving 2 depots, 20 customers and 4 vehicles.

The experimental results have shown that when solving the MDVRP using the internal PSO only (i.e. the depot locations are not optimized in advance), the optimal performance is obtained when adjusting the particle inertia weights using a top chord DC method (0.6–0.02), setting the heuristic mechanism parameter to IHR = 0.3, and applying the LSS enhancement mechanism with a mutation local search method and a savings algorithm. In addition, it has been shown that the use of the external PSO to optimize the multi-depot locations reduces the average total routing distance by 13.16% on average over the six considered test cases. Moreover, in the real-world case study, the two-tier PSO scheme reduces the total vehicle routing distance by around 18%. Overall, the experimental results show that the proposed two-tier PSO framework provides a computationally efficient and effective means of solving practical MDVRPs. In particular, the results confirm the importance of optimizing the depot locations prior to solving the MDVRP.

This study has assumed that the vehicles provide a delivery service only. However, in real-world situations, it is often necessary for the vehicles to provide both a delivery service and a pick-up service. This has significant implications for both the vehicle routing planning process and the capacity constraint applied in the optimization process. Hence, future studies will extend the proposed two-tier PSO scheme to the more realistic case of MDVRPs with both pick-up and delivery requirements. In addition, this study has considered the VRP to be a symmetric VRP, which can be treated as a complete undirected graph problem. In other words, it has been assumed that in making deliveries, each vehicle departs from and returns to the same depot. However, in practical cases, the vehicles may actually start from and return to different depots. Consequently, future studies will extend the proposed PSO framework to the more realistic case of asymmetric vehicle route planning.