Enhanced intelligent water drops algorithm for multi-depot vehicle routing problem

The intelligent water drop algorithm is a swarm-based metaheuristic algorithm, inspired by the characteristics of water drops in the river and the environmental changes resulting from the action of the flowing river. Since its appearance as an alternative stochastic optimization method, the algorithm has found applications in solving a wide range of combinatorial and functional optimization problems. This paper presents an improved intelligent water drop algorithm for solving multi-depot vehicle routing problems. A simulated annealing algorithm was introduced into the proposed algorithm as a local search metaheuristic to prevent the intelligent water drop algorithm from getting trapped into local minima and also improve its solution quality. In addition, some of the potential problematic issues associated with using simulated annealing that include high computational runtime and exponential calculation of the probability of acceptance criteria, are investigated. The exponential calculation of the probability of acceptance criteria for the simulated annealing based techniques is computationally expensive. Therefore, in order to maximize the performance of the intelligent water drop algorithm using simulated annealing, a better way of calculating the probability of acceptance criteria is considered. The performance of the proposed hybrid algorithm is evaluated by using 33 standard test problems, with the results obtained compared with the solutions offered by four well-known techniques from the subject literature. Experimental results and statistical tests show that the new method possesses outstanding performance in terms of solution quality and runtime consumed. In addition, the proposed algorithm is suitable for solving large-scale problems.


Introduction
Managing a fleet of vehicles outsourced for the distribution of a specific number of products to a set of customers with specific supply and demand requirements, is considered an important challenge in dealing with distribution problems. The challenge in this case, is not only restricted to making decisions on the number of vehicles to be dispatched on the road, but in deciding how a customer receives a service and which customer receives services first, based on assigned priority [1]. This type of problem can be presented as a vehicle routing problem (VRP) modeled by using graph theory. The expected performance metrics involve determining the optimal PLOS  sequence of customers to be visited by each vehicle, which satisfy the criteria such as travel time, the length of route and the cost involved in the operation [2]. The VRP optimization problem is widely-studied with several attractive solutions and different implementation techniques proposed in the literature [3,4,5,6]. Similarly, several variants of the VRP, of which the design concepts are based on the operational mechanism and mathematical modeling of the problem's diverse conditions in real-world applications, have been scrutinized recently by researchers [7,8,9,10,11]. Some of these variants of the VRP include capacitated vehicle routing problems (CVRP), heterogeneous fleet vehicle routing problems (HVRP), multi-depot vehicle routing problems (MDVRP), periodic vehicle routing problems (PVRP), stochastic vehicle routing problems (SVRP), and vehicle routing problems with time windows (VRPTW) [12,11]. The classical VRP is the most studied version of the VRP's constrained optimization problem, drawing substantial interest in literature. In the classical VRP, there is only one depot and all vehicles start and end their routes in that depot [13,14]. However, in present day real-world applications the classical VRP has limited practical application due to the complex structural orientation of the real-world problem. The VRP and its variants, including the multi-depot vehicle routing problem (MDVRP) discussed in this paper, are considered to be NP-hard problems, since they cannot be solved in polynomial time [15,16]. In other words, finding an optimal solution for the VRP usually is very difficult to achieve, requiring an extensive amount of computational effort.
Since VRP is an NP-hard problem, several methods are proposed for solving VRPs and its variants. These methods include the exact algorithms, heuristics and metaheuristics techniques. However, using the exact algorithms are usually not effective, due to their time consuming characteristics when applied, which become obvious when the number of nodes increases and connectivity causes a drastic increase in the computational time required. Examples of the exact solutions are column generation [17,18], integer programming [19], branch-and-bound and branch-and-cut algorithms and set-covering-based solution methods [11]. As a result of the underperformances of the exact algorithms with large graphs, a large number of classical heuristics or evolutionary computing methods and metaheuristics algorithms, has been proposed to solve the VRPs. This include simulated annealing [20], Tabu search [21], genetic algorithms [22], ant algorithms [23], neural networks [24], particle swarm optimization [25], symbiotic organisms search algorithm [26], and intelligent water drop algorithms [27,28].
In the MDVRP, more than one depot exists and each customer is visited by one of the vehicles based at one of the several depots. The MDVRP has more practical applications than the classical VRP, as it is inclined toward offering better services to multiple customers from several depots. Although the MDVRP is said to be more robust in its application, it presents a greater challenge to the decision makers in determining which customer is visited by a vehicle from which depot without exceeding the actual capacity constraints of each vehicle. The inherent complexity of the MDVRP makes it difficult to solve the problem on a larger scale. Therefore, the MDVRP can then rather be viewed as a clustering problem, in which case grouping is performed to cluster the customers based on distance between the customers and the depots. Therefore, the problem can be solved in two phases simultaneously. The first phase is the allocation of customers to depots and the second phase is the routing of customers allocated to the same depots through routes. To adequately handle large problems, a reasonable approach would be to further decompose the two phases into the following sub-problems, namely clustering, routing, scheduling and optimization. Although most of the existing classical heuristic algorithms, such as genetic algorithms [22] and simulated annealing [20] guarantee near-optimal solutions, this is not always the case when solving MDVRP with large problem sizes, where for example the number of possible solutions would increase with every slightest increase in the number of decision variables.
Consequently, the computational effort is most likely to increase likewise, thereby supporting the view that classical based heuristics are not computationally efficient in solving MDVRPs with large problem sizes. However, several hybrid evolutionary algorithms have proven to perform better, in terms of being able to produce better solution quality than their classical counterparts [15,29,30,31], especially in cases of large sized problems. Therefore, this serves as a motivation for proposing a hybrid metaheuristic algorithm that combines both the intelligent water drop (IWD) algorithm and simulated annealing (SA) local search metaheuristic, denoted here as IWD-SA to solve the MDVRP.
IWD is a graph-based metaheuristic algorithm, which makes it suitable for modeling and solving most VRPs and their variants including the MDVRP. The IWD algorithm is a very simple and effective population-based optimization technique, using a constructive approach in finding optimal solution for a given problem [32,33]. The IWD is inspired by the natural phenomena, based on the idea of water drops and their interactions with soils in river beds. The process is such that each water drop would construct a solution by traversing in the problem search space, at the same time modifying its environment. IWD has found applications in dealing with a wide range of optimization problems which include the well-known travelling salesman problem [34], VRP and its variants [29] and software quality assurance testing [35]. Results from different researchers' work, have shown that the IWD algorithm compete favorably well with the other state-of-the-art metaheuristic algorithms [27,28,29].
This paper proposes the development of a hybrid metaheuristic algorithm based on the IWD algorithm and the SA local search heuristic. Since the SA algorithm has been applied to solve a number of optimization problems with fairly good results in most cases [36], it was selected for its high and better objective values, and for the SA algorithm's ability to move from the current solution to the neighborhood solutions. This invariably helps the search process escape from the local minima in its search for the global optima by using the specified acceptance probability criteria to either accept or reject solutions with worse objective values. However, the computation of the acceptance probability function consumes a large number of system resources, more specifically CPU time. This is as a result of the exponential calculation required to determine the probability of acceptance or rejection of a new solution. Therefore, approximating the calculation of this function without compromising the decision rule and solution quality, can significantly improve the performance of the framework in terms of cost of execution. In addition, this paper aims to find an even more efficient and effective search strategy and an optimal set of performance parameters to achieve better results and faster convergence speed, especially with MDVRP benchmark problems with graphs ranging from 50 up to 360 nodes. In addition, the comparative analysis of the hybrid methods proposed in this paper is also novel. Therefore, this can be considered as the primary motivation for developing the proposed hybrid algorithm, which consists of the Graph-based IWD and SA search algorithms. The proposed hybrid algorithm is enhanced with the capability of improved IWD and SA algorithms to explore and exploit the solution search space of the MDVRP in a more efficient and effective way. The proposed algorithm is evaluated on thirty three (33) MDVRP benchmark instances. The experimental results obtained demonstrate the efficiency of this method over other compared algorithms from the existing literature [37,30,10].
The main technical contributions of this paper are the • design and implementation of an improved IWD based SA algorithm for solving MDVRP; • employment of some new adjustments features to SA probability of acceptance criteria function to improve the overall computation speed of the proposed IWD-SA algorithm; • comparative analysis of the proposed IWD-SA method presented in this paper, in terms of computational cost incurred by the exponential probability function execution of the SA technique; and • detailed evaluation and statistical analysis of results obtained by the proposed IWD-SA method against its classical approach (that is the standard IWD) and other existing techniques from literature.
This article is structured as follows: The paper starts with the brief discussion on the basics of MDVRP, problem formulation and related works. This is followed by a discussion of the two metaheuristic algorithms, namely IWD and SA, the presentation of hybrid methods, and their applications to solve MDVRP. Presented next are extensive computational results and discussions on the comparative analysis of the obtained results using statistical methods. Finally, the concluding remarks and future directions for the study are presented.

Multi-depot vehicle routing problem
The MDVRP was presented as a least cost problem, with the objective of finding routes with the least cost from each designated depot to a set of geographically located customers [36]. In modeling the MDVRP, certain assumptions regarding the vehicle routing plan, capacity, and customers are considered. Firstly, each route begins and ends at the same depot. Secondly, each customer is served exactly once a by vehicle. Thirdly, the total demand on each route is less than, or equal to the capacity of the vehicle assigned to that route. Finally, customers' demand can be met. An example illustration of the MDVRP with 2 depots and 15 customers is shown in Fig 1. In solving the MDVRP, three decision making processes are involved [38]. The first is clustering, which deals with the grouping of customers to be served by the same depot, based on the distance of each customer from the servicing depots. The second is routing, which is the assignment of customers of the same depot to several routes in such a way that the capacity constraint of the vehicle is not violated. The last is scheduling, which handles the delivery sequence of each route in every depot.
The main objective of the MDVRP is to minimize the total delivery distance or time spent in attending to each customer. As shown in Fig 1, the two rectangular boxes represent the actual depots, while the circles represent the actual customers to be visited. If The cost matrix satisfies the triangle inequality when- In particular, this is the case of planer problems for which the nodes are points is the Euclidean distance. The triangle inequality is also satisfied if d i,i+1 is the length of a shortest path from i to i + 1 on G.
For the given MDVRP, the distance between customer i and depot h, is represented as . The calculated Euclidean distance between customers and depots can then be used to make grouping decision of assigning customers to specific depots and routing decision of assigning customers on the same link to several routes. According to Ho et al. [38], a feasible solution can easily be generated based on three steps namely, clustering, routing, and scheduling. The clustering process is achieved by computing the distance matrix, since the main objective of the MDVRP is to minimize the total time spent in serving each customer. Therefore, customers are clustered according to the nearest depots to them. Considering the example given in Fig 1, with two available depots say h A and h B , each customer is assigned to exactly one depot based on the following conditions: represents the distance between customer c i and depot h k . The scheduling process entails sequencing the assignments of vehicle to customers based on the customers' proximity to each other, such that no customer is visited more than once. In the proposed method, the parallelism capability of the IWD algorithm is employed to solve the MDVRP using similar steps of clustering and routing techniques described above.

Related work
Literature reviewed shows that several metaheuristic approaches have been proposed and used to satisfactorily solve the MDVRP with efficient approximation of the optimal solution. Pisinger and Ropke [39] proposed a unified heuristic, which was able to solve five different variants of the vehicle routing problem by using the Adaptive Large Neighborhood Search (ALNS). The authors employed a number of insertion and removal heuristics, to intensify and diversify the ALNS search process. In their study, five different variants of the VRP were considered, namely VRPTW, CVRP, MDVRP, the site dependent vehicle routing problem (SDVRP) and the open vehicle routing problem (OVRP). Experimental results from the study [39] show that the proposed algorithm improved 183 best known solutions out of the 486 benchmark tests.
Vidal et al. [30] developed a hybrid genetic algorithm framework for solving the MDVRP and periodic VRP. The proposed hybrid algorithm combines the advantages of the exploration breadth of population-based evolutionary search process, the aggressive improvement capabilities of the neighborhood-based methods, and the advanced population diversity management scheme to improve its solution quality and robustness. The proposed method was tested on standard VRP benchmark instances. Experimental results conducted demonstrate the superiority of the hybrid algorithm in terms of computational efficiency and solution accuracy. Since the new method was able to efficiently identify the best known solutions, including the optimal solutions and new solutions for all the tested benchmark instances, it was concluded that the algorithm was better than the compared techniques.
Juan et al. [10] combined biased randomization with iterated local search to develop a hybrid method for solving the MDVRP with a limited number of identical vehicles per depot. Biased randomization is the use of non-symmetric probability distributions to generate randomness. In [10], two biased-randomized processes were employed at different stages of the iterated local search framework in order to assign customers to depots and also to improve routing solutions. One advantage of their method is that the hybrid algorithm utilizes only one parameter and as such, does not require any parameter tuning. The proposed hybrid algorithm by Juan et al. [10] was tested on a standard benchmark with great success, which makes it an interesting approach for solving MDVRP.
In a related work by Teymourian et al. [29], a framework consisting of enhanced intelligent water drops and cuckoo search algorithms was proposed for solving the capacitated vehicle routing problem. The framework employs four state-of-the-art algorithms to solve the CVRP. These algorithms include an improved IWD algorithm, advanced cuckoo search algorithm, two local search hybrid algorithms, and a post-optimization hybrid algorithm. The proposed hybrid approach takes merit of the capabilities of the two metaheuristics (that is, the improved IWD and advance cuckoo search algorithms) to explore and exploit the solution search spaces. The performance of the proposed algorithms was evaluated on two well-known benchmark instances from literature and their results were compared with some other state-of-the-art algorithms from literature. The experimental results show that proposed methods can effectively cope with large-scale problems, where in most cases the local search hybrid algorithm yields the best gained solutions in the literature.

Metaheuristic solution to MDVRP
The two metaheuristic algorithms employed to develop the proposed hybrid algorithm, are discussed in this section. The IWD metaheuristic algorithm is discussed firstly, followed by the description of the approximate SA algorithm.

Intelligent water drops algorithm
The IWD algorithm is a population-based evolutionary metaheuristic algorithm influenced by the movement of natural water drops finding its way to the river, lakes, or seas. The idea of the IWD algorithm was first proposed by Shah-Hosseini [32] in 2009. Since then, the algorithm has been applied to solve several optimization problems, such as the n-queen puzzle multidimensional knapsack problem [33], multilevel thresholding of gray-level images [33], multiobjective job shop scheduling in scheduling system [40], optimum reservoir operation in water resources systems [41], robot path planning in robotics [42], economic load dispatch problem in power systems [43], feature selection with rough set [44], search and selection optimization processes [45], and examination time-tabling scheduling problem [46]. (Reference could be made to [47] for a comprehensive summary of the various problems that have been successfully solved using the IWD algorithm).
The IWD algorithm problem solving approach is modeled in the form of a graph, G = (V, E), where V and E denote sets of nodes and edges. The structure of the graph depends on the problem representation. The orientation of the problem for which the IWD algorithm is supposed to optimize and find a solution, is usually viewed based on the assumption that there exists a node source from which the water drop moves through a selected path to the next unvisited node. The paths through which the water drops traverses have some loads of soils. Therefore, the choice of selecting a specific path by the IWD is dependent on the amount of soil present on the unvisited path and the path with less soil is usually selected. However, during this process the velocity of the water drops may change, depending on the quantity of soil it is able to offset or accumulate along the selected path of movement. The whole process is repeated iteratively and the best path is updated periodically until the best solution or global solution is found and updated subsequently, after which the algorithm is terminated. The main algorithm procedure and the flowchart are presented in algorithm listing 1 and   Enhanced intelligent water drops algorithm for multi-depot vehicle routing problem There are two types of parameter settings defined for the IWD algorithm, namely static and dynamic parameter settings. Some examples of the static parameters include termination criteria, which determine when the algorithm should be terminated, the initial soil paths and velocity update parameters, which are constant throughout the iterative execution of the algorithm. On the other hand, examples of the dynamic parameters include a list of visited nodes, and the initial amount of soil load on water drop. This type of parameters changes its values with every increment in iteration steps.

Intelligent water drops algorithm for solving MDVRP
To successfully apply the IWD algorithm to solve the MDVRP, the structure or environment of the MDVRP which is denoted by G = (N, E) must be defined. Applying the steps presented in this section help minimize the total dispatching cost of each assigned vehicle that is supposed to service a customer. To start with, a graph with n nodes and n(n−1)/2 directed edges are constructed. This is used by the IWDs as input to construct solutions for the optimization process. In this case a node denotes a customer, while an edge denotes a route to a customer. Every IWD begins its journey by starting from the first or initial node i and terminates it by visiting the end node i + 1 on the graph. Therefore, some of the important factors and parameters employed by the IWD algorithm to solve the optimization process of the MDVRP are discussed.

Parameter initialization
The first step in the solution construction phase is the initialization of all the static and dynamic parameters such as the number of water drops, which in the current case is equal to the number of vehicles, each edge's soil and each IWD's velocity. The static parameters remain constant throughout the program run, and examples of the static parameters include number of IWDs, the initial value of soil, maximum number of iterations, velocity and soil updating parameters. The dynamic parameters are usually initialized at the start of the IWD search process and are subsequently updated throughout the search process. Examples of the dynamic parameters include a list of nodes visited by IWD (which in this case is vehicle k), initial velocity of IWD, and the initial soil load on IWD.
Route building. Route design by the IWD algorithm is such that at each iteration phase of the algorithm, the individual IWDs build a solution for the MDVRP by moving from customer i to the next customer i + 1 according to a defined selection rule, which is explained below. This process is usually evaluated based on the amount of soil present at each customer's node and the route distance to that node. Each IWD is able to keep track of the list of nodes visited. It is equally important to mention here that route construction for all the available customers is done in parallel by the IWDs. For each of the IWD iterations, only one customer is selected by using the probability distribution function given in Eq 1 below. In the proposed IWD implementation, the route building process executed by the IWDs is based on a parallelism processes adopted by the algorithm to construct its solution.
Solution representation. Building a solution for the MDVRP as stated earlier, consists of two stages, namely clustering of customers and routing for each cluster. To retain the simplicity of the IWD algorithm, a straightforward solution representation scheme is employed, which has also been used in existing research [29,48]. For the solution representation mechanisms, a single string is used to represent the clustering and routing stages. The clustering process is formulated by computing the distance matrix given in Eq 8, which is based on the steps described in the multi-depot vehicle routing problem section above. In Fig 3, a solution representation of a typical MDVRP solution scheme is illustrated, where three vehicles (k = 3) serve ten customers (n = 10), and the depot in this case is represented by 0. In this solution representation, customers 2, 3 and 5 are assigned to the same vehicle, which is vehicle 1, and the routing sequence is the same as the customer order, which is in the sequence of customer 2 first, customer 3 second, followed by customer 5. The second vehicle (vehicle 2) visits customers 7 and 4, for which the routing sequence is customer 7 first, then customer 4. The same goes for vehicle 3 that visits customers 1, 9, and 10, with the routing sequence the same as the customer order. The resulting clusters as shown in Fig 3 consist of some list of customers and the advantage of the clustering process is to generate optimum routes for the allocation of different vehicles to customers. In addition, since clustering in this case can be treated as a local search process, the SA algorithm can be employed as a local search metaheuristic to determine the best cost required for each vehicle to travel to a customer in the same cluster. Therefore, a fitness function, which is defined as the total distance covered by each vehicle per cluster, can be used to evaluate the quality of IWD solution. The implementation of SA is detailed in subsequent section below.
The main goal of the solution construction phase is to generate a population of IWD feasible solutions for MDVRP by implementing the IWD hybrid method described in algorithm listing 3 below. At the initial stage of the algorithm implementation, each IWD starts from a depot and selects the next customer to visit within its mapped cluster that consists of a list of candidate customers. The vehicle constraints are updated after the first visit and if these constraints are not violated, the second customer is selected by the IWD. More so, considering the VRP solution representation presented in Fig 3, each individual IWD generate a cluster of feasible solution for the MDVRP by simulating a vehicle and constructing its routes by iteratively selecting customers based on the vehicle's predefined constraints of not exceeding its capacity and visiting each customer once, until all the customers in the feasible solution search spaces (as shown in Fig 4) are serviced. Each IWD follows a similar approach to construct a feasible route for all vehicles. The construction phase is completed by the transition of all IWDs through the feasible solutions, namely that once the vehicle constraints are satisfied and all the customers having been visited, the process is terminated with the IWDs returning back to depot. The solution construction phase is modeled based on an undirected graph with edges and nodes. Therefore, each IWD transverses along the edges from its nodes according to a certain defined probability function explained in the next section.
Selection rule. The selection rule is employed by IWD to select an optimum path i + 1 to a customer or node according to the probability P IWD i i þ 1 ð Þ. In the proposed optimization method, each IWD utilizes a probability distribution function assigned along each edge starting from the current node i through to all other nodes, i + 1, which do not violate the assigned constraints of the problem under consideration. The selection mechanism for choosing an edge connected to the next node by the IWD, is expressed as follows: Let an IWD be at node i, then the probability distribution denoted by P IWD i , which is required by IWD to select an edge to move from node i to node i + 1, is calculated using a fitness function as expressed in Eq (1) as follows.
where the function SM(.) is the saving matrix proposed by [49], which is used to compute the distance traveled by the IWDs along the edges to visit the nodes or as in the current case, the distance traveled by the vehicles for serving the customers. Here, the saving matrix is constructed for every two customers i and i + 1, on the same link path to the given depot h. The parameter ε is a very small positive number assigned to prevent singularity or possible division by zero. The function g(soil(i,i+1)) serves as shift function that moves the soil through an edge joining any two nodes, i+1 towards a positive value. The function is given as follows: where C denotes the set of nodes that IWD is not allowed to visit. Therefore, the edge selection procedure between two nodes i + 1 and k from node can be summarized based on the following conditions:  Enhanced intelligent water drops algorithm for multi-depot vehicle routing problem then select edge arbitrarily to visit any of the connected node Set of visited nodes. An updated list of visited paths C created earlier on, is maintained by IWD as a means of keeping track of all the visited nodes. Iteratively, each node is first evaluated on the basis of whether it has been visited or not, before a decision is taken. If it is confirmed that the node has been visited previously, the node is deleted. The check is performed in order to prevent the IWD from traversing the same path twice. Relating this technique to the MDVRP for example, having customer i on the route of vehicle k from depot h, and after i has been visited by k, it is removed completely from that route and any subsequently visit by k to i is declared taboo for a specific number of iterations. This condition holds for as long as k is only allowed to serve i just once.
Local soil and velocity update. As IWD move from node i to node i + 1, its velocity is updated as follows: where a v , b v and c v are the IWD velocity updating parameters and vel IWD (t) is the previous velocity of the IWD.
The time taken by IWD to move from node i to node i + 1 is computed as follows: The max(.,.) returns the maximum value between its arguments'. This value is then used to threshold the negative velocities to a very small positive number ε v . The function d(.) represents the distance taken by IWD to move from node i to node i + 1.
The amount of soil on the link will be reduced, since IWDs carries some soil as it traverses from node i to node i + 1. Therefore, the soil carried by the IWD as it moves along the link from node i to node i + 1 is updated as follows: where a s , b s , and c s are the IWD soil updating parameters, ρ is a small positive number (0<ρ<1).
Fitness function. The main reason for determining the fitness function is to increase the chances of finding a global-best solution and to also improve the convergence speed of the IWD algorithm. The fitness function determines the ranking of the individual solution obtained by calculating the total length of the constructed route traversed by each IWD in the iterations. The solution with the minimum route length among all the IWD constructed routes is then taken as the best solution. Since this is a minimization case, the route length denoted by T IWD can be expressed as follows: Therefore, the fitness function can be defined as follows: where n is the total number of nodes or customers and the function d(.) is the Euclidean distance between customer is and customer i + 1.
Global update. To prevent IWD from plunging into local minima, the amount of soil on each of the current iteration's best solution with the minimum route length T M , was subsequently updated as follows: However, if at the end of each iteration process, T M is found to be shorter than the current best solution, which is denoted by T B , the best route is updated as follows: Termination condition. The program is terminated once there is no further improvement on the global soil updating, that is after a number of successive iterations have been performed. This also corresponds with the value of the constant parameter referred to as the maximum number of iteration, which is set to 100 in the case of the proposed method.

Simulated annealing
Since its introduction as a solution method into the field of optimization techniques, SA algorithm has been used to solve numerous optimization problems, either on the basis of it being used as a classical algorithm or as part of a hybrid implementation when combined with other metaheuristics. As a local search metaheuristic, the application of SA has produced fairly good results in comparisons with other heuristic based algorithms [50,51]. The SA algorithm is intensely studied in the literature and in some cases it is specifically applied to solve the VRP and its variants problems [20].
The reason for the introduction of SA into IWD was to enhance the existing solution quality of the results produced by the IWD. Another reason for the combination of these two algorithms was to develop a robust local search technique that prevents the IWD from getting stuck into local minima as earlier discussed in the clustering process presented above. SA can be viewed as a search process that can always attempt to move from one current solution to another solution in its neighborhood solutions. Therefore, it has the potential of providing better objective values for the IWD to solve the MDVRP efficiently. Unlike the hill climbing method, SA is able to escape from getting trapped into local minima by allowing worse moves (lesser quality) or uphill steps to be taken at random some of the time. SA choice of selecting best solution is based on its movement procedure, which is such that if the anticipated move is better than its current position, then SA will always take it. If the move is worse, it will be accepted based on some probability selection criteria. For the MDVRP, the SA procedure begins by considering the solution construction of each IWD, given as T IWD i ji ¼ 1; 2; . . . ; n, for a set of given customers with an updated solutions T IWD iþ1 created by randomly switching the orders of two customers. The switching order is implemented based on the concept of 2-opt exchanges of paths between different customers' routes [52,53]. In the 2-opt operator, all possible pairwise exchanges of customers inside an individual vehicle route are evaluated to see if a shorter route distance can be obtained by exchanging the order in which the customers are visited by the vehicle. To briefly illustrate this concept as shown in Fig 5, the following nodes are considered {i, i+1} and {j, j+1} and the process is represented as follows: {i, i+1}{j, j+1}!{i,j}{i+1,j+1}. In Fig 5, the rectangle represents a depot, while the circles denote customers or nodes.
The cost function or fitness function, which represents the quality of solution T IWD i (that is the total distance or route length travelled by the vehicles in a cluster) is denoted by f ðT IWD i Þ. The relative change in cost Δf between T IWD i and T IWD iþ1 is expressed as Since the fitness functions of the solution are computed iteratively, an approximate evaluation of the function that at the same time provides a good indication of the solution quality can be employed. It is equally important to note that a simplified evaluation of the fitness function is necessary to reduce the computation time of the SA.
Beginning with the initial solution, only the solution which has a smaller fitness value than the previous solution is accepted by the algorithm. In other words, a solution is only accepted with the fitness value of f ðT IWD iþ1 Þ < f T IWD i À Á . However, accepting or rejecting a new solution with higher fitness values for T IWD iþ1 is based on the probability of acceptance criteria given as follows: In the case of large problem sizes, instead of using Eq 16, Eq 17 can be employed to improve the performance of the SA: Similarly, Eq 17 can be approximated to Eq 18. This approximation method is applied in this case to reduce the computational effort or time of the SA procedure.
where t τ is the parameter temperature at the τ th instance of accepting a new solution. The probability of accepting a new solution is a function of both the temperature of the system and the difference in the fitness value. It has been noted that the probability of accepting a worse Enhanced intelligent water drops algorithm for multi-depot vehicle routing problem solution decreases as the temperature deteriorates. This means that as the temperature reduces to zero, only better solutions will be accepted. In this paper, the following cooling schedule is adopted (Eq 19): where, α denotes the rate at which the temperature is lowered each time a new solution T IWD iþ1 is discovered. The SA procedure is presented in the algorithm listing 2. Therefore, in this paper, two versions of the hybrid methods are proposed namely, the hybrid method that combines IWD algorithm and SA algorithm using Eq 17, which is denoted here as IWD-SA, and the second proposed hybrid method that combines IWD algorithm with approximated SA algorithm using Eq 18, which is denoted here as IWD-ASA.

The proposed hybrid method
The hybrid IWD-SA algorithm introduces the SA probability of acceptance criteria in determining between the current best cost and the new solution cost, which of them to choose as the solution to be updated. The algorithm achieves this technique by computing and comparing iteratively the quality of solution obtained by both the old and new IWD's drops, which revolve around T M and T B . Therefore, the acceptance rule is evaluated based on two conditions namely, the fitness function or quality of solution, and the environmental temperature. The acceptance criteria enable the IWD process to easily escape from being trapped into local minima, thereby increasing the rate of the algorithm exploration and convergence. The step by step description of how the combined IWD algorithm and SA method work together to achieve the aforementioned process, is presented in Algorithm listing 3. The illustration of these steps is also illustrated using the flowchart shown in Fig 6. Algorithm

Computational experiment
This section presents the various experimental assessments of the IWD algorithms demonstrated in three phases. The first phase is the experiment phase that combines the computational effort of both the IWD algorithm and SA to find the minimum best cost through the optimization of the MDVRP, which is achieved by using Eq 17. The second phase of the experiment demonstrates the results obtained by combining IWD algorithm and the approximate version of the SA algorithm using Eq 18. The last phase of the experiment demonstrates the application and effectiveness of IWD algorithm to solve the MDVRP. The objective of this evaluation is to report the variations in the results emanating from the improvements of the IWD algorithm based on the three empirical scenarios. The results in each method are compared with the best known solution from Cordeau's instances taken from Cordeau et al., [37].
The three experiments were set up and executed under the same processing conditions. The three case scenarios of the IWD and SA algorithms were executed in Windows 7 OS, Intel Celeron CPU@1.8GHz with 3.88GB of RAM. MATLAB R2014b was used as the programming Enhanced intelligent water drops algorithm for multi-depot vehicle routing problem language. The MDVRP instances results described in Cordeau et al. [37], Pisinger and Ropke [39], Vidal et al. [30] and Juan et al. [10], were benchmarked to evaluate the performance of the three methods proposed in this paper. Each instance of the problem was run 10 times for a maximum number of 100 iterations. The parameters used for the implementation of the proposed methods are as presented in Table 1. In the course of the experiment, it was observed that the choice of parameter selection sometimes influences the quality of solution obtained, the higher the values of the parameters being selected, the more computation time it consumes in some cases. This is true in the case of the simulated annealing with an exponential calculation of the probability of acceptance criteria that is computationally expensive. In the analysis of the three methods described in this paper, both the static and dynamic parameters are configured as shown in the parameter Table 1 for the IWD and Table 2 for the SA algorithm. The static parameters are set using similar theoretical values from the work presented in [32].
The first experiment conducted aimed to compare between IWD algorithm and the two versions of the SA, that is, SA with the exponential calculation of the acceptance probability function (IWD-AS) and with an approximate version of the exponential probability of acceptance criteria (IWD-ASA). The results of these tests are summarized in Table 3. Moving column wise in Table 3, the best known solution (BKS) for the Cordeau MDVRP instances is shown, while the average results for the 10 runs is denoted by 'AVE', the best solution found is denoted by 'Best', the standard deviation is denoted by 'Std. Dev', and the average CPU time in seconds is denoted by 'Time' are equally shown. Table 3 shows the results of the 33 Cordeau MDVRP benchmark instances, which were ran on the three employed algorithms, namely IWD, IWD-AS, and IWD-ASA. As can be seen in Table 3, the IWD was the least performed method compared to the results of the best known solutions, while the IWD-AS and IWD-ASA were able to find optimal solutions in all the 33 test instances. The superiority of IWD-AS and IWD-ASA are also demonstrated in some cases where their results outperformed the best known solution, indicated in bold and negative signs in the corresponding 'Best' and percentage gap 'Gap' columns.

Experimental evaluation of proposed algorithms
Also, to statistically determine the differences among the three algorithms' solution qualities, two measurement tests were computed, namely the resulting gaps between the three proposed algorithms' best solution found, and the best known solutions of the other techniques. The second test is the student's test known as the t-test, which was used to compare the average results of the IWD-SA and IWD-ASA. The gaps were calculated using Eq 20.
While the t statistics was calculated using the formula given in Eq 21.
S " where " X 1 = Average value for the set of IWD results Enhanced intelligent water drops algorithm for multi-depot vehicle routing problem " X 2 = Average value for the other methods results S 1 = Average value for the other methods results S 1 = Standard deviation of the other methods including IWD-SA and IWD-ASA n 1 = frequency of the IWD results for the 10 runs n 2 = frequencies of the other methods for the 10 runs The essence of the t-test is to investigate whether there is any significant difference among the average results of the IWD and the other two methods (i.e. IWD-SA and IWD-ASA). In each case the two-tailed P-values obtained were less than 0.0001 and by conventional criteria these differences are considered to be extremely statistically significant. A confidence interval at the 95% confidence level (t 0.05 = 18) was stated. From the results of the t-tests displayed in Table 3, it can be seen that there is statistically significant differences between the average results of the IWD and the other two techniques. In all the test cases, it was obvious that the IWD is outperformed by the IWD-SA and IWD-ASA based on the compared average values of the three techniques. This significant difference in the performance of the other two methods can be attributed to the local search and hill climbing characteristics of the SA introduced into the IWD. These two features of the SA are efficient mechanisms which can assist the IWD to escape being trapped into local minima, and it also increases the explorative and exploitative power of the IWD within the given solution space. One other factor that tends to improve the performance of the hybrid IWD is the efficient calculation of the cost function without compromising the quality of the solution by the SA. This is very important, since the problem cost function is computed at every iteration run of the algorithm. In the current implementation, this bottleneck was avoided by evaluating the cost function based on the difference between the current solution and the neighborhood solution.
Similarly, in the implementation of the proposed methods, it was observed that the calculation of the exponential probability function alone took approximately one third of the whole algorithm execution time. This however, justifies the initial claim that the computation of the acceptance probability function consumes a copious number of system resources. More Enhanced intelligent water drops algorithm for multi-depot vehicle routing problem Enhanced intelligent water drops algorithm for multi-depot vehicle routing problem specifically, CPU time resulting from the exponential calculation involved in determining the probability of accepting or rejecting a new solution. Similar observations were also pointed out in [54,55]. Thus, this compelled the researcher to further compare the execution times of the two SA probabilities of acceptance criteria models presented in Eqs 17 and 18. From the experimental results presented in Table 3, it is obvious that IWD-SA took longer CPU time than IWD-ASA. The average execution times of 1,145.72 and 949.66 seconds were recorded for IWD-SA and IWD-ASA respectively, while an average of 927.21 seconds was recorded for the IWD algorithm. The variations in the CPU execution time in seconds for the two methods on some selected problem instances are demonstrated in Figs 7, 8 and 9. However, there were no noticeable significant differences in the solution quality of the two SA techniques aside from the observed wide variation between their processing speeds.

Experimental evaluation of the algorithms with other methods
This section presents the best solution obtained by each of the proposed methods for comparison with some other approaches available in literature. The gaps in percentages for each of the best results obtained by the different methods are presented for evaluation in Tables 4-6 and similarly, illustrated visually in Figs 10-12. The proposed methods were able to achieve 19 matches out of the 33 existing best known solutions, compared to the results reported in [10], in which 12 matches out of the 33 previous best known solutions were recorded. Similarly, in comparing the average percentage gaps of the three proposed methods, it can be seen that, they are reasonably low, except for the IWD technique. For example, the average gaps between IWD-SA and other techniques are obtained as -0.10% for comparison between IWD-SA and Cordeau et al., [37], 0.17% for comparison between IWD-SA and Pisinger and Ropke [39], 0.27% for comparison between IWD-SA and Vidal et al. [30], and 0.01% for comparison between IWD-SA and Juan et al. [10] for all the 33 instances.

Statistical and convergence behavior analysis
In addition to the previous comparison made, two statistical tests were carried out, using the results presented in Tables 4, 5 and 6 in order to obtain more rigorous and fair conclusions regarding the performance and convergence behavior of each of the compared algorithms. Firstly, Friedman's non-parametric test was conducted on the five different algorithms and the result is presented in Table 7. Secondly, an application was conducted of post hoc analysis with Wilcoxon signed-rank tests using IWD, IWD-SA, and IWD-ASA as controlled algorithms the result of which is presented in Table 8. Friedman's non-parametric test for multiple comparisons conducted indicated that there was a statistically significant difference among the results of the five (5) algorithms. The test was conducted by taking into account the confidence interval which was stated at the 95% confidence level. As shown in Table 7, Friedman's test revealed that there were statistically significant differences among the techniques whilst running, X 2 (4) = 82.68, p-value = 0.000 for the comparison between IWD and other methods; X 2 (4) = 37.60, p-value = 0.000 for the comparison between IWD-SA and other methods; and X 2 (4) = 37.52, p-value = 0.000 for the comparison between IWD-ASA and other methods. Since the Friedman test result was statistically significant, there was a need to still run the post hoc tests, so that it could be ascertained exactly which of the algorithms performed better.
To properly evaluate the statistical significance of the efficient performance of the proposed IWD, IWD-SA and IWD-ASA algorithms, the post hoc tests have been conducted further by using each of the three methods as the controlled algorithm in all the comparisons made. However, the post hoc analysis with Wilcoxon signed-rank tests conducted with a Bonferroni correction applied, resulted in a new significance level set at p < 0.0125 (i.e. 0.05/4, since 4 comparisons were made). Analyzing the test statistics data shown in Table 8, it can be seen that at p < 0.0125 significance level, the proposed methods produced results that were significantly better than Cordeau et al. [37] and Vidal et al. [30]. However, there were no statistically significant differences between the results obtained by the two proposed methods (that is, IWD-SA and IWD-ASA) (see Pisinger and Ropke [39], and Juan et al. [10]), since their resulting p-values were greater than the computed p-value of 0.0125. Therefore, it can be concluded that the two proposed methods are statistically significantly better than the IWD and Cordeau et al. [37] for solving MDVRP problems at 95% confidence level. Table 7 reports the average ranks returned by Friedman's non-parametric test for the five set of algorithm on 33 MDVRP instances. As shown in Table 7, Vidal et al. [30] yielded the best performance among all the other methods, while two of the current methods namely,  Evaluating the different methods proposed in this paper, is based on each algorithm's solution qualities and speed of convergence, as shown in Figs 16 and 17. These two figures depict that both the IWD-SA and IWD-ASA competed favorably well in terms of solution quality and speed of convergence compared to the standard IWD, despite the additional computational cost incurred by the IWD-SA. This further support the significant contributions of the different strategies employed to implement the proposed algorithms. Some of these main strategies include the use of efficient cost function, approximation of SA's probability of acceptance criteria and the use of SA as an alternative local search technique to intensify the entire search process, and at the same time to accelerate the convergence speed of the proposed hybrid algorithms.
In summarizing the whole analysis, from the computational result obtained in Table 3, it is concluded that there is a significant difference in the CPU execution time between IWD-SA and IWD-ASA, with the approximate version of the SA having the least execution time. This reduction in the algorithms' run time indirectly resulted in accelerating the convergence speed of the IWD-SA without compromising its solution quality. Therefore, it can be stated that the approximate method of the hybrid algorithm (IWD-ASA) performed better than its IWD-SA counterpart in terms of CPU time. After analyzing the results presented in Tables 4, 5, and 6, it is concluded that the two proposed hybrid algorithms achieved high quality solutions compared to other methods described in the literature (19 out of 33). Specifically, it is also concluded that the two hybrid methods performed well in MDVRP instances with up to 360 nodes (customers). In summary, as shown in all the results evaluated, the proposed methods can be employed as alternative hybrid metaheuristics for solving other similar MDVRP problems.

Conclusion and future research
In this paper the multi-depot vehicle routing problem has been studied. Presented are two hybrid metaheuristics that incorporate both the exponential and approximate version of the Enhanced intelligent water drops algorithm for multi-depot vehicle routing problem simulated annealing into the basic intelligent water drops algorithm. This study utilized these two methods of the SA solution to develop a better alternative algorithm for the MDVRP NPhard problem. The proposed algorithms (IWD-SA and IWD-ASA) were tested on two sets of Cordeau's MDVRP benchmark instances, covering instances P01-P23 and Pr01-Pr10. Two different statistical analyses were also conducted to verify the performances of the proposed  Enhanced intelligent water drops algorithm for multi-depot vehicle routing problem A statistical analysis was also conducted to compare the performances of the proposed methods. Results of the analysis revealed that there are no significant statistical differences in terms of solution qualities obtained by the two SA methods. However, statistical results also Enhanced intelligent water drops algorithm for multi-depot vehicle routing problem Enhanced intelligent water drops algorithm for multi-depot vehicle routing problem  Table 8. Application of post hoc analysis with Wilcoxon signed-rank tests using IWD, IWD-SA, and IWD-ASA as controlled algorithms.

Test Statistics b
IWD-Cordeau et al. [37] IWD-Pisinger and Ropke [39] IWD-Vidal et al. [30] IWD-Juan et al. [  revealed that the exponential implementation of the SA consumed more CPU time than the approximate method. For example, the IWD-SA took an average of 1,145.72 seconds to run all the 33 benchmark instances considered in the implementation, in comparison with 949.66 and 927.21 seconds for IWD-ASA and IWD. The basic IWD was outperformed by the two IWS-SA and IWD-ASA hybrid methods, which is a good indication that SA actually contributed to improving the solution quality of the standard IWD algorithm. Comparisons with other methods show that there are significant statistical differences among the compared algorithms from the literature with the proposed IWD-SA based methods, while the classical IWD Enhanced intelligent water drops algorithm for multi-depot vehicle routing problem was outperformed by all the other algorithms using the same parameter settings with the SA variants. Therefore, these results indicate that both IWD-SA and IWD-ASA are more effective alternative candidate algorithms than the compared methods, for solving multi-depot vehicle routing problem. As a direction for further research, the proposed methods could be further improved and subsequently applied to solve other variants of capacitated vehicle routing problems. Also a better approach that could be exploited is to build a look-up table for a set of values over the range Δf/t to be used as alternative to the exponential calculation of the probability of acceptance criteria.