A Hybrid Heuristic Based on a Particle Swarm Algorithm to Solve the Capacitated Location-Routing Problem With Fuzzy Demands

In this paper, the capacitated location-routing problem with fuzzy demands (CLRP-FD) is considered, which simultaneously solves two problems: locating the facilities and designing the vehicle routes among the established facilities and customers. In the CLRP-FD, the capacities of the employed vehicles and established facilities cannot exceed, and the demands of the customers are assumed to be triangular fuzzy variables. To model the CLRP-FD, a fuzzy chance constrained programming approach is designed using fuzzy credibility theory. To solve this problem, a hybrid particle swarm optimization (HPSO) algorithm, which includes a stochastic simulation and a local search strategy based on the variable neighborhood search algorithm, is proposed. Finally, the influence of the value of the dispatcher preference index (DPI) on the total distribution cost is analyzed by conducting numerical experiments. To evaluate the efficiency of the proposed HPSO and the performance of the CLRP-FD model, the results obtained using the HPSO were compared with their corresponding lower bound provided using the CPLEX solver. Moreover, we also evaluated the performance of HPSO through computational experiments on some well-known benchmark CLRP instances. The numerical results show that the proposed HPSO is competitive, and can give a satisfactory solution in a reasonable amount of time.


I. INTRODUCTION
In modern logistics decisions, facility locations and vehicle routing designs are vitally important to supply chain management because the dependable, efficient, and flexible decisions on the depots and the delivery routes can not only save distribution cost and time but also enhance the company's ability to compete. Thus, a large number of researchers have been devoted to the facility location problem (FLP) [1] and vehicle routing problem (VRP) [2]- [4]. FLP is an example of a strategic level decision in which a long-term decision is made to select the locations of the facilities. VRP is an example of a tactical level decision in which vehicle routes are designed according to the demands of the customers and The associate editor coordinating the review of this manuscript and approving it for publication was Jenny Mahoney. other factors. However, many researchers have indicated that solving of the FLP without considering the VRP could be immoderate and lead to a suboptimal solution [5], [6], and vice versa. As a combination of FLP and VRP, the location routing problem (LRP) simultaneously includes the strategic and tactical levels to overcome this drawback [7] and is conductive to building a flexible and efficient logistics system.
Since both FLP and VRP are well-known NP-hard combinatorial problems, the LRP must necessarily be an NP-hard problem [8], [9]. Overall, when the LRP instance size increases, the computing time increases exponentially. Thus, exact methods, such as the Branch-and-Cut method [5], [9] and Branch-and-Bound algorithm, cannot optimally solve large-size LRP instances in a reasonable amount of computing time. On the other hand, the wide range of applications of the LRP, such as waste collection [9], [10], newspaper delivery, bill delivery [11], blood shipments [12] and perishable food distribution [13], often leads to large-size instances. Under such conditions, many researchers focus their attentions on heuristics. Heuristics have been applied to many fields. Approaches such as simulated annealing heuristic [16]- [18], particle swarm optimization [19], genetic algorithm [20], ant colony algorithm [21], and hybrid artificial intelligence and robust optimization [22] have been used to solve the LRP and other supply chain management problems [14], [15], and they can provide suboptimal solutions.
The capacitated location-routing problem (CLRP), an important variant of the standard LRP, is constrained by the capacities of the vehicles and depots. CLRP is much closer to reality than the standard LRP and has been attracting more and more attention worldwide from researchers. The aim of the CLRP is to optimally determine the locations of the depots, assign customers to each established depot as well as detect a set of distribution routes for a set of identical vehicles in such a way that all customers with known demands can be serviced by one vehicle exactly once. Moreover, every vehicle should start and end its route at the same given depot. Due to the difficulty of the CLRP, exact methods are usually not good at solving intermediate and large instances. Therefore, heuristic and meta-heuristic approaches, which can find approximate solutions in a reasonable amount of computing time, are suitable methods for solving intermediate and large CLRP instances [3]. In addition to the basic characteristics of each solution method, the corresponding experimental results are also illustrated and studied in the literature. Nagy and Salhi surveyed the meta-heuristics and categorized them into four main groups, namely, sequential, clustering, iterative and hierarchical methods [3], [23]. Interested researchers can refer to the review paper [23] for more details about the mathematical formulations and solution methods for the LRP.
The CLRP requires that the demand of each customer should be known clearly before distribution. However, in most of the real-world situations, the demands of the customers can be vague or uncertain. These uncertain customers' demands are considered in the capacitated location routing problem with fuzzy demands (CLRP-FD) [24]. In CLRP-FD, the demand of each customer is assumed to be imprecise, but according to experience, it can be approximately a given amount or within a given range, for example, approximately 40 units, or between 30 and 80 units. For these cases, fuzzy logic can be worthwhile, and a fuzzy chance constrained programming model can be formulated based on fuzzy credibility theory. The CLRP-FD was first modeled and solved by Mehrjerdi and Nadizadeh in 2013 [24]. In [24], Mehrjerdi and Nadizadeh presented a fuzzy chance constrained programming approach with a credibility measurement to model the CLRP-FD, which was then solved by a greedy clustering method with a stochastic simulation. Furthermore, Nadizadeh and Nasab considered a dynamic CLRP-FD, in which the facility location problem and vehicle routing problem are solved on a time horizon [21]. To solve the dynamic CLRP-FD formulation, they developed a hybrid heuristic algorithm that had four phases, including the stochastic simulation and a local search method [21]. In addition, Nadizadeh and Kafash also studied the fuzzy capacitated location-routing problem with simultaneous pickup and delivery (FCLRP-SPD), in which both pickup and delivery demands of customers were considered as fuzzy variables, and developed a greedy clustering method with four phases to solve this problem [25]. Without considering the capacities of the depots, Ghaffari-Nasab et al. considered an LRP with fuzzy demands and proposed a fuzzy chance-constrained programming model based on fuzzy credibility theory [17]. They used a hybrid simulated annealing-based heuristic to solve this problem, in which stochastic simulation was employed to assess the believability of a solution. Zarandi et al. proposed a CLRP with time windows under uncertainty, in which the demands of the customers and the travel time were assumed to be fuzzy variables [26]. They also developed a simulation-embedded simulated annealing algorithm [26] to solve this problem. Fazayeli et al. developed a mixed-integer nonlinear mathematical fuzzy model for the LRP in a multimodal transportation network with time windows and fuzzy demands, in which the capacities of the depots were not considered [27]. Moreover, they proposed a genetic algorithm with two parts and a new chromosomes structure [27] to solve this model.
To the authors' knowledge, the CLRP-FD has significantly fewer investigations, and there is still plentiful space for the improvement of the corresponding model and solution methods. In this paper, the CLRP-FD is also considered, which is formulated using the approach of fuzzy chance constrained programming with credibility measurement. In this problem, the actual demand of a customer can be known only when the vehicle arrives at the customer [24], [28]. Thus, an additional distance arises when the vehicles implement preplanned distribution routes, which are designed using a given solution approach. Owing to the uncertain demands of the customers, a customer might not be served by a vehicle that arrives there and does not have sufficient capacity. In such a situation, the vehicle is assumed to return to the depot for loading, then go back to the customer and resume its service along the remainder of the preplanned route [24]. This strategy results in an additional distance, which also greatly impacts the distribution cost and efficiency. Thus, the additional distance is also minimized in this work when minimizing the locating and routing costs.
To solve the proposed CLRP-FD model, a hybrid particle swarm optimization (HPSO) algorithm with three phases, including a local search method and stochastic simulation, is presented. The motivation of using HPSO in this work is threefold. First, the PSO, a swarm-based meta-heuristic algorithm, has provided good solutions for many combinatorial optimization problems (e.g., multidimensional knapsack problem [29], data allocation problem [30] and personnel assignment problem [31]) within a reasonable amount of computing time, especially many different variants of the FLP, LRP and VRP have been solved using PSO-based approaches in the literature (e.g. VRP with stochastic demands [32], VRP with dynamic requests [2], multi-workshop facility layout problem [33] and LRP with stochastic demands [19]). Second, a meta-heuristic is a general algorithmic framework for finding solutions to optimization problems within which local heuristics can be guided and adapted to effectively explore the solution space [3]. Third, to obtain better results, there should still be some opportunities to improve the PSO by introducing some heuristic techniques into it. The authors are intensively motivated by these factors to use the hybrid PSO to solve the proposed CLRP-FD model. This paper makes the following contributions. First, a chance-constrained programming model of CLRP with fuzzy demands is proposed. Second, a hybrid particle swarm optimization algorithm with a stochastic simulation is presented for solving the CLRP-FD, in which an effective initialization method based on the greedy cluster algorithm and a local search method based on the path-relinking algorithm and variable neighborhood search algorithm [34] are employed to provide good initial solutions to the algorithm and strengthen the exploitation in the promising solution space, respectively. Finally, some experiments are conducted to evaluate the operational gains obtained using the proposed CLRP-FD model and to analyze the performance of the proposed HPSO algorithm. As shown later, the proposed algorithm can effectively produce good solutions to the CLRP-FD.
The rest of this paper is organized as follows: Fuzzy credibility theory and some assumptions are introduced in the section II. The details of CLRP model with fuzzy demands are stated in the section III. The proposed HPSO algorithm to solve the CLRP-FD model is discussed in the section IV. Then, computational results and analysis are given in the section V. Finally, conclusions and some future researches are presented in the section VI.

II. FUZZY CREDIBILITY THEORY
The concept of a fuzzy set was first proposed by Zadeh [35]. To measure a fuzzy event, Kaufmann [36] and Zadeh [37] successively proposed the term fuzzy variable and the corresponding possibility measure theory, respectively. Although possibility measure theory has been widely used, it does not obey the law of truth conservation and is inconsistent with the laws of excluded middle and contradiction [38]. Because a self-dual measure is very necessary and important to the possibility measure, Liu [39] proposed the concept of a credibility measure for defining self-duality. In this study, the CLRP with fuzzy demand is formulated based on a credibility measure. The related concepts and definitions used in this paper are briefly introduced as follows: Let and φ be non-empty and empty set, respectively. P is the power set of (P{ }), and each element in P, denoted as A, is called an event. Pos{A} is the occurring possibility of event A. The following four axioms are the basic of credibility theory, and all the concepts of credibility theory can be obtained from them [39].
Axiom 4: If i is a non-empty set,and the set function Pos i {}(i = 1, 2, · · · , n) satisfies above three axioms, and = 1 × 2 × · · · × n , then for each A ∈ P( ), Note that Pos = Pos 1 ∧ Pos 2 ∧ · · · Pos n . If Pos satisfies the first three axioms, it is named possibility measure. The triple set ( , P( ), Pos) is a possibility space.With the axiom 4, the basic of credibility theory was built.
Definition 1: Let ( , P( ), Pos) be a possibility space, and A is a set in P( ), then the necessity measure and credibility measure of A are defined by Nec{A} = 1 − Pos{A c } (A c is the complement of set A) and Cr{A} = 1 2 (Pos(A)+Nec(A)), respectively.
Definition 1 shows that the credibility of a fuzzy event is the average of its possibility and necessity [24]. The credibility measure is self-dual [24]. A fuzzy event may fail even though its possibility is equal to 1, and occur even though its necessity is 0 [24]. However, the fuzzy event must occur if its credibility is equal to 1, and fail if its credibility is 0 [39]. Now, a triangular fuzzy variabled = (d 1 , d 2 , d 3 ) is considered as the demand of a customer whered is described by its left boundary d 1 and its right boundary d 3 . The dispatchers can subjectively estimate the value of d 1 and d 3 based on their experience and available information, and the actual demand of each customer will be between d 1 and d 3 . The value of d 2 corresponding to a grade of membership of 1 can also be subjectively estimated. The membership function is shown as Figure 1. Let r be the actual demand of a customer, the possibility Pos(d ≥ r), necessity Nec(d ≥ r) and credibility Cr(d ≥ r) are easy to obtain as follows [40]: VOLUME 8, 2020

III. THE CLRP-FD MODEL
In this section, the fuzzy chance-constrained programming model of CLRP-FD is presented. The following assumptions were considered for it: • Each vehicle and potential depot have their distinct capacities.
• Each vehicle can be assigned to one and only one distribution route at most.
• The demand of each customer can be distributed once by exactly one vehicle at most.
• Each vehicle starts and ends at a single depot.
• For each vehicle, the total demand distributed by it cannot exceed its capacity.
• For each depot, the total demand of all customers allocated to it cannot exceed its capacity.
• The demand of customer i can be estimated in triangular Let Q be the capacity of a vehicle. It is easy to know that the remaining capacity of the vehicle is Q h = Q − h i=1d i after distributing the demands of h customers. Based on the rules of fuzzy arithmetic, it is easy to deduce that Q h , is also a triangular fuzzy number. Then, the credibility that the demand of the h + 1th customer (d h+1 ) is not greater than Q h can be given as follows: Let P be the capacity of a depot. Similarly, after allocating h customers to a depot, its remaining capacity equals , which is also a triangular fuzzy number. Then, the credibility that the demand of the h + 1th customer (d h+1 ) is not greater than P h can be given as follows: Obviously, the possibility of serving the next customer is greater if the remaining capacity of the vehicle is high and the demand of the next customer is low [21]. The value of Cr (Cr ∈ [0, 1]), calculated using formula (5), shows the level of certainty about the ability of a vehicle to serve a new customer [26]. When Cr = 0, the vehicle cannot serve the next customer and should return to the depot [26]. When Cr = 1, the dispatcher will definitely make sure that the remaining vehicle capacity can serve the next customer. In this study, a threshold is defined as DPI, where DPI ∈ [0, 1] is the dispatcher preference index, which can affect the solution considerably [26]. We decide whether the vehicle serves the next customer by comparing the value of DPI with the credibility obtained using the equation (5). The vehicle should serve the h + 1th customer if Cr ≥ DPI ; otherwise, it should return to the depot to load sufficient goods and, then, return to the customer. The process does not terminate until all of the customers are serviced.
Similarly, the value of Cr (Cr ∈ [0, 1]), calculated using formula (7), shows the level of certainty with which a depot is able to serve a new customer [26]. If the remaining capacity of the depot is high and the demand of the next customer is low [21], the possibility of allocating the new customer to the depot is greater. When Cr = 1, the manager of the depot is absolutely certain that the remaining capacity of the depot can provide service to the next customer. When Cr = 0, the manager of the depot will definitely ensure that the next customer cannot be allocated to the depot due to insufficient available capacity. In this study, another threshold is defined as API , where API ∈ [0, 1] is the manager preference index, which can also affect the solution considerably. The manager decides whether the depot serves the next customer by comparing the value of Cr with the value of API . The h + 1th customer will be allocated to the current depot if Cr ≥ API ; otherwise, he/she will be assigned to another open depot. This process does not terminate until all of the customers are allocated [21].
In this paper, an additional distance is considered, which is produced by the ''failure'' of the planned route and has non-neglected impact on the routing cost. As stated above, the parameters DPI and API can considerably affect the length of the planned routes and the additional distance. For example, a lower value of DPI , which shows that the dispatcher wants to use the the capacity of the vehicle the best he can, increases the possibility of failure that a vehicle arrives at a customer and is unable to provide service for it. Thus, the lower value of DPI results in shorter planned routes and longer additional distances. In contrast, a higher value of DPI shows that the dispatcher desires to serve every customer along the planned route. Hence, a higher value of the DPI results in higher planned routes and lower additional distances. In addition, the higher value of DPI requires more vehicles, and it could lose some of customers because there are not enough vehicles to serve all the customers. Therefore, the cost of lost opportunity due to lack of service to customers is also considered in this paper. As a result, the sensitive parameters DPI and API significantly influence the sum lengths of the planned routes and additional distance and should be determined properly [24]. In this work, the additional distance caused by service failures is estimated using a stochastic simulation.
The CLRP-FD can be modeled as a directed complete graph G(V , E), with a vertex set V (Set of all candidate depots and customers) and an edge set E = { i, j |i, j ∈ V and i = j}. Before presenting the mathematical formulation of the problem, the following notations are defined.

Set and parameters:
J : Set of customers. I : Set of candidate depots. V : Set of all candidate depots and customers, demand of customer i. Q: Capacity of each vehicle. It is assumed that all of the vehicles are homogeneous.
Cost of the lost opportunity due to the lack of providing service to customer j.
Decision Variables: 1, if vehicle k goes directly from node i to node j 0, otherwise w jk = The sequence number of customer j in route k, which are auxiliary variables for eliminating sub-tour in route k [21].
AD k = Additional distance traveled by vehicle k.
Given the above defined assumptions, parameters, and decision variables, the corresponding model of the CLRP-FD based upon the credibility theory is formulated as follows: Subject to: i∈V k∈K i∈I j∈J The objective function (8) is to minimize the total costs which consists of sum of the fixed costs of opening depots, the fixed costs of using vehicles, the travel costs, the costs of lost opportunities, respectively [21]. The objective function (9) is to minimize the total additional distances. AD k is only related to the designed route of vehicle k, and the objective function (9) is calculated by the simulation procedure presented in subsection IV-D. Constraints (10) and (11) assure that all of the customers are visited within vehicle capacity and are allocated within depot capacity with a confidence level, respectively [21]. Constraints (12) ensure that each customer is served by one vehicle at most and that each customer should have one predecessor at most. Constraints (13) guarantee that each vehicle can be used only once at most, where N is the number of customers served by vehicle k. Constraints (14) state the elimination of the sub-tour for each vehicle route. Constraints (15) ensure that every vehicle will exit the customer after visiting it, and that each vehicle starts from and ends at a single depot. Constraints (16) represent that a customer must be allocated to a depot if they are connected by a route. Constraints (17) ensure the limitation on the driving distance for each employed vehicle. Constraints (18) illustrate that each customer must be allocated to exact one depot. Constraints (19) and (20) show two decision variables' relation. Constraints (21)-(23) and (24) specify the binary variables and auxiliary variables, respectively. Constraints (25) specify the continuous variables.

IV. THE HPSO FOR THE CLRP-FD
Based on the PSO and stochastic simulation, the hybrid heuristic algorithm for the CLRP-FD is proposed in this section. In a PSO-based algorithm, initially, a group of particles are initialized and every particle represents a possible solution. In the solution space, every particle moves with a given velocity from the current position to the next possible position. A successful designation of a PSO-based algorithm for the CLRP-FD depends mainly on three key issues, which are (1) creating the initial set of particles; (2) designing a desirable mapping between the particles in the PSO and the CLRP-FD solutions; and (3) moving every particle with a given velocity. In this study, a greedy cluster algorithm is used to design an efficient initialization approach and enhance the efficiency of the proposed algorithm. In addition, a pathrelinking algorithm and variable neighborhood search (VNS) algorithm are implemented to update the positions of the particles, which not only effectively avoids transforming the solutions between continuous space and discrete space [32], but also strengthens the exploitation on the promising search space.
The basic steps of the HPSO are addressed in the following subsections. Its pseudocode is depicted in Algorithm 1. It takes as parameters the maximum number of iterations, G, the number of the particles, H , and the range of the velocity, [V min , V max ]. After initializing the population in the first step (line 1), the fitness of each individual that is included in the initial population is evaluated (line 2), and the best individual in the current population is also found (line 3). The best solution found for the CLRP-FD, gbest, and the velocity of each particle, v i , are initialized in line 4 and lines 5-7, respectively. The iterations of the HPSO are executed in lines 8-27. Each HPSO iteration includes the calculation of the average velocity of each particle (line 10), the updating of the thresholds (R 1 and R 2 ) of the average velocity of particle p i (line 10), a position updating strategy (lines [11][12][13][14][15][16][17][18][19] based on the mixed path-relinking algorithm (see subsection IV-C1 for the details) and the variable neighborhood search (see subsection IV-C2 for the details), resetting the best individual (lines 20-21) and saving the best so far solution (lines [22][23]. During the optimization processes, the mixed path-relinking algorithm and the variable neighborhood search strengthen the exploitation in the neighborhood of each particle.

A. SOLUTION REPRESENTATION
A solution representation should determine the assigned customers to each vehicle, the depots to be established and the sequence of customers to be served by a specific vehicle starting from and ending at a depot [26]. Given m candidate depots and n customers, a solution representation is presented using all of the elements within the sets I = {1, 2, · · · , m} and J = {m + 1, m + 2, · · · , m + n}, where the m depots and n customers are numbered with the elements in sets I and J , respectively. The first number in a solution is always in the set I , which indicates the first depot under consideration [17]. In the solution representation, the customers assigned to each depot are between it and the next different depot. If there are no customers after a depot, it indicates that this candidate depot is not open. Each route of a given depot starts by serving the first customer, and the other customers allocated to that depot are assigned to the current route from left to right [17], [41]. If the credibility of having sufficient capacity for serving the next customer falls below the dispatcher preference index [17], the current route is terminated, and the next number in the solution representation will be a depot if there are unsigned customers. If the next number in the solution representation is still this depot, a new route will start from it to service the remaining customers allocated to it; otherwise, a new route will start from a new different depot. It can be verified that this solution representation always gives a CLRP-FD solution without violating the capacity constraints of the vehicles and depots [17], [41].
To further clarify the above description, a feasible solution of a CLRP-FD instance with m = 5 and n = 20 is shown in Figure 2. It can be seen that the candidate depot, depot 1, is not open, and the 20 customers are assigned to 4 open depots (2, 3, 4 and 5). The delivery routes related to the open depots are shown in Figure 3.

B. INITIAL SOLUTION GENERATION
Initialization of a solution plays a significant role for the performance of a meta-heuristic because the quality of the  final solution produced by meta-heuristic algorithms usually relies on the quality of the initial solution [17]. In this work, a greedy clustering algorithm is designed to generate a relatively good-quality initial solution. This greedy clustering algorithm has 4 main steps, which are clustering the customers, calculating the gravity center of each cluster, assigning the customers to the depots and constructing the routes. In the first step, a greedy search algorithm is conducted to cluster the customers. In the second step, the gravity center of each cluster is calculated for selecting the depots. Based on the consideration of the distance between the gravity centers of the clusters and the depots as well as the capacities of the depots, the clusters generated in the first step are assigned to the open depots in the third step. In the fourth step, routes such as in the traveling salesman problem (TSP) are constructed. The details of each step are presented in the following.
Step 1 (Clustering the Customers): Based on the consideration of the fuzzy demands and intra-distances of the customers and the capacities of the vehicles, the customers are clustered using a greedy search method in this step. First, to generate a cluster, a customer from the non-clustered customers is randomly selected. Then, the nearest customer to the last selected customer of the current cluster is selected. Considering the DPI value and the credibility of the customer, if the closest customer has not been clustered and its fuzzy demand satisfies the remaining capacity of the vehicle, it will be assigned to the cluster. Otherwise, it will be excluded from the cluster, and the greedy search method will search for another different customer, which is closest to the last added member of the cluster among the non-clustered customers. If the remaining capacity of the vehicle cannot serve a new customer, a new cluster will be launched, and the non-clustered customers will be assigned to the new cluster until all of the customers are assigned. Obviously, the capacity of the vehicle can be fully utilized using this procedure.
Step 2 (Calculation of the Gravity Center of Each Cluster): To establish the depots among the candidate sites, the gravity center of each cluster is firstly calculated using equation (26), where n c is equal to the number of customers assigned to cluster C, and (X c , Y c ) and (x j , y j ) represent the coordinates of the gravity center of cluster C and customer j, respectively. Second, the sum of the Euclidean distances between n customers and each candidate depot is calculated using equation (27). In (27), U i is the sum of the distances from depot i to all of the customers, and (x i , y i ) and (x j , y j ) represent the coordinates of depot i and customer j, respectively. Since most of the state-of-the-art heuristic algorithms for solving CLRP instances in the existing literature used Euclidean distance, this study also employed Euclidean distance to conveniently compare the computational results obtained by HPSO and other algorithms. If the actual traffic distance is used, it will be closer to the realistic situation, but the type of distance and calculation method cannot affect the optimization performance of the algorithm.
Third, the candidate depots are sorted in descending order according to the amount of factor W i , which is calculated by equation (28).
Step 3 (Allocating the Customers to the Depots): To assign the clusters to the ranked depots, the Euclidean distance between the gravity center of each cluster and the top-ranked depot is calculated first in this step. Afterward, the unassigned clusters are ranked in ascending order according to the distance between their gravity centers and the top-ranked depot. Algorithm 1 Pseudocode of the Hybrid Particle Swarm Optimization (HPSO) Input: G = number of iterations, H = size of the population, V max = maximum velocity of each particle, V min = minimum velocity of each particle, Output: gbest = Best solution found for the CLRP-FD. 1: Initialize a population H by the proposed greedy cluster algorithm; 2: Evaluate each individual p i ∈ H (i = 1, 2, · · · , H ); 3: pbest := the best particle in the current population; 4: gbest := pbest; 5: for i := 1 to H do 6: v i := rand(V min , V max ); 7: end for 8: for g := 1 to G do 9: for i := 1 to H do 10: Calculate Average_Vi, R 1 and R 2 . Average_Vi is the average velocity of particle p i . R 1 and R 2 are the thresholds of average velocity of particle p i . 11: if R 1 ≤ Average_Vi & Average_Vi < R 2 then 12: p i := Path_Relinking(p i , pbest); 13: else 14: if Average_Vi ≥ R 2 then 15: p i := Path_Relinking(p i , gbest); 16: else 17: p i := BVNS(p i ); 18: end if 19: end if 20: if f (p i ) < f (pbest) then 21: pbest := p i ; 22: if f (pbest) < f (gbest) then 23: gbest := p i ; 24: end if 25: end if 26: end for 27: end for return gbest The top-ranked cluster is assigned to the top-ranked depot if Cr ≥ API . The clusters in the sorted list should be allocated one by one to the top-ranked depot considering the above relation. The process of allocating clusters to a depot will be terminated when the depot does not have enough capacity to serve a new cluster. In this situation, the allocating procedure will be repeated for the next-ranked depot until all of the clusters are allocated [21].
Step 4 (Constructing the Routes): In this step, the routing problem is solved using the greedy search method for each cluster and its relevant depot. First, the nearest customer to each depot selected in step 3 is routed. Then, a process similar to that employed in step 1 is used to select the closest customer to the last added customer of the current route. If the remaining capacity of the current vehicle is not sufficient to serve a new customer, a new route will be added, and the non-routed customers will be assigned to the new route until all of the customer allocated to the given depots are routed.

C. POSITION UPDATING
In the PSO, the position of each individual particle is adjusted according to its velocity, which is updated iteratively through its personal best position (i.e., the best position found so far by the particle) and the best position found by the particles in its neighborhoods. As designed before, the position of each individual particle in this study is represented by a d−dimensional vector in problem space x i = (x i1 , x i2 , · · · , x id ), i = 1, 2, · · · , H (H is the population size, and d is the dimension of the vector) [32]. In a classic PSO, the new position and velocity of each particle can be calculated by Equations (29) and (30) [42], respectively. where where t is the loop counter, c 1 and c 2 are the acceleration coefficients, and rand 1 and rand 2 are random values generated in the interval (0,1). If Equations (29) and (30) are employed in this study, the solutions of the proposed model (8)- (25) should be transformed to continuous space. To avoid this transformation, Equation (29) is not used at all. As mentioned before, in this study, the position of each particle is updated using a Path Relinking Algorithm [19], which is an intensification strategy and is usually used to explore trajectories between solutions. Therefore, the Path Relinking algorithm follows the idea of the movement of an individual particle in a PSO; that is, when the path relinking algorithm is used to adjust the position of a particle, the particle can either move along its own path or move toward its previous optimal solution, or move toward the global optimal solution, where the roles of the starting solution and target solution are played by the particle (current solution) and the global best particle of the swarm or the current best solution of the particle. In this study, the average velocity of each particle is used to decide which target solution it follows. The average velocity of each particle, Average_V i (i = 1, 2, · · · , H ), is calculated using Equation (32). If Average_V i < R 1 (R 1 is a real number), then path relinking will not perform, and the particle will follow its own path, which is implemented using a variable neighborhood search algorithm (which will be described later); otherwise, the particle will perform a path relinking with its current best solution if R 1 ≤ Average_V i < R 2 (R 2 is also a real number) or with the global best particle if Average_V i ≥ R 2 . The two real numbers R 1 and R 2 are calculated using the following Equations (33) and (34), where G is the maximum iterations, t is the current iteration, and V max and V min are the upper and lower bounds for the velocity of each particle. If the value of the velocity of a particle is not between these bounds in some iterations, a value between these bounds will be randomly selected to initialize it. With an increase in the number of iterations, we want to increase the possibility of performing path relinking and decrease the possibility of performing the variable neighborhood search, and therefore, R 1 and R 2 have large values at the beginning of the algorithm, and the values are reduced during the iterations.

1) PATH RELINKING ALGORITHM
Path-relinking is a search intensification strategy proposed by Glover [43], [44]. As a major enhancement to heuristic search methods for solving combinatorial optimization problems, its hybridization with other meta-heuristics has led to significant improvements in both the solution quality and running time [43]. Because the ideas behind the path-relinking algorithm and the movement of an individual particle are very similar [32], the path-relinking algorithm is employed in this paper to improve the current solution and make it close to the target solution. In the path relinking strategy, the current solution plays the role of the starting solution, and the best particle of the swarm or the current best solution of the particle plays the role of the target solution [32]. If the Hamming distance (the Hamming distance between two solutions is the number of positions at which the corresponding elements are different), HD, between the starting solution and the target solution is not equal to 0, then the trajectories between the two solutions are explored by a cycling exchange operator until the value of HD is decreased to 0. Algorithm 2 shows the pseudocode of the path-relinking algorithm employed in this study. First, the best solution produced by this optimization procedure is initialized in line 2, and the Hamming distance between the starting solution x and the target solution x t is calculated in line 3. Then, the trajectories between the two solutions are explored by a cycling exchange operator (lines 5-31) if HD > 0. If the variable Temp is equal to 1, the exchanging procedure is implemented from the left-hand side of the solutions (lines 6-14); otherwise, it is implemented from the right-hand side of the solutions (lines [15][16][17][18][19][20][21][22][23][24]. The solution x obtained by the exchange operator replaces the best solution x * if its fitness is better than the one of x * (lines [25][26][27]. The starting solution x is also updated using the solution x in line 28. To find solutions with better quality, the roles of the starting and target solutions are interchanged (line 29) at each iteration of the while HD > 0 do 6: if Temp = 1 then //Comparison between x and x t starts from their left hand side 7: jst := the first left position in which the elements located in x and x t are different; // Comparing the elements in x and x t , the first position from the left hand side of x and x t is returned by jst, where x and x t have different elements 8: Temp = 0; 9: for i := 1 to jst do 10: x i := x t i ; 11: end for 12: for i := jst + 1 to n do 13: x i := x i ; 14: end for 15: else//Comparison between x and x t starts from their right hand side 16: jst := the first right position in which the elements located in x and x t are different; // Comparing the elements in x and x t , the first position from the right hand side of x and x t is returned by jst, where x and x t have different elements 17: Temp = 1; 18: for i := n to jst do 19: x i := x t i ; 20: end for 21: for i := jst − 1 to 1 do 22: x i := x i ; 23: end for 24: end if 25: if f (x ) < f (x * ) then 26: x * := x ; 27: end if 28: x := x ; // Update the starting solution x using x 29: x ←→ x t ; // Interchange the roles of the starting and target solution 30: HD := HD − 1; 31: end while 32: return x * ; 33: end function path-relinking algorithm. The Hamming distance is updated in line 30, and the procedure repeats until HD = 0. Finally, the best solution obtained by this optimization procedure is returned in line 32.

2) VARIABLE NEIGHBORHOOD SEARCH
As mentioned previously, if Average_V i < R 1 for individual particle i, the Variable Neighborhood Search (VNS) VOLUME 8, 2020 Algorithm 3 Pseudocode of the Variable Neighborhood Search Algorithm 1: function VNS(x, t max , q max ) 2: for t := 1 to t max do 3: for q := 1 to q max do 4: x ←Shake(x, q); //Generating a point x at random from the kth neighborhood of x, i.e., x ∈ N q (x) 5: x ←FirstImprovement(x , q); //Local search; Finding a point x from the qth neighborhood of x where f (x ) ≤ f (x ) 6: if f (x ) < f (x) then 7: x = x ; algorithm is applied as a local search strategy to improve its position. Given the CLRP-FD defined on a solution space X , a subset N (x) ⊆ X is called a neighborhood of solution x where x ∈ X . {N q }(q = 1, 2, · · · , q max ) is the set of neighborhood structures selected to be explored during the search procedure, and N q (x) is the set of solutions in the qth neighborhood of x. In this study, the VNS algorithm shown in Algorithm 3 is used in the following way. The VNS starts with an initial solution x and its neighborhoods set {N q }. At the tth iteration, a solution x is randomly generated from its qth neighborhood (line 4). Then, as a local search procedure, the first improvement is applied to x , and a solution x from the qth neighborhood of x is obtained (line 5). If x is better than x in terms of fitness, the current solution x is replaced by its better neighbor x (line 7), and the neighborhood counter is initialized to 1 (line 8). The process consecutively investigates all of the neighborhoods. Otherwise, the neighborhood counter is increased by one (line 10), which indicates that a higher-order neighborhood will be explored next [43].
The variable neighborhood search algorithm shown in Algorithm 3 combines stochastic and deterministic changes in the neighborhoods. The random selection of a point x from the qth neighborhood (line 4) represents the stochastic phase. A local search based on the first improvement heuristic represents the deterministic phase. The first improvement heuristic starts at an initial solution x, and then, it makes a move as soon as a direction of decent within the neighborhood N q (x) is found. This procedure is summarized in Algorithm 4.
In this paper, the VNS algorithm employs 6 neighborhood structures, N 1 , N 2 , · · · , N 6 . To describe them, the sequence codes shown in Figure 2 are used as an example. In total, 6 different operators are implemented on it to obtain 6 corresponding neighborhood structures, which are presented in Figures 4-8. i := 0; x := x; 3: while i < N q (x) and f (x) ≥ f (x ) do 4: i ← i + 1; 5: x ← arg min{f (x), f (x i )}, x i ∈ N q (x) 6: end while 7: return x; 8: end function Neighborhood N 1 . We randomly select two elements, which are located in position i and j (j > i), e.g., 13 located in position 9, and 4 located in position 24. The element located in position i is moved to position j, and the other elements from i+1 to j are separately moved to the position previous to them. The new sequence obtained using this shifting operator is shown in Figure 4.
Neighborhoods N 2 and N 3 . The approach of generating neighborhoods N 2 and N 3 is similar to the approach of generating N 1 . The difference is that the sequential elements from position i to i+l are selected and moved to the positions from j − l to j correspondingly (j > i and j / ∈ {i, i + 1, · · · , i + l}), and the other elements from i + l + 1 to j are separately moved forward l + 1 positions. According to the numerical experiments, l is set to 2 and 3 in this study. Here, l = 2 and l = 3 correspond to neighborhoods N 2 and N 3 , respectively. For example, if i = 9, j = 24 and l = 2, then the sequential elements 13, 24 and 23 are moved to the positions j − 2, j − 1 and j, respectively. The result for neighborhood N 2 is shown in Figure 5.
Neighborhoods N 4 and N 5 . Neighborhoods N 4 and N 5 are obtained by randomly exchanging one pair of customers and depots, respectively. The swapping between a customer and a depot is not allowed in N 4 and N 5 . For example, new sequences are produced by exchanging customers 13 and 25 in Figure 6, and depots 2 and 5 in Figure 7.
Neighborhood N 6 . Neighborhood N 6 is obtained using the inverse operator. We randomly select two points, e.g. 11 and 14. A new sequence shown in Figure 8 is obtained by reversing the subsequence between these two points [3].
Note that the solution obtained by conducting the local search in these above neighborhoods can be infeasible. In this case, a new solution must be encoded again to ensure its feasibility [45]. If the first position is selected, then the other position must be a depot to maintain feasibility [45].

D. STOCHASTIC SIMULATION
As stated previously, each customer's demand is assumed to be a triangular fuzzy number, and it cannot be directly considered as a deterministic number [21]. To reveal the actual demand of each customer, a stochastic simulation is employed [21], [24]. For each feasible planned route, additional distances due to route failures are also provided using this simulation experiment. Its main steps are summarized as follows:     Step1: Generating the actual demand of each customer by using the following processes: (1) randomly generate a real number r in the interval between the left and right bounds of the triangular fuzzy number representing demand of the customer, and compute its membership µ [21]; (2) randomly generate a real number α in interval [0, 1], if α ≤ µ, the actual demand of the customer is adopted as r.
Otherwise, randomly generate r and µ again until the relation α ≤ µ is fulfilled; (3) repeat (1) and (2) until every customers have an actual simulation demand [21], [24]. Step2: Move along each distribution route generated by the proposed HPSO and calculate the additional distances because of route failure in terms of actual demands. Step3: Repeat the above two steps M times (In this study, M is equal to 400 which is a proper value suggested in literature [24]), and compute the average additional distances which are considered as the additional distances.

V. NUMERICAL EXPERIMENTS
In order to study the performance of the proposed HPSO for solving the CLRP-FD, a number of computational experiments were carried out. Matlab 2017a was used for the implementation of the algorithm. All the experiments were carried out on a personal laptop equipped with Inter core i5 4200M CPU at 2.50GHz and 8GB of RAM under Windows 10 professional system.

A. SENSITIVITY ANALYSIS AND PARAMETER SETTINGS IN HPSO
The quality of the solution generated using HPSO can be affected by the parameters, namely, the number of iterations, G, the population size, H , acceleration coefficients, c 1 and c 2 , and w 1 and w 2 , which are used to control the searching method. For simplicity, the parameters in HPSO were fixed in the initialization stage and did not change during the running of the algorithm. This method is usually called parameter tuning in the literature. Because we would like to make the computational results as reproducible as possible, the runs in all of the experiments were limited to a single set of parameter settings generated using parameter tuning.
To illustrate how the parameters affect the overall performance of the HPSO, an analysis between the approximate optimal solution and the crucial parameters (H , c 1 , c 2 , w 1 and w 2 ) was conducted on 5 benchmark CLRP instances (Gaskell67-21 × 5, Gaskell67-32 × 5-2, Gaskell67-36 × 5, 20-5-1a, 50-5-1a). For the detailed reasons why we used CLRP instances instead of CLRP-FD instances, please refer to Subsection V-C2. When investigating the relationship between one (or one pair) of parameters and the approximate optimal solution, the other parameters were kept unchanged. Each instance was independently tested 10 runs.

1) NUMBER OF ITERATIONS
The number of iterations plays an important role in HPSO. A large number of iterations can waste computing time, while a small number of iterations can stop the search process early and cannot obtain an optimal/near optimal solution. Table 1 presents the results of each tested instance over 10 runs, including the names of the instances, number of candidate depots (m), number of customers (n), optimal or best-known objective value (BKOV) reported in the literature, average of the best objective values over all of the runs (BOV), how many runs with which the best-known solution was found (Found), average percentage deviation from the best-known solution, APD = BOV−BKOV BKOV × 100%, average CPU time in seconds (Time), and average iterations (Aiters) in which the best-known solution was found.
From the results shown in Table 1, we can observe that all 5 tested instances obtained their optimal or near optimal solutions within 1500 iterations. Because large instances can require a large number of iterations, we set G = 2000 in this study. To obtain a high-quality solution without wasting a large amount of computing time, a stopping criterion was also used in the HPSO. If the best solution found thus far could not be improved after G/4 iterations, the search process of the HPSO was terminated.

2) POPULATION SIZE
For setting the parameter H , we considered 6 different population sizes in the range [30,50] with a step of 5. Table 2 presents the average results of the 5 tested instances over 10 runs, including the population size H , the average number of runs that produced the best-known objective values, the average CPU time in seconds, and the average value of APD. Table 2 shows that the average number of runs and the average values of APD are relatively better when H is equal to or greater than 40. In addition, the larger the value of H is, the smaller the value of APD, but the more CPU time the HPSO needs. Considering both the solution quality and computing time, we set H = 40 in this study.

3) PARAMETERS w 1 AND w 2
The parameters w 1 and w 2 can not only decide which target solution a particle follows but also determine whether the VNS algorithm or path relinking is performed. According to Formulas (33) and (34), we can determine that the greater the values of w 1 and w 2 are, the more the VNS is performed during the first several iterations of the HPSO. To set w 1 and w 2 , we considered 5 different combinations for their values, following the literature [32]. Table 3 shows that when w 1 = 0.8 and w 2 = 0.9, the average value of APD and the average number of runs that produced the best-known objective values are comparatively better, and they are 0.38% and 5, respectively. Hence, we set w 1 = 0.8 and w 2 = 0.9 in this study. Average results of the 5 tested CLRP instances over 10 runs obtained using different population size.

TABLE 3.
Average results of the 5 tested CLRP instances over 10 runs obtained using different combinations of w 1 and w 2 .

TABLE 4.
Average results of the 5 tested CLRP instances over 10 runs obtained using different combinations of c 1 and c 2 .

4) ACCELERATION COEFFICIENTS c 1 AND c 2
As shown in Formulas (30) and (31), c 1 and c 2 are used to compute the velocity of each particle, which has a large influence on the performance of the HPSO. Table 4 presents the average results of the 5 instances over 10 runs that were obtained from the HPSO using different combinations of c 1 and c 2 , where both c 1 and c 2 were greater than 2, which followed the literature [32].
From the results presented in Table 4, we can determine that the first combination is relatively competitive, with which HPSO could produce the best-known solutions in 5 out of 10 runs, and it obtained a lower average percentage deviation within lower computational time.
According to the above discussion and experimental experience, throughout the experiment, the following parameters were used: the maximal number of iterations G = 2000, the number of particles H = 40, the acceleration coefficients c 1 = c 2 = 2.05, and parameters w 1 = 0.8, and w 2 = 0.9.

B. SENSITIVITY ANALYSIS ON PARAMETER DPI
In these experiments, a set of instances is generated using the approach presented in [21]. A small size instance with 30 customers and 5 candidate depots, a medium size instance with 50 customers and 5 candidate depots, and a large size instance with 100 customers and 7 candidate depots were considered. The coordinates of all of the customers and depots were generated randomly in the interval [100 × 100] [21]. In addition, the fuzzy demands of the customers, which are triangular fuzzy numbers such asd = (d 1 , d 2 , d 3 ), were generated randomly [21]. In this way, d 1 , d 2 and d 3 were generated within the intervals [10,35], [36,60] and [61, 110], respectively [21]. The data for the three test instances are shown in Table 5.
The values of DPI and API were also set following the approach presented in the literature [21]. The values of DPI varied from 0.1 to 1 with a step of 0.1. For the sake of simplicity, the value of API was always considered to be 1, which reduces the number of investigation cases. Tables 6-8 give the average computational results of 10 runs for the three test instances, respectively. The columns included in these three tables are the dispatcher preference index (DPI), the planned routes, the additional distances, the routing costs, the lost opportunity costs, the depot costs, the vehicle costs, the total costs (which consists of the routing costs, depot costs, vehicles costs and lost opportunity), and the CPU time obtained by the proposed algorithm. The results of Tables 6-8 are also depicted in Figures 9-11. As shown in Tables 6-8 and Figures 9-11, the total cost reaches its minimum when DPI=0.6. This value is the same as the value presented in the literature [17], [21], [24]. Figures 9-11 show the tendencies of the planned routes, additional distances, routing costs and total costs, which vary with DPI from 0.1 to 1. As shown in Figures 9-11, lower values of DPI have shorter planned routes and longer additional distances. This finding indicates that lower values of   DPI have a tendency to use all of the vehicle capacity, which results in increased total additional distances because the vehicles visit the customers but do not have enough capacity to serve them. A higher value of DPI denotes a tendency to use less vehicle capacity, which results in longer planned routes and shorter additional distances due to ''route failure''. Moreover, Tables 6-8 also show that the cost of employing vehicles is increased with increases in the value of DPI. To ensure higher service to customers at a higher value of the DPI, fewer customers are considered in each vehicle route, and thereby, more vehicles are employed. Although the cost of lost opportunity should also be increased with the increase in the DPI because of the unavailability of vehicles, its values shown in Tables 6-8 are always equal to ''0'' because enough vehicles can be employed by each candidate depot in this study. Considering the total cost, the proper value of DPI should be approximately 0.6.

C. PERFORMANCE OF THE PROPOSED HPSO 1) EFFECTIVENESS OF THE POSITION UPDATING STRATEGY IN THE HPSO
The main difference between the traditional PSO and HPSO is the position updating strategy employed in this paper, which consists of a path relinking algorithm and VNS. To study the effectiveness of the proposed position updating mechanism, we compared 2 PSOs with the HPSO on six benchmark CLRP instances. PSO1: the proposed HPSO with the path relinking algorithm only to update the position. This is to evaluate whether the VNS employed in the HPSO enhances the exploration to improve the overall performance of the HPSO. PSO2: the proposed HPSO with the VNS only to update the position. This is to evaluate the effectiveness of the path relinking algorithm.
The mean and standard deviation (SD) values of the best objective values produced by PSO1, PSO2 and HPSO over 10 runs are presented in Table 9. It can clearly be seen that the HPSO outperforms PSO1 and PSO2 in every instance with regard to the mean values. PSO1 provided the worst mean values, and the difference in HPSO versus PSO1, and PSO2 versus PSO1 is quite significant. This finding indicates that the VNS can effectively improve the performances of HPSO and PSO2. Although the VNS is mainly for local exploitation, 5 different neighborhoods were allowed to explore multiple promising areas in the search space, and this strategy also enhanced the exploration ability. Therefore, HPSO and PSO2 achieve better performance than PSO1. In addition, compared with PSO2, the position    updating mechanism is more effective in guiding the search toward promising areas in the search space because it takes advantage of the best particle in the current population and the global best particle during the evolution.
According to the results shown in Table 9, we can conclude that HPSO is more stable. To support this observation, we present the box plots of the three PSOs on the six selected instances in Figure 12. In Figure 12, the box plots of the HPSO and PSO2 are enlarged for clearly observing the variability of the results. It is no doubt that HPSO has the best stability.

2) COMPARISON WITH SOFTWARE CPLEX
Because of the fuzzy demands of the customers and the sites of all customers and depots are randomly generated in this study, it is not possible to compare the HPSO results shown in Tables 6-8 with previous work. To evaluate the efficiency of the HPSO, the solutions shown in Tables 6-8 are compared with their lower bounds. In this study, the assumption of each customer's demand is relaxed for obtaining the lower bound, and each customer's demand equals its left boundary d 1 . Thus, the constraints (10) and (11) in the model of CLRP-FD are replaced by the constraints (35) and (36), Apparently, compared with the case of fuzzy demands, the total demand of the customers is decreased when the left boundary of the fuzzy demands is considered to be the deterministic demands of the customers [21]. Moreover, the utilization of depots and vehicles should be decreased because of the reduced demands of all of the customers. Thus, the total cost in the case of deterministic demands provides a lower bound to the cost in the case of fuzzy demands. Table 10 shows the results of three test instances, with their lower bounds. The first column gives the sizes of the test instances, and columns 2-4 show the quality of the results, including the total cost obtained using the HPSO, lower bounds obtained by the commercial solver CPLEX12.8, and the gap between them, where Gap = the total cost obtained using HPSO−Lower bound Lower bound ×100%. The last two columns show the computing time of the HPSO and CPLEX solution.
As shown in Table 10, the lower bounds provided by CPLEX for the three test instances are 704.59, 821.4 and 1008.56, respectively, which have not been proved optimally within 3 hours. In other words, CPLEX did not optimally solve these three test instances within 3 hours. Therefore, the proposed HPSO is more efficient than the CPLEX solver in terms of the computing time, which can provide a satisfactory solution within a shorter CPU time.  Table 11). To use the standard CLRP test instances in the case of CLRP-FD, the following steps that were proposed in [21] are implemented on the CLRP instances to change them into CLRP-FD instances: Step 1: Based on the deterministic demand of each customer, we generate the corresponding triangular fuzzy demand as follows: (1)   Clearly, the optimal solution of the CLRP instance will be a lower bound for the optimal solution of its corresponding CLRP-FD instance, obtained using the above steps. Table 11 summarizes the computational results of 8 standard CLRP instances and their corresponding CLRP-FD instances. The first column of Table 11 gives the names of the test instances. The second column shows the solutions of the standard CLRP instances obtained using the CPLEX solver within 3 hours. The next 5 columns present the solutions of the corresponding CLRP-FD instances obtained using HPSO with DPI=0.6 and API=1. The Gap between the HPSO solution and its lower bound is presented in the last column. The values of gap indicate that the total cost of the CLRP-FD instances is at least 28.97% larger than that of its corresponding CLRP instances. Therefore, the decision makers should cautiously make the decisions in this situation, which lacks sufficient information to make a precise judgment.

3) COMPARISON WITH THE STATE-OF-THE-ART HEURISTIC ALGORITHMS
To thoroughly evaluate the performance of HPSO, we also compared it with seven state-of-the-art algorithms presented in the literature, which are listed as follows.
• Greedy randomized adaptive search procedure (GRASP) [46] • Memetic algorithm with population management (MA | PM) [47] • Lagrangean relaxation and granular tabu search method (LRGTS) [48] • Two-phase hybrid heuristic algorithm (2-Phase HGTS) [49] • Multiple ant colony optimization algorithm (MACO) [50] • Hybrid genetic algorithm (HGA) [51] • Granular variable tabu neighborhood search (GVTNS) [52] Table 12 provides a comparison of 27 benchmark CLRP instances between HPSO and the above popular heuristic algorithms. The first 11 instances are from the Barreto dataset, and the remaining 16 instances are from the Prins dataset. Table 12 presents the name of the instance, the number of candidate depots (m), the number of customers (n), the best known objective value (BKOV), the best objective value (BOV) obtained by each algorithm, and the average computational time (Time) of each algorithm for solving each instance. In Table 12, the figures shown in bold indicate the best results among the compared algorithms. From Table 12, it can be seen that GRASP, MA|PM, LRGTS, 2-Phase HGTS, MACO, HGA, GVTNS, and HPSO obtained optimal results in 7, 10,5,11,17,20,14, and 18 compared instances, respectively. HPSO did not produce optimal solutions for 9 of the 27 instances compared to the best-known results presented in the literature. It is evident that except for HGA, HPSO produced better solutions than the other six compared algorithms. Therefore, we can conclude that HPSO has better performance than (or similar performance with) the other six popular algorithms in terms of the solution quality. We should also note that the state-of-the-art hybrid genetic algorithm presented in [51] incorporates additional heuristic components. Apart from the selection and crossover operators, the mutation operator in the HGA is very complicated, which employs the add operator, swap operator and multiple local searches. One of the crucial things in the optimization process of HGA is to determine the proper crossover and mutation rate. Compared with the HGA, the proposed HPSO has a simpler framework from the point of implementation.
In Table 12, it can be seen that HPSO spent much more computational time in solving all of the compared instances. We should note that the different hardware, compilers, and programming languages have a large impact on the heuristic algorithms' effectiveness. Table 13 presents the implemented platforms and programming languages of the comparison algorithms. It is possible that the slower MATLAB coder and larger number of generations and local search iterations resulted in a larger computational time for the HPSO.
We also conducted statistical tests using the statistical software SPSS to determine whether the above observations between the HPSO and the compared algorithms are statistically significant. The BOVs of each instance presented in Table 12 were considered for the statistical tests.
The results of the nonparametric Friedman test are shown in Table 14, where the significance level α = 0.05. Table 14 presents the average ranks calculated for the compared algorithms, the Friedman χ 2 F statistic values, p-values, and the decision about acceptance or rejection. As shown in    with 7 degrees of freedom was 14.07), and the null hypothesis was rejected at the significance level of α = 0.05, which indicates that a statistically significant difference exists among the compared algorithms. The Wilcoxon signed-ranks nonparametric test was also used to investigate the differences between the two compared algorithms. The results of the Wilcoxon tests are presented in Table 15, which include the p-values and the acceptance or rejection of the null hypothesis for each pair of the compared algorithms.
As the results in Table 15 show, HPSO did not always confirm statistically significant differences between the HPSO and the seven compared algorithms. Statistically significant differences were observed between HPSO and GRASP, MA|PM, and LRGTS. This finding indicates that HPSO outperformed GRASP, MA|PM, and LRGTS in solving the instances presented in Table 12. The null hypotheses were accepted in comparisons of HPSO with 2-phase HGTS, MACO, GVTNS, and HGA. This finding indicates that the performance of the HPSO in solving the instances in Table 12 is similar with 2-phase HGTS, MACO, GVTNS and HGA.
The above observation illustrates that HPSO can perform well on most of the CLRP instances and is competitive with the existing state-of-the-art algorithms.

VI. CONCLUSION
In this paper, the CLRP-FD was studied. Based on credibility theory, a fuzzy chance constrained programming model was formulated for this problem. For the purpose of solving this model, a hybrid heuristic algorithm based on the PSO was proposed. In the HPSO, the path relinking algorithm is used to update the position of each particle, and local search is implemented using VNS to improve the quality of the solution. In the solution process of the HPSO, stochastic simulation was employed to estimate the additional distances that resulted from fuzzy demands and route failures for every planned route. To verify the validity of the presented model and evaluate the performance of the proposed HPSO algorithm, three numerical instances with different sizes, which were generated according to the approach shown in the literature, were solved. The experimental results showed that the total distribution cost can be greatly influenced by the value of the dispatcher preference index (DPI). In addition, to investigate the efficiency of the proposed HPSO, the lower bounds of the solutions were computed using the optimization solver, CPLEX 12.8. Results revealed that CPLEX could not optimally solve the tested instances within 3 hours, but the proposed HPSO could obtain satisfactory solutions in a reasonable computation time. To further illustrate the performance of the proposed model, 8 CLRP-FD instances were generated based on the corresponding standard CLRP instances with some changes. The Gap between the HPSO solution with DPI=0.6 and API=1 and its lower bound is at least 28.97%, which indicates that the decision makers should cautiously make the decisions in the absence of sufficient information. We also compared HPSO with seven state-ofthe-art algorithms presented in the literature, and the results showed that HPSO has better or similar performance with these state-of-the-art heuristic algorithms.
To extend this work, more future work can be done as follows: (1) considering other parameters of the CLRP as fuzzy variables, such as the travel time and time windows, (2) developing other variants of the CLRP-FD, such as the CLRP-FD with open routing tours and the CLRP-FD with time windows, (3) designing other meta-heuristic algorithms to solve the CLRP-FD and its variants, and (4) integrating exact technique/optimizer like CPLEX into the HPSO to solve subproblems in an optimal way. From July 1988 to May 1994, he was an Assistant Professor with the Department of Management Engineering, Shanghai Institute of Mechanical Technology. Since June 1994, he has been a Professor of management science and engineering with the University of Shanghai for Science and Technology. He is the author of eight books and more than 300 articles. His research interests include operations research, management system engineering, and meta-heuristic algorithms.
Prof. Ma's awards and honors mainly include the Shanghai Shuguang Scholar Award and the Shanghai Yucai Award from the Shanghai Municipal Education Commission.
ZIYING ZHANG was born in Hubei, China, in 1978. He received the Ph.D. degree in materials science from Fudan University, Shanghai, China, in 2013.
Since September 2013, he has been an Associate Professor with the School of Materials Engineering, Shanghai University of Engineering Science. He has participated more than ten academic research projects as a Principal Investigator/Director or a Key Participant. He is a Special Expert on Development in Northern Jiangsu Province of China. His recent research interests include the corrosion and protection of metallic materials, the theoretical research and computational simulation of the structure and properties of IOL materials, and battery electrochemistry.