# A Hybrid Ant Colony Optimization Algorithm for Multi-Compartment Vehicle Routing Problem

Authors: Ning Guo; Bin Qian; Rong Hu; Huai P. Jin; Feng H. Xiang
Journal: Complexity
 (2020)

---

## Abstract

The multi-compartment vehicle routing problem (MCVRP) has been applied in fuel or food delivery, waste collection, and livestock transportation. Ant colony optimization algorithm (ACO) has been recognized as an efficient method to solve the VRP and its variants. In this paper, an improved hybrid ant colony optimization algorithm (IHACO) is proposed to minimize the total mileage of the MCVRP. First, a probabilistic model is designed to guide the algorithm search towards high-quality regions or solutions by considering both similar blocks of customers and customer permutations. Then, a heuristic rule is presented to generate initial individuals to initialize the probabilistic model, which can drive the search to the high-quality regions faster. Moreover, a new local search using the geometry optimization is developed to execute exploitation from the promising regions. Finally, two types of variable neighborhood descent (VND) techniques based on the speed-up search strategy and the first move strategy are devised to further enhance the local exploitation ability. Comparative numerical experiments with other algorithms and statistical analyses are carried out, and the results show that IHACO can achieve better solutions.


---

## Body

## 1. Introduction

The multi-compartment vehicle routing problem (MCVRP) is a variant of the well-known vehicle routing problem (VRP) [1]. Compared with the standard VRP, the MCVRP can simultaneously transport different types of products, which have to be placed in separate compartments. Since the customers probably have different demands for various products whose individual characteristics are so different that they cannot be mixed, each vehicle is divided into several compartments with a certain capacity for satisfying customers’ needs. There are many industrial fields in which the multi-compartment vehicles are employed. The first example of the MCVRP arises in the process of supplying the fuels, where a lot of vehicles or ships with some tanks of various capacities are used to settle this problem [1–7]. Another example is the transportation of food, where different degrees of refrigeration goods are stored in different compartments in one vehicle [8–11]. In addition, the MCVRP has also been widely applied in the waste collection, where separate compartments are required to convey different kinds of garbage [12–14]. Despite the practical applications, the solution of MCVRP has been rarely studied. Thus, this paper aims to minimize the total mileage for such a problem.Since the VRP is an NP-hard problem and it can be reduced to the MCVRP, the latter issue is also NP-hard. Some operational research algorithms were presented to resolve small-scale problems for the MCVRP [1, 15–20]. However, with the increasing problem size, these algorithms are usually of limited use owing to their excessive time-consuming nature. Then, a variety of meta-heuristic algorithms were proposed to obtain desirable solutions within a reasonable calculation time [21–24]. There are many intelligent algorithms for solving nonlinear problems [25–27]. The ant colony optimization algorithm (ACO) is one of the highly efficient meta-heuristic or intelligent algorithms. The ACO is inspired by the intelligent behavior of ants in searching for food, during which the ants leave the so-called pheromone on the road. The ants exchange information with other ants through the pheromone. The higher the number of pheromones on a path, the more likely the ants are to go after that path, and vice versa. In this way of the pheromone positive feedback mechanism, the ants can find the shortest path to the food location.Because the ACO has a strong global exploration and distributed computing capability, it has been successfully applied to handle different kinds of optimization problems. While solving the VRP and the complex problems derived from it, ACO has been usually adopted by many researchers. The existing results usually use only pheromone concentration matrix (hereinafter referred to as pheromone matrix) to construct a probabilistic model to generate new individuals. The pheromone matrix learns and accumulates similar blocks of customers from every outstanding individual. Here any two successive visited customers in an individual compose one block. Mavrovouniotis and Yang [28] proposed an ACO with immigrant schemes for the dynamic vehicle routing problem. In this ACO, only the pheromone matrix was utilized to get information of each block, and its initial value was set to 1/Cnn, where Cnn is the corresponding path length of a solution constructed by a nearest neighbor heuristic. Schyns [29] presented an ACO to settle the responsive dynamic vehicle routing problem with time window. The pheromone matrix was adopted merely and initialized to 1/nL0, where n is the number of customers and L0 is the length of the route formed on the basis of the FIFO (first in first out) principle. Besides, a local search basic heuristic was developed to improve the algorithm performance. Wang et al. [30] designed a novel ACO called AMR to solve the vehicle routing problem. In AMR, only the pheromone matrix was used and its initial value was a small positive constant denoted by phmin. In addition, the saving algorithm was introduced to enhance the performance of AMR. Huang et al. [31] employed an ACO combined with a local search method in the resolution of the routing problem with bridge inspection tasks. Only the pheromone matrix was applied, and the initial pheromone was set to 1/N⋅Jψini, where N is the total number of nodes and Jψini is the cost of the initial solution obtained by the first ant. Moreover, ACO has been rarely used to solve the MCVRP. Reed et al. [21] added a technique of clustering and local search of 2-opt to ACO in their study. Only the pheromone matrix was utilized and initialized to 1/nLneig, where n is the customers’ number and Lneig is the total length of the nearest neighbor path. Motivated by the work of Reed et al. [21], Abdulkader et al. [22] proposed a hybridized ant colony (HAC) algorithm with local search schemes. The pheromone matrix was adopted alone, and the initial value was set to 1/L, where L is the length of a path which was generated randomly.As it can be seen from above literatures, the probabilistic model only consists of the pheromone matrix in which only the relationship between two consecutive customers is considered. However, as mentioned earlier, the ACO generates new individuals based on the probabilistic model which stores the information of excellent individuals. Obviously, the more comprehensive the information of excellent individuals are learned, the easier it is to guide the algorithm to reach different high-quality regions as soon as possible. Hence, only the information accumulation of similar blocks in the customer sequence is incomplete and limited. Besides, although pheromone matrices in different literatures above are initialized to different values, the values in the pheromone matrix of each literature are the same. This means that, at the beginning of the algorithm, there is no valuable information in the pheromone matrix. With the iteration of the algorithm, when the information in the pheromone matrix accumulates to a certain amount, it will guide new individuals towards high-quality areas. Invisibly, the efficiency of the algorithm is reduced in this way. Moreover, it is also found that the introduction of local search on the basis of the effective neighborhood operation can further obtain higher performance of the algorithm. Hence, in this research, an improved hybrid ant colony optimization algorithm (IHACO) is devised to solve the MCVRP. The salient features of this IHACO are as follows:(1)
A new matrix is designed to save the whole order information of customers. Together with the pheromone matrix, the probabilistic model of IHACO is formed to record the worthy information of optimal individuals. Compared with the previous work, this design can learn and retain valuable information more comprehensively.(2)
Several initial individuals produced by a heuristic rule are used to initialize the probabilistic model so that some valuable information can be utilized at the beginning of the algorithm. Then, the speed of reaching high-quality solution regions will be accelerated, and thus, the efficiency of the algorithm is raised.(3)
Different from the majority of current local searches where 2-opt is commonly used within the same route, a new local search is proposed by adopting the geometry optimization method to perform narrow but productive exploitation.(4)
In the IHACO’s local search, three techniques are devised to enhance the exploitation capability from many aspects: two kinds of VND are used to increase the search depth of the algorithm, a speed-up search strategy based on the problem’s properties is used to accelerate the evaluation process of individuals, and a first move strategy is used to enrich search regions.The remainder of this paper is organized as follows: Section2 introduces the mathematical model for the MCVRP; Section 3 describes the proposed IHACO for details; Section 4 presents computational experiments and analyses; and finally, Section 5 gives some conclusions and suggestions of future work.


## 2. Problem Description and Model Formulation

The MCVRP is defined on an undirected graphG=N,A with a set of vertices N=0,1,…,n and a set of arcs A=i,j:i,j∈N,i≠j. Denote dij as the distance from i to j. Node 0 in N represents the depot, and Nc=N/0 are nodes of customers. The depot has an available homogeneous fleet. These identical vehicles start from the depot, visit the customers distributed in the service area, and go back to the depot in the end. All vehicles are equipped with multiple compartments whose number equals to the number of products. Each customer i∈Nc has the quantity qip to be carried for the product p, and the compartment capacity reserved for this product does not exceed Qp. Moreover, each customer is visited exactly once by only one vehicle.To formulate the model, some sets, parameters, and variables adopted in the proposed problem are listed in Table1.Table 1
Notations applied in the MCVRP model.
SymbolDescription of symbolNNode setN=0∪NcNcCustomer subsetPProduct setdijDistance between nodesi and jqipQuantity of productp to be carried from node iQpCapacity of the vehicle compartment dedicated to productpxijEquals 1 if the vehicle travels from nodei to j and is zero otherwiseQipTotal quantity of productp while the carried vehicle departing from node iThen, the integer linear program of the MCVRP can be described as follows:(1)minimize⋅∑i∈N∑j∈Ndijxij,subject to(2)∑j∈Nxij=1,∀i∈Nc,(3)∑i∈Nxij=1,∀j∈Nc,(4)∑i∈Ncx0i=∑j∈Ncxj0=1,(5)qip≤Qip≤Qp,∀i∈Nc,p∈P,(6)Qip−Qjp+Qpxij≤Qp−qjp,∀i∈Nc,p∈P,(7)xij∈0,1,∀i,j∈N,i≠j.Equation (1) is the objective function indicating the total mileage for the problem. Constraints (2) and (3) imply that one customer is visited just once by only one vehicle. Constraint (4) guarantees that each vehicle starts and ends its trip at the depot. Constraint (5) ensures that the total quantity delivered by each vehicle after serving customer i is less than its capacity for this product. Constraint (6) imposes the connectivity requirement of the feasible routes. Both inequalities (5) and (6) are subtour elimination constraints. Constraint (7) describes the domains of the decision variable.


## 3. IHACO for the MCVR

This section describes the implementation of the proposed IHACO in detail, including the solution representation, establishment, and initialization of the probabilistic model, formation of the individuals, update of probability matrices, speed-up VND-based local search, and the complete IHACO process.


### 3.1. Solution Representation

Usually, the encoding rules of ACO for solving VRP are that each ant represents one vehicle, the path of each ant represents the driving path of each vehicle, the depot is represented by the number 0, and the customer is represented by a natural number [32]. In this paper, the decimal coding based on customer arrangement is still used, but the path of each ant which is denoted πant represents the path of all customers served by the depot, and the depot is a separator for different vehicles. For example, the path sequence of the ant i is πanti=0,1,5,3,7,0,6,2,8,4,0 which means that vehicle 1 leaves the depot 0, visits customers 1, 5, 3, and 7 in turn, and then returns to the depot 0; vehicle 2 leaves the depot 0, visits customers 6, 2, 8, and 4 in turn, and finally returns to the depot 0.


### 3.2. Probabilistic Model and Its Initialization

The performance of ACO has strong connections with the probabilistic model. A better choice of the model is critical to improve the efficiency of the algorithm. In this research, two probability matrices are employed to form the probabilistic model. Besides, a few good individuals produced by a heuristic method are used to initialize these two probability matrices. This initialization method makes the algorithm converge to the high-quality solutions faster.


#### 3.2.1. Probabilistic Model

In the traditional ant colony algorithm, the probabilistic model usually refers to ants’ pheromone matrix which is obtained by calculating the pheromone concentrationτij. This pheromone value represents the probability that the customer j appears immediately after the customer i. In other words, the pheromone matrix only reflects similar blocks of customers presented in the chosen individuals. But it is not enough to accumulate this kind of information alone. Counting the chance of customer j appearing before or in the current position is also important. Hence, a new probability matrix based on the order of customers called sequence information matrix is established. Then, the probabilistic model in this paper contains a pheromone matrix and a sequence information matrix.In a subpath visited by one vehicle, letskj be the appearance times that customer j before or in position k in the chosen sequences. Next, the specific probability value in the sequence information matrix is determined by(8)μkj=skj∑l∈Nksk,l,where Nk is the set of all available customers not already served by the current vehicle until position k. Here is an example to compute μkj. Suppose that the path taken by three ants are πant1=0,1,2,3,0,4,0, πant2=0,2,3,4,0,1,0, πant3=0,3,4,0,1,2,0, and πant4=0,3,2,1,4,0, and then, skj is given as follows:(9)skj5×5=−−−−−−3121−3432−4443−4444.Since the first location is fixed as the depot, the corresponding probability values ofμkj for the second position of each customer in N1=1,2,3,4 are μ21=3/3+1+2+1=0.43, μ22=1/3+1+2+1=0.14, μ23=2/3+1+2+1=0.29, and μ24=1/3+1+2+1=0.14.Then, the pheromone matrix is normalized by formula (10). Finally, let σkj denote the probability of customer j at position k which is calculated by formula (11). Based on this value, the ant selects the next customer to visit:(10)ζij=τij∑l∈Nkτi,l,(11)σkj=μkj+ζij,i=2,3,…,n.


#### 3.2.2. Initialization of the Probabilistic Model Based on a Heuristic Rule

As mentioned in the introduction, normally, the probabilistic model (that is the pheromone matrix in ACO) is initialized to the same value at the beginning of the algorithm. In order to achieve good solutions more quickly, this paper adopts the initial solutions to initialize the proposed probabilistic model but not the same value.Sort the distance between all customers and the depot from small to large, and take the topt=roundn/10 customers into set cfirst. Let cfirst=c1first,c2first,…,ctfirst denote the first customer to be visited by the t-th ant, R=R1,R2,…,Rt the complete path passed by the t-th ant (also the t-th initial solution), and L=L1,L2,…,Lt the length of the route taken by the t-th ant. The procedure of the probabilistic model initialization is given as follows:Step 1: sett=1 and τij0=zerosn×n and skj0=zerosn×n.Step 2: the ant starts from the depot0 and visits the first customer ctfirst. Then, the ant selects the remaining customers for access based on the nearest neighbor principle. This process forms the t-th initial solution Rt.Step 3: calculate the route lengthLt corresponding to Rt.Step 4: setskjt=skjt−1+skjRt.Ift=1, then set τij1=1/L1; otherwise, setτijt=τijt−1+1/Lt1.Step 5: sett=t+1.Step 6: ift≤roundn/10, then go to step 2.Step 7: outputskjt and τijt.


### 3.3. Ant Path Construction

In this study, each ant produces a complete solution for the MCVRP. The ant begins from the depot and continuously chooses customers to its tour as long as the constraint of capacity is satisfied. If there are no feasible customers to select, the ant comes back to the depot and restarts another trip.Letm denote the amounts of ants and t the amounts of the first customers selected. Set t=0,1,…,m. To guarantee the diversity of individuals, the first customer is picked randomly from n−t customers. For example, the first ant randomly selects one of the n customers as the first one to visit. Then, the first customer to be visited by the second ant will be randomly selected from the remaining n−1 customers. If m>n, when t=n+1, reset t to zero.Afterwards, each ant proceeds to add customers to its route based on the probability function (12). The probability of choosing customer j as the next one in the t-th iteration is given by formula (13):(12)δij=σkjαηijβ.In formula (12), σkj contains information about how often the customer j appears in the current position k, ηij is the inverse of distance between customer i and j, and α and β denote the importance of σkj and ηij, respectively:(13)pijt=δij∑l∈Nkδil,ifj∈Nk,0,otherwise.Setq0 be a constant value and q a random variable which is uniformly distributed in [0, 1]. If q>q0, select the next customer in accordance with the formulation of pijt; otherwise, select the customer who has the maximum value of δij [22].


### 3.4. Update Mechanism of Probability Matrices

The probabilistic model consists of two probability matricesskj and τij. They are updated by the best solution found at present. Let ρ0≤ρ<1 be the pheromone persistence, then equation (14) is used to update the pheromone concentration in the path. Meanwhile, equation (15) is adopted to update the probability matrix based on the order of customers:(14)τijt+1=ρ×τijt+1/Lbestt,if arci,j∈Lbestt,ρ×τijt,otherwise.In equation (14), Lbestt is the total length corresponding to the optimal ant path in the t-th iteration:(15)skjt+1=skjt+skjπbestt.In equation (15), skjπbestt is the probability matrix skj obtained by the t-th optimal individual πbestt. We refer to Section 3.2.1 for more details.


### 3.5. Speed-Up VND-Based Local Search

Generally, for the combinatorial optimization problems, e.g., MCVRP, there are many other high-quality solutions near the high-quality solutions found so far. The algorithm in this section performs the problem-dependent local search on the global optimal solutionπbestgen, for the purpose of realizing the detailed search of the area near πbestgen. In this paper, the operations of insert, swap, and intersection elimination based on VND, the speed-up search strategy, and the first move strategy are adopted.


#### 3.5.1. Basic and Pipe VND

The VND is a variant of the variable neighborhood search. The basic VND [33] explores the list of neighborhoods in a sequential way one by one. Once an improvement is detected in any neighborhood structure, the search goes back to the ﬁrst neighborhood from the list. Until all the neighborhoods encounter the local optimum, the search ends. However, the pipe VND [33] executes continuous search in the same neighborhood if an improvement is obtained. In this paper, both these VNDs are applied.


#### 3.5.2. Operation of Insert and Swap

Assume an ant’s tour is represented asπantold=1,…,i−1,i,i+1,…,j−1,j,j+1,…,n. The insert operator selects one customer like i and inserts it before or after another customer like j. Then, the tour becomes πantnew_insert=1,…,i−1,i+1,…,j−1,i,j,j+1,…,n or πantnew_insert=1,…,i−1,i+1,…,j−1,j,i,j+1,…,n. The swap operator selects two customers like i and j on this tour and swaps them. Then, πantnew_swap=1,…,i−1,j,i+1,…,j−1,i,j+1,…n.


#### 3.5.3. Operation of Intersection Elimination

Since the objective function of the MCVRP is constructed to minimize the total length of the path, the optimal path of each vehicle should not be crossed.Theorem 1.
In the route of every vehicle, if two lines intersect, the total path length of the vehicle is larger than that of the total length without intersection.Proof.
As shown in Figure1, there is one vehicle route denoted by πant_vehicleold=0,1,3,2,4,0. For any two given intersecting line segments, such as lines L13 and L24 which intersect at point u, the end points of the two segments that can form the opposite side are connected like L12 and L34. Next, the new path without intersection is obtained and denoted by πant_vehiclenew=0,1,2,3,4,0. According to the theorem that the sum of any two sides in a triangle is greater than the third side, in Δ1u2, L1u+L2u>L12; in Δ3u4, Lu3+Lu4>L34. Let Ltotalold be the total length of vehicle path with intersecting lines, then Ltotalold=L01+L1u+Lu3+L32+L2u+Lu4+L40. Let Ltotalnew be the total length after eliminating intersecting lines, then Ltotalnew=L01+L12+L23+L34+L40. As a result, Ltotalnew−Ltotalold=L12−L1u+L2u+L34−Lu3+Lu4<0.
Fromπant_vehicleold and πant_vehiclenew defined above, it can be seen that if two paths intersect, there are four customers involved. According to the order of customers visited by the vehicle, only exchange positions of the second and third customers and then cross paths can be eliminated. In this paper, the geometric method is used to determine whether the crossing occurs.
Suppose the coordinates of customers 1, 3, 2, and 4 in Figure1 are x1,y1, x2,y2, x3,y3, and x4,y4, and then, the linear equations of line segments L13 and L24 are shown in the following equations:(16)x=x1+λx2−x1,y=y1+λy2−y1,(17)x=x3+ψx4−x3,y=y3+ψy4−y3,(18)x2−x1λ+x3−x4ψ=x3−x1,y2−y1λ+y3−y4ψ=y3−y1.
In formulas (16) and (17), λ is the slope of the line where L13 is located, and ψ is the slope of the line where L24 is located.
IfL13 and L24 intersect, then formula (18) has a solution. In other words, inequality (19) should be satisfied. Meanwhile, it is necessary to satisfy the inequality (22) in order to make two line segments intersect within the line segment (this is the case of path intersection shown in Figure 1), rather than the vertex of the line segment or its extension line. However, equations (20) and (21) are used to calculate λ and ψ, respectively:(19)Δ=x2−x1x3−x4y2−y1y3−y4≠0,(20)λ=1Δ×x3−x1x3−x4y3−y1y3−y4,(21)ψ=1Δ×x2−x1x3−x1y2−y1y3−y1,(22)0<λ<1,0<ψ<1.Figure 1
Route of one vehicle.


#### 3.5.4. Speed-Up Search Strategy

Letinsertπ,i,j, swapπ,i,j, and elim_intersectionπ,i,j represent the operations of inserting, swapping, and uncrossing the i-th and j-th elements in solution π, respectively. Let IJinsertinfeasible denote the set of i,j whose neighborhood is infeasible (the infeasibility here means that the vehicle capacity constraint is not satisfied) after executing insertπ,i,j, and IJswapinfeasible is the set of i,j whose neighborhood is infeasible after executing swapπ,i,j. Because the intersection elimination operation is only carried out within the customers served by each vehicle, there is no neighborhood that violates the capacity constraint. Let Ninsertπ,i denote the feasible neighborhood set based on insert operation for i customers in solution π, Nswapπ,i the feasible neighborhood set based on swap operation for i customers in solution π, and Nelim_intersectionπ,k the feasible neighborhood set based on intersection elimination operation for customers serviced by the k-th vehicle in solution π. Set nk to be the number of customers visited by the k-th vehicle, then l=∑i=1kni−nk+k+1 is the position of the first customer visited by the k-th vehicle in solution π. Hence, there are as follows:(23)Ninsertπ,i=πi,j=insertπ,i,j|i=1,2,…,n−1,j=i+1,…,n,andi,j∉IJinsertinfeasible,(24)Nswapπ,i=πi,j=swapπ,i,j|i=1,2,…,n−1,j=i+1,…,n,andi,j∉IJswapinfeasible,(25)Nelim_intersectionπ,k=πi,j=elim_intersectionπ,i,j|i=l,j=i+1,…,i+nk−1.It can be seen from (23)–(25) that Ninsertπ,i is a feasible neighborhood set that any customer i inserts before or after any possible customer j within the same vehicle or between different vehicles. Nswapπ,i is a feasible neighborhood set that any customer i swaps with any possible customer j within the same vehicle or between different vehicles. Nelim_intersectionπ,k is a feasible neighborhood set that customers involved in the cross-line can be interchanged within the k-th vehicle.For the MCVRP, if customerj is too far away from customer i, the total route distance after the insert and swap operation will be greater than before, which makes the two operations meaningless. Hence, in order to avoid the search of invalid areas, the distances between all possible customers j and customer i are arranged from small to large, and just, the first S points are taken out to insert and swap. Thus, the more compact neighborhood sets Ninsertnarrow1π,i and Nswapnarrow1π,i can be obtained:(26)Ninsertnarrow1π,i=πi,j=insertπ,i,j|i=1,2,…,n−1,j∈S,andi,j∉IJinsertinfeasible,Nswapnarrow1π,i=πi,j=swapπ,i,j|i=1,2,…,n−1,j∈S,andi,j∉IJswapinfeasible.For the sake of speeding up the searching process and making the neighborhood search focus on high-quality areas, the nature of the problem can be utilized to further decrease invalid operations and reduce the feasible neighborhood set. LetLπ be the value of the objective function which is the total mileage of the route corresponding to solution π and dπi,πj be the distance between element πi which is located in position i of π and element πj which is located in position j of π.ForNinsertnarrow1π,i, if customeri is inserted before customer j, Lπi,j−Lπ=dπi,j−dπ=dπi−1,πi+1+dπj−1,πi+dπi,πj−dπi−1,πi+dπi,πi+1+dπj−1,πj; if customer i is inserted after customer j, Lπi,j−Lπ=dπi,j−dπ=dπi−1,πi+1+dπj,πi+dπi,πj+1−dπi−1,πi+dπiπi+1+dπj,πj+1. For any i,j, if dπi,j−dπ>0, then Lπi,j>Lπ. i,j in this case can be directly excluded without calculating Lπi,j. If Lπi,j<Lπ, then Lπi,j can be obtained directly from Lπi,j=Lπ+dπi,j−dπ instead of starting from the first customer of π to calculate. Therefore, a more compact and effective neighborhood set Ninsertnarrow2π,i can be obtained on the basis of Ninsertnarrow1π,i:(27)Ninsertnarrow2π,i=πi,j=insertπ,i,j|i=1,2,…,n−1,j∈S,i,j∉IJinsertinfeasible,anddπi,j−dπ<0.ForNswapnarrow1π,i, Lπi,j−Lπ=dπi,j−dπ=dπi−1,πj+dπj,πi+1+dπj−1,πi+dπi,πj+1−dπi−1,πi+dπi,πi+1+dπj−1,πj+dπj,πj+1. And for any i,j, if dπi,j−dπ>0, then Lπi,j>Lπ. i,j in this case can be directly excluded without calculating Lπi,j. If Lπi,j<Lπ, then Lπi,j can be obtained directly from Lπi,j=Lπ+dπi,j−dπ instead of starting from the first customer of π to calculate. Therefore, a more compact and effective neighborhood set Nswapnarrow2π,i can be obtained on the basis of Nswapnarrow1π,i:(28)Nswapnarrow2π,i=πi,j=swapπ,i,j|i=1,2,…,n−1,j∈S,i,j∉IJswapinfeasible,anddπi,j−dπ<0.ForNelim_intersectionπ,k, we can use the method described in Section 3.5.3 to determine whether two paths intersect. If so, according to Theorem 1, we know Lπi,j−Lπ=dπi,j−dπ=dπi−1,πj+dπi,πj+1−dπi−1,πi+dπj,πj+1<0. Then, Lπi,j can be obtained directly from Lπi,j=Lπ+dπi,j−dπinstead of starting from the first customer of π to calculate.LetSP_FindBestNinsertnarrow2π,i, SP_FindBestNswapnarrow2π,i, and SP_FindBestNelim_intersectionπ,k be the optimal neighborhoods in searching and outputting Ninsertnarrow2π,i, Nswapnarrow2π,i, and Nelim_intersectionπ,k, respectively. These three optimal neighborhoods utilize the above speed-up search method in the local search of IHACO. Meanwhile, the pipe VND described in Section 3.5.1 is also applied. The corresponding procedures are given as follows:Step 1: seti=1 and πbestprevious=πbestgen.Step 2:πlocal_besti=SP_FindBestNinsertnarrow2πbestgen,i.Step 3: ifLπlocal_besti<Lπbestgen, then πbestgen=πlocal_besti.Step 4: ifi<n−1, then i=i+1 and go to step 2.Step 5: ifi=n−1, then πbestcurrent=πbestgen.Ifπbestcurrent<πbestprevious, go to step 1.Step 6: seti=1, and πbestprevious=πbestgen.Step 7:πlocal_besti=SP_FindBestNswapnarrow2πbestgen,i.Step 8: ifLπlocal_besti<Lπbestgen, then πbestgen=πlocal_besti.Step 9: ifi<n−1, then i=i+1 and go to step 7.Step 10: ifi=n−1, then πbestcurrent=πbestgen.Ifπbestcurrent<πbestprevious, go to step 6.Step 11: setk=1 and πbestprevious=πbestgen.Step 12:πlocal_bestk=SP_FindBestNelim_intersectionπbestgen,k.Step 13: ifLπlocal_bestk<Lπbestgen, then πbestgen=πlocal_bestk.Step 14: ifk<K, then k=k+1 and go to step 12.Step 15: ifk=K, then πbestcurrent=πbestgen.Ifπbestcurrent<πbestprevious, go to step 11.Step 16: outputπbestgen.


#### 3.5.5. First Move Strategy

In order to reach more different regions, a first move strategy is adopted inSP_FindBestNinsertnarrow2π,i, SP_FindBestNswapnarrow2π,i, and SP_FindBestNelim_intersectionπ,k. Here is an example to illustrate the strategy. While the first neighbor πfirsti in Ninsertnarrow2π,i which can improve πbestgen is gained, SP_FindBestNinsertnarrow2π,i terminates and outputs πfirsti. Furthermore, the sequence π in SP_FindBestNinsertnarrow2π,i is replaced by πbestgen for the purpose of implementing the search instantly from the promising areas discovered by SP_FindBestNinsertnarrow2π,i. Obviously, the algorithm will search in the vicinity of varied suboptimal solutions, but not limited to the vicinity of the optimal solution, thus expanding the scope of the search area. The changed procedures with the above strategy, named SP_FindFirstMoveNinsertnarrow2π,i, SP_FindFirstMoveNswapnarrow2π,i, and SP_FindFirstMoveNelim_intersectionπ,k, are only different from SP_FindBestNinsertnarrow2π,i, SP_FindBestNswapnarrow2π,i, and SP_FindBestNelim_intersectionπ,k in step 2, step 7, and step 12. The pseudocodes of the changed steps are as follows:Step 2:πlocal_besti=SP_FindFirstMoveNinsertnarrow2πbestgen,i.Step 7:πlocal_besti=SP_FindFirstMoveNswapnarrow2πbestgen,i.Step 12:πlocal_bestk=SP_FindFirstMoveNelim_intersectionπbestgen,k.


#### 3.5.6. Procedure of Local Search

LetSPVND_FindFirstMoveNinsertnarrow2π, SPVND_FindFirstMoveNswapnarrow2π, and SPVND_FindFirstMoveNelim_intersectionπ denote SP_FindFirstMoveNinsertnarrow2π, SP_FindFirstMoveNswapnarrow2π, and SP_FindFirstMoveNelim_intersectionπ with the pipe VND method. Moreover, the basic VND described in Section 3.5.1 is also adopted. Then, the procedure of the IHACO’s local search is given in the following:Step 1. Perturbation phase.Setπ0=πbestgen;Forkk=1 to KPRandomly selecti and ji≠j;π1=insertπ0,i,j;π0=π1;End.Step 2: exploitation phase.Setloop=0;Repeatπ2=SPVND_FindFirstMoveNinsertnarrow2π1;π3=SPVND_FindFirstMoveNswapnarrow2π2;Ifπ3<π2 then π1=π3 andgo toπ2=SPVND_FindFirstMoveNinsertnarrow2π1;elseπ4=SPVND_FindFirstMoveNelim_intersectionπ3;Ifπ4<π3 then π1=π4 andgo toπ2=SPVND_FindFirstMoveNinsertnarrow2π1;elseloop++;end;end;untilloop=1.Step 3: ifLπ4<Lπbestgen, then πbestgen=π4.Step 4: outputπbestgen.In the above procedure, step 1 is a perturbation phase, which is helpful to prevent the cycling search and overcome the local optima. The special devised step 2 performs exploitation from the region obtained by step 1.


### 3.6. IHACO

Denotegen as a generation and maxgen as the maximum generation. Based on the above subsections, the proposed IHACO is given in the following:Step 1: parameter initialization. Set the crucial parameterm, maxgen, α, β, ρ, ρ0, S, and KP.Step 2: probabilistic model initialization. Setgen=1//Section 3.2.2Step 3: ant path construction. Generate the population at generationgen by sampling σkj and then evaluate its individuals to obtain the global best individual πbestgen//Section 3.3Step 4: local exploitation. Execute three types of local search toπbestgen and update πbestgen//Section 3.5Step 5: update mechanism of probability matrices. Calculate the objective function valueLπbestgen corresponding to the best individual πbestgen. Update πij and skj based on Lπbestgen and πbestgen//Section 3.4Step 6: setgen=gen+1. If gen≤maxgen, go to step 3.Step 7: output the best solution and its objective value found so far.As it can be seen from the above steps, not only does IHACO takes advantage of the new probabilistic model-based global search to execute exploration to obtain promising regions within solution space but also it makes use of the speed-up VND-based local search with the first move strategy to perform exploitation from these promising regions. On account of both global and local search are well balanced, the IHACO is expected to acquire good results.


## 3.1. Solution Representation

Usually, the encoding rules of ACO for solving VRP are that each ant represents one vehicle, the path of each ant represents the driving path of each vehicle, the depot is represented by the number 0, and the customer is represented by a natural number [32]. In this paper, the decimal coding based on customer arrangement is still used, but the path of each ant which is denoted πant represents the path of all customers served by the depot, and the depot is a separator for different vehicles. For example, the path sequence of the ant i is πanti=0,1,5,3,7,0,6,2,8,4,0 which means that vehicle 1 leaves the depot 0, visits customers 1, 5, 3, and 7 in turn, and then returns to the depot 0; vehicle 2 leaves the depot 0, visits customers 6, 2, 8, and 4 in turn, and finally returns to the depot 0.


## 3.2. Probabilistic Model and Its Initialization

The performance of ACO has strong connections with the probabilistic model. A better choice of the model is critical to improve the efficiency of the algorithm. In this research, two probability matrices are employed to form the probabilistic model. Besides, a few good individuals produced by a heuristic method are used to initialize these two probability matrices. This initialization method makes the algorithm converge to the high-quality solutions faster.


### 3.2.1. Probabilistic Model

In the traditional ant colony algorithm, the probabilistic model usually refers to ants’ pheromone matrix which is obtained by calculating the pheromone concentrationτij. This pheromone value represents the probability that the customer j appears immediately after the customer i. In other words, the pheromone matrix only reflects similar blocks of customers presented in the chosen individuals. But it is not enough to accumulate this kind of information alone. Counting the chance of customer j appearing before or in the current position is also important. Hence, a new probability matrix based on the order of customers called sequence information matrix is established. Then, the probabilistic model in this paper contains a pheromone matrix and a sequence information matrix.In a subpath visited by one vehicle, letskj be the appearance times that customer j before or in position k in the chosen sequences. Next, the specific probability value in the sequence information matrix is determined by(8)μkj=skj∑l∈Nksk,l,where Nk is the set of all available customers not already served by the current vehicle until position k. Here is an example to compute μkj. Suppose that the path taken by three ants are πant1=0,1,2,3,0,4,0, πant2=0,2,3,4,0,1,0, πant3=0,3,4,0,1,2,0, and πant4=0,3,2,1,4,0, and then, skj is given as follows:(9)skj5×5=−−−−−−3121−3432−4443−4444.Since the first location is fixed as the depot, the corresponding probability values ofμkj for the second position of each customer in N1=1,2,3,4 are μ21=3/3+1+2+1=0.43, μ22=1/3+1+2+1=0.14, μ23=2/3+1+2+1=0.29, and μ24=1/3+1+2+1=0.14.Then, the pheromone matrix is normalized by formula (10). Finally, let σkj denote the probability of customer j at position k which is calculated by formula (11). Based on this value, the ant selects the next customer to visit:(10)ζij=τij∑l∈Nkτi,l,(11)σkj=μkj+ζij,i=2,3,…,n.


### 3.2.2. Initialization of the Probabilistic Model Based on a Heuristic Rule

As mentioned in the introduction, normally, the probabilistic model (that is the pheromone matrix in ACO) is initialized to the same value at the beginning of the algorithm. In order to achieve good solutions more quickly, this paper adopts the initial solutions to initialize the proposed probabilistic model but not the same value.Sort the distance between all customers and the depot from small to large, and take the topt=roundn/10 customers into set cfirst. Let cfirst=c1first,c2first,…,ctfirst denote the first customer to be visited by the t-th ant, R=R1,R2,…,Rt the complete path passed by the t-th ant (also the t-th initial solution), and L=L1,L2,…,Lt the length of the route taken by the t-th ant. The procedure of the probabilistic model initialization is given as follows:Step 1: sett=1 and τij0=zerosn×n and skj0=zerosn×n.Step 2: the ant starts from the depot0 and visits the first customer ctfirst. Then, the ant selects the remaining customers for access based on the nearest neighbor principle. This process forms the t-th initial solution Rt.Step 3: calculate the route lengthLt corresponding to Rt.Step 4: setskjt=skjt−1+skjRt.Ift=1, then set τij1=1/L1; otherwise, setτijt=τijt−1+1/Lt1.Step 5: sett=t+1.Step 6: ift≤roundn/10, then go to step 2.Step 7: outputskjt and τijt.


## 3.2.1. Probabilistic Model

In the traditional ant colony algorithm, the probabilistic model usually refers to ants’ pheromone matrix which is obtained by calculating the pheromone concentrationτij. This pheromone value represents the probability that the customer j appears immediately after the customer i. In other words, the pheromone matrix only reflects similar blocks of customers presented in the chosen individuals. But it is not enough to accumulate this kind of information alone. Counting the chance of customer j appearing before or in the current position is also important. Hence, a new probability matrix based on the order of customers called sequence information matrix is established. Then, the probabilistic model in this paper contains a pheromone matrix and a sequence information matrix.In a subpath visited by one vehicle, letskj be the appearance times that customer j before or in position k in the chosen sequences. Next, the specific probability value in the sequence information matrix is determined by(8)μkj=skj∑l∈Nksk,l,where Nk is the set of all available customers not already served by the current vehicle until position k. Here is an example to compute μkj. Suppose that the path taken by three ants are πant1=0,1,2,3,0,4,0, πant2=0,2,3,4,0,1,0, πant3=0,3,4,0,1,2,0, and πant4=0,3,2,1,4,0, and then, skj is given as follows:(9)skj5×5=−−−−−−3121−3432−4443−4444.Since the first location is fixed as the depot, the corresponding probability values ofμkj for the second position of each customer in N1=1,2,3,4 are μ21=3/3+1+2+1=0.43, μ22=1/3+1+2+1=0.14, μ23=2/3+1+2+1=0.29, and μ24=1/3+1+2+1=0.14.Then, the pheromone matrix is normalized by formula (10). Finally, let σkj denote the probability of customer j at position k which is calculated by formula (11). Based on this value, the ant selects the next customer to visit:(10)ζij=τij∑l∈Nkτi,l,(11)σkj=μkj+ζij,i=2,3,…,n.


## 3.2.2. Initialization of the Probabilistic Model Based on a Heuristic Rule

As mentioned in the introduction, normally, the probabilistic model (that is the pheromone matrix in ACO) is initialized to the same value at the beginning of the algorithm. In order to achieve good solutions more quickly, this paper adopts the initial solutions to initialize the proposed probabilistic model but not the same value.Sort the distance between all customers and the depot from small to large, and take the topt=roundn/10 customers into set cfirst. Let cfirst=c1first,c2first,…,ctfirst denote the first customer to be visited by the t-th ant, R=R1,R2,…,Rt the complete path passed by the t-th ant (also the t-th initial solution), and L=L1,L2,…,Lt the length of the route taken by the t-th ant. The procedure of the probabilistic model initialization is given as follows:Step 1: sett=1 and τij0=zerosn×n and skj0=zerosn×n.Step 2: the ant starts from the depot0 and visits the first customer ctfirst. Then, the ant selects the remaining customers for access based on the nearest neighbor principle. This process forms the t-th initial solution Rt.Step 3: calculate the route lengthLt corresponding to Rt.Step 4: setskjt=skjt−1+skjRt.Ift=1, then set τij1=1/L1; otherwise, setτijt=τijt−1+1/Lt1.Step 5: sett=t+1.Step 6: ift≤roundn/10, then go to step 2.Step 7: outputskjt and τijt.


## 3.3. Ant Path Construction

In this study, each ant produces a complete solution for the MCVRP. The ant begins from the depot and continuously chooses customers to its tour as long as the constraint of capacity is satisfied. If there are no feasible customers to select, the ant comes back to the depot and restarts another trip.Letm denote the amounts of ants and t the amounts of the first customers selected. Set t=0,1,…,m. To guarantee the diversity of individuals, the first customer is picked randomly from n−t customers. For example, the first ant randomly selects one of the n customers as the first one to visit. Then, the first customer to be visited by the second ant will be randomly selected from the remaining n−1 customers. If m>n, when t=n+1, reset t to zero.Afterwards, each ant proceeds to add customers to its route based on the probability function (12). The probability of choosing customer j as the next one in the t-th iteration is given by formula (13):(12)δij=σkjαηijβ.In formula (12), σkj contains information about how often the customer j appears in the current position k, ηij is the inverse of distance between customer i and j, and α and β denote the importance of σkj and ηij, respectively:(13)pijt=δij∑l∈Nkδil,ifj∈Nk,0,otherwise.Setq0 be a constant value and q a random variable which is uniformly distributed in [0, 1]. If q>q0, select the next customer in accordance with the formulation of pijt; otherwise, select the customer who has the maximum value of δij [22].


## 3.4. Update Mechanism of Probability Matrices

The probabilistic model consists of two probability matricesskj and τij. They are updated by the best solution found at present. Let ρ0≤ρ<1 be the pheromone persistence, then equation (14) is used to update the pheromone concentration in the path. Meanwhile, equation (15) is adopted to update the probability matrix based on the order of customers:(14)τijt+1=ρ×τijt+1/Lbestt,if arci,j∈Lbestt,ρ×τijt,otherwise.In equation (14), Lbestt is the total length corresponding to the optimal ant path in the t-th iteration:(15)skjt+1=skjt+skjπbestt.In equation (15), skjπbestt is the probability matrix skj obtained by the t-th optimal individual πbestt. We refer to Section 3.2.1 for more details.


## 3.5. Speed-Up VND-Based Local Search

Generally, for the combinatorial optimization problems, e.g., MCVRP, there are many other high-quality solutions near the high-quality solutions found so far. The algorithm in this section performs the problem-dependent local search on the global optimal solutionπbestgen, for the purpose of realizing the detailed search of the area near πbestgen. In this paper, the operations of insert, swap, and intersection elimination based on VND, the speed-up search strategy, and the first move strategy are adopted.


### 3.5.1. Basic and Pipe VND

The VND is a variant of the variable neighborhood search. The basic VND [33] explores the list of neighborhoods in a sequential way one by one. Once an improvement is detected in any neighborhood structure, the search goes back to the ﬁrst neighborhood from the list. Until all the neighborhoods encounter the local optimum, the search ends. However, the pipe VND [33] executes continuous search in the same neighborhood if an improvement is obtained. In this paper, both these VNDs are applied.


### 3.5.2. Operation of Insert and Swap

Assume an ant’s tour is represented asπantold=1,…,i−1,i,i+1,…,j−1,j,j+1,…,n. The insert operator selects one customer like i and inserts it before or after another customer like j. Then, the tour becomes πantnew_insert=1,…,i−1,i+1,…,j−1,i,j,j+1,…,n or πantnew_insert=1,…,i−1,i+1,…,j−1,j,i,j+1,…,n. The swap operator selects two customers like i and j on this tour and swaps them. Then, πantnew_swap=1,…,i−1,j,i+1,…,j−1,i,j+1,…n.


### 3.5.3. Operation of Intersection Elimination

Since the objective function of the MCVRP is constructed to minimize the total length of the path, the optimal path of each vehicle should not be crossed.Theorem 1.
In the route of every vehicle, if two lines intersect, the total path length of the vehicle is larger than that of the total length without intersection.Proof.
As shown in Figure1, there is one vehicle route denoted by πant_vehicleold=0,1,3,2,4,0. For any two given intersecting line segments, such as lines L13 and L24 which intersect at point u, the end points of the two segments that can form the opposite side are connected like L12 and L34. Next, the new path without intersection is obtained and denoted by πant_vehiclenew=0,1,2,3,4,0. According to the theorem that the sum of any two sides in a triangle is greater than the third side, in Δ1u2, L1u+L2u>L12; in Δ3u4, Lu3+Lu4>L34. Let Ltotalold be the total length of vehicle path with intersecting lines, then Ltotalold=L01+L1u+Lu3+L32+L2u+Lu4+L40. Let Ltotalnew be the total length after eliminating intersecting lines, then Ltotalnew=L01+L12+L23+L34+L40. As a result, Ltotalnew−Ltotalold=L12−L1u+L2u+L34−Lu3+Lu4<0.
Fromπant_vehicleold and πant_vehiclenew defined above, it can be seen that if two paths intersect, there are four customers involved. According to the order of customers visited by the vehicle, only exchange positions of the second and third customers and then cross paths can be eliminated. In this paper, the geometric method is used to determine whether the crossing occurs.
Suppose the coordinates of customers 1, 3, 2, and 4 in Figure1 are x1,y1, x2,y2, x3,y3, and x4,y4, and then, the linear equations of line segments L13 and L24 are shown in the following equations:(16)x=x1+λx2−x1,y=y1+λy2−y1,(17)x=x3+ψx4−x3,y=y3+ψy4−y3,(18)x2−x1λ+x3−x4ψ=x3−x1,y2−y1λ+y3−y4ψ=y3−y1.
In formulas (16) and (17), λ is the slope of the line where L13 is located, and ψ is the slope of the line where L24 is located.
IfL13 and L24 intersect, then formula (18) has a solution. In other words, inequality (19) should be satisfied. Meanwhile, it is necessary to satisfy the inequality (22) in order to make two line segments intersect within the line segment (this is the case of path intersection shown in Figure 1), rather than the vertex of the line segment or its extension line. However, equations (20) and (21) are used to calculate λ and ψ, respectively:(19)Δ=x2−x1x3−x4y2−y1y3−y4≠0,(20)λ=1Δ×x3−x1x3−x4y3−y1y3−y4,(21)ψ=1Δ×x2−x1x3−x1y2−y1y3−y1,(22)0<λ<1,0<ψ<1.Figure 1
Route of one vehicle.


### 3.5.4. Speed-Up Search Strategy

Letinsertπ,i,j, swapπ,i,j, and elim_intersectionπ,i,j represent the operations of inserting, swapping, and uncrossing the i-th and j-th elements in solution π, respectively. Let IJinsertinfeasible denote the set of i,j whose neighborhood is infeasible (the infeasibility here means that the vehicle capacity constraint is not satisfied) after executing insertπ,i,j, and IJswapinfeasible is the set of i,j whose neighborhood is infeasible after executing swapπ,i,j. Because the intersection elimination operation is only carried out within the customers served by each vehicle, there is no neighborhood that violates the capacity constraint. Let Ninsertπ,i denote the feasible neighborhood set based on insert operation for i customers in solution π, Nswapπ,i the feasible neighborhood set based on swap operation for i customers in solution π, and Nelim_intersectionπ,k the feasible neighborhood set based on intersection elimination operation for customers serviced by the k-th vehicle in solution π. Set nk to be the number of customers visited by the k-th vehicle, then l=∑i=1kni−nk+k+1 is the position of the first customer visited by the k-th vehicle in solution π. Hence, there are as follows:(23)Ninsertπ,i=πi,j=insertπ,i,j|i=1,2,…,n−1,j=i+1,…,n,andi,j∉IJinsertinfeasible,(24)Nswapπ,i=πi,j=swapπ,i,j|i=1,2,…,n−1,j=i+1,…,n,andi,j∉IJswapinfeasible,(25)Nelim_intersectionπ,k=πi,j=elim_intersectionπ,i,j|i=l,j=i+1,…,i+nk−1.It can be seen from (23)–(25) that Ninsertπ,i is a feasible neighborhood set that any customer i inserts before or after any possible customer j within the same vehicle or between different vehicles. Nswapπ,i is a feasible neighborhood set that any customer i swaps with any possible customer j within the same vehicle or between different vehicles. Nelim_intersectionπ,k is a feasible neighborhood set that customers involved in the cross-line can be interchanged within the k-th vehicle.For the MCVRP, if customerj is too far away from customer i, the total route distance after the insert and swap operation will be greater than before, which makes the two operations meaningless. Hence, in order to avoid the search of invalid areas, the distances between all possible customers j and customer i are arranged from small to large, and just, the first S points are taken out to insert and swap. Thus, the more compact neighborhood sets Ninsertnarrow1π,i and Nswapnarrow1π,i can be obtained:(26)Ninsertnarrow1π,i=πi,j=insertπ,i,j|i=1,2,…,n−1,j∈S,andi,j∉IJinsertinfeasible,Nswapnarrow1π,i=πi,j=swapπ,i,j|i=1,2,…,n−1,j∈S,andi,j∉IJswapinfeasible.For the sake of speeding up the searching process and making the neighborhood search focus on high-quality areas, the nature of the problem can be utilized to further decrease invalid operations and reduce the feasible neighborhood set. LetLπ be the value of the objective function which is the total mileage of the route corresponding to solution π and dπi,πj be the distance between element πi which is located in position i of π and element πj which is located in position j of π.ForNinsertnarrow1π,i, if customeri is inserted before customer j, Lπi,j−Lπ=dπi,j−dπ=dπi−1,πi+1+dπj−1,πi+dπi,πj−dπi−1,πi+dπi,πi+1+dπj−1,πj; if customer i is inserted after customer j, Lπi,j−Lπ=dπi,j−dπ=dπi−1,πi+1+dπj,πi+dπi,πj+1−dπi−1,πi+dπiπi+1+dπj,πj+1. For any i,j, if dπi,j−dπ>0, then Lπi,j>Lπ. i,j in this case can be directly excluded without calculating Lπi,j. If Lπi,j<Lπ, then Lπi,j can be obtained directly from Lπi,j=Lπ+dπi,j−dπ instead of starting from the first customer of π to calculate. Therefore, a more compact and effective neighborhood set Ninsertnarrow2π,i can be obtained on the basis of Ninsertnarrow1π,i:(27)Ninsertnarrow2π,i=πi,j=insertπ,i,j|i=1,2,…,n−1,j∈S,i,j∉IJinsertinfeasible,anddπi,j−dπ<0.ForNswapnarrow1π,i, Lπi,j−Lπ=dπi,j−dπ=dπi−1,πj+dπj,πi+1+dπj−1,πi+dπi,πj+1−dπi−1,πi+dπi,πi+1+dπj−1,πj+dπj,πj+1. And for any i,j, if dπi,j−dπ>0, then Lπi,j>Lπ. i,j in this case can be directly excluded without calculating Lπi,j. If Lπi,j<Lπ, then Lπi,j can be obtained directly from Lπi,j=Lπ+dπi,j−dπ instead of starting from the first customer of π to calculate. Therefore, a more compact and effective neighborhood set Nswapnarrow2π,i can be obtained on the basis of Nswapnarrow1π,i:(28)Nswapnarrow2π,i=πi,j=swapπ,i,j|i=1,2,…,n−1,j∈S,i,j∉IJswapinfeasible,anddπi,j−dπ<0.ForNelim_intersectionπ,k, we can use the method described in Section 3.5.3 to determine whether two paths intersect. If so, according to Theorem 1, we know Lπi,j−Lπ=dπi,j−dπ=dπi−1,πj+dπi,πj+1−dπi−1,πi+dπj,πj+1<0. Then, Lπi,j can be obtained directly from Lπi,j=Lπ+dπi,j−dπinstead of starting from the first customer of π to calculate.LetSP_FindBestNinsertnarrow2π,i, SP_FindBestNswapnarrow2π,i, and SP_FindBestNelim_intersectionπ,k be the optimal neighborhoods in searching and outputting Ninsertnarrow2π,i, Nswapnarrow2π,i, and Nelim_intersectionπ,k, respectively. These three optimal neighborhoods utilize the above speed-up search method in the local search of IHACO. Meanwhile, the pipe VND described in Section 3.5.1 is also applied. The corresponding procedures are given as follows:Step 1: seti=1 and πbestprevious=πbestgen.Step 2:πlocal_besti=SP_FindBestNinsertnarrow2πbestgen,i.Step 3: ifLπlocal_besti<Lπbestgen, then πbestgen=πlocal_besti.Step 4: ifi<n−1, then i=i+1 and go to step 2.Step 5: ifi=n−1, then πbestcurrent=πbestgen.Ifπbestcurrent<πbestprevious, go to step 1.Step 6: seti=1, and πbestprevious=πbestgen.Step 7:πlocal_besti=SP_FindBestNswapnarrow2πbestgen,i.Step 8: ifLπlocal_besti<Lπbestgen, then πbestgen=πlocal_besti.Step 9: ifi<n−1, then i=i+1 and go to step 7.Step 10: ifi=n−1, then πbestcurrent=πbestgen.Ifπbestcurrent<πbestprevious, go to step 6.Step 11: setk=1 and πbestprevious=πbestgen.Step 12:πlocal_bestk=SP_FindBestNelim_intersectionπbestgen,k.Step 13: ifLπlocal_bestk<Lπbestgen, then πbestgen=πlocal_bestk.Step 14: ifk<K, then k=k+1 and go to step 12.Step 15: ifk=K, then πbestcurrent=πbestgen.Ifπbestcurrent<πbestprevious, go to step 11.Step 16: outputπbestgen.


### 3.5.5. First Move Strategy

In order to reach more different regions, a first move strategy is adopted inSP_FindBestNinsertnarrow2π,i, SP_FindBestNswapnarrow2π,i, and SP_FindBestNelim_intersectionπ,k. Here is an example to illustrate the strategy. While the first neighbor πfirsti in Ninsertnarrow2π,i which can improve πbestgen is gained, SP_FindBestNinsertnarrow2π,i terminates and outputs πfirsti. Furthermore, the sequence π in SP_FindBestNinsertnarrow2π,i is replaced by πbestgen for the purpose of implementing the search instantly from the promising areas discovered by SP_FindBestNinsertnarrow2π,i. Obviously, the algorithm will search in the vicinity of varied suboptimal solutions, but not limited to the vicinity of the optimal solution, thus expanding the scope of the search area. The changed procedures with the above strategy, named SP_FindFirstMoveNinsertnarrow2π,i, SP_FindFirstMoveNswapnarrow2π,i, and SP_FindFirstMoveNelim_intersectionπ,k, are only different from SP_FindBestNinsertnarrow2π,i, SP_FindBestNswapnarrow2π,i, and SP_FindBestNelim_intersectionπ,k in step 2, step 7, and step 12. The pseudocodes of the changed steps are as follows:Step 2:πlocal_besti=SP_FindFirstMoveNinsertnarrow2πbestgen,i.Step 7:πlocal_besti=SP_FindFirstMoveNswapnarrow2πbestgen,i.Step 12:πlocal_bestk=SP_FindFirstMoveNelim_intersectionπbestgen,k.


### 3.5.6. Procedure of Local Search

LetSPVND_FindFirstMoveNinsertnarrow2π, SPVND_FindFirstMoveNswapnarrow2π, and SPVND_FindFirstMoveNelim_intersectionπ denote SP_FindFirstMoveNinsertnarrow2π, SP_FindFirstMoveNswapnarrow2π, and SP_FindFirstMoveNelim_intersectionπ with the pipe VND method. Moreover, the basic VND described in Section 3.5.1 is also adopted. Then, the procedure of the IHACO’s local search is given in the following:Step 1. Perturbation phase.Setπ0=πbestgen;Forkk=1 to KPRandomly selecti and ji≠j;π1=insertπ0,i,j;π0=π1;End.Step 2: exploitation phase.Setloop=0;Repeatπ2=SPVND_FindFirstMoveNinsertnarrow2π1;π3=SPVND_FindFirstMoveNswapnarrow2π2;Ifπ3<π2 then π1=π3 andgo toπ2=SPVND_FindFirstMoveNinsertnarrow2π1;elseπ4=SPVND_FindFirstMoveNelim_intersectionπ3;Ifπ4<π3 then π1=π4 andgo toπ2=SPVND_FindFirstMoveNinsertnarrow2π1;elseloop++;end;end;untilloop=1.Step 3: ifLπ4<Lπbestgen, then πbestgen=π4.Step 4: outputπbestgen.In the above procedure, step 1 is a perturbation phase, which is helpful to prevent the cycling search and overcome the local optima. The special devised step 2 performs exploitation from the region obtained by step 1.


## 3.5.1. Basic and Pipe VND

The VND is a variant of the variable neighborhood search. The basic VND [33] explores the list of neighborhoods in a sequential way one by one. Once an improvement is detected in any neighborhood structure, the search goes back to the ﬁrst neighborhood from the list. Until all the neighborhoods encounter the local optimum, the search ends. However, the pipe VND [33] executes continuous search in the same neighborhood if an improvement is obtained. In this paper, both these VNDs are applied.


## 3.5.2. Operation of Insert and Swap

Assume an ant’s tour is represented asπantold=1,…,i−1,i,i+1,…,j−1,j,j+1,…,n. The insert operator selects one customer like i and inserts it before or after another customer like j. Then, the tour becomes πantnew_insert=1,…,i−1,i+1,…,j−1,i,j,j+1,…,n or πantnew_insert=1,…,i−1,i+1,…,j−1,j,i,j+1,…,n. The swap operator selects two customers like i and j on this tour and swaps them. Then, πantnew_swap=1,…,i−1,j,i+1,…,j−1,i,j+1,…n.


## 3.5.3. Operation of Intersection Elimination

Since the objective function of the MCVRP is constructed to minimize the total length of the path, the optimal path of each vehicle should not be crossed.Theorem 1.
In the route of every vehicle, if two lines intersect, the total path length of the vehicle is larger than that of the total length without intersection.Proof.
As shown in Figure1, there is one vehicle route denoted by πant_vehicleold=0,1,3,2,4,0. For any two given intersecting line segments, such as lines L13 and L24 which intersect at point u, the end points of the two segments that can form the opposite side are connected like L12 and L34. Next, the new path without intersection is obtained and denoted by πant_vehiclenew=0,1,2,3,4,0. According to the theorem that the sum of any two sides in a triangle is greater than the third side, in Δ1u2, L1u+L2u>L12; in Δ3u4, Lu3+Lu4>L34. Let Ltotalold be the total length of vehicle path with intersecting lines, then Ltotalold=L01+L1u+Lu3+L32+L2u+Lu4+L40. Let Ltotalnew be the total length after eliminating intersecting lines, then Ltotalnew=L01+L12+L23+L34+L40. As a result, Ltotalnew−Ltotalold=L12−L1u+L2u+L34−Lu3+Lu4<0.
Fromπant_vehicleold and πant_vehiclenew defined above, it can be seen that if two paths intersect, there are four customers involved. According to the order of customers visited by the vehicle, only exchange positions of the second and third customers and then cross paths can be eliminated. In this paper, the geometric method is used to determine whether the crossing occurs.
Suppose the coordinates of customers 1, 3, 2, and 4 in Figure1 are x1,y1, x2,y2, x3,y3, and x4,y4, and then, the linear equations of line segments L13 and L24 are shown in the following equations:(16)x=x1+λx2−x1,y=y1+λy2−y1,(17)x=x3+ψx4−x3,y=y3+ψy4−y3,(18)x2−x1λ+x3−x4ψ=x3−x1,y2−y1λ+y3−y4ψ=y3−y1.
In formulas (16) and (17), λ is the slope of the line where L13 is located, and ψ is the slope of the line where L24 is located.
IfL13 and L24 intersect, then formula (18) has a solution. In other words, inequality (19) should be satisfied. Meanwhile, it is necessary to satisfy the inequality (22) in order to make two line segments intersect within the line segment (this is the case of path intersection shown in Figure 1), rather than the vertex of the line segment or its extension line. However, equations (20) and (21) are used to calculate λ and ψ, respectively:(19)Δ=x2−x1x3−x4y2−y1y3−y4≠0,(20)λ=1Δ×x3−x1x3−x4y3−y1y3−y4,(21)ψ=1Δ×x2−x1x3−x1y2−y1y3−y1,(22)0<λ<1,0<ψ<1.Figure 1
Route of one vehicle.


## 3.5.4. Speed-Up Search Strategy

Letinsertπ,i,j, swapπ,i,j, and elim_intersectionπ,i,j represent the operations of inserting, swapping, and uncrossing the i-th and j-th elements in solution π, respectively. Let IJinsertinfeasible denote the set of i,j whose neighborhood is infeasible (the infeasibility here means that the vehicle capacity constraint is not satisfied) after executing insertπ,i,j, and IJswapinfeasible is the set of i,j whose neighborhood is infeasible after executing swapπ,i,j. Because the intersection elimination operation is only carried out within the customers served by each vehicle, there is no neighborhood that violates the capacity constraint. Let Ninsertπ,i denote the feasible neighborhood set based on insert operation for i customers in solution π, Nswapπ,i the feasible neighborhood set based on swap operation for i customers in solution π, and Nelim_intersectionπ,k the feasible neighborhood set based on intersection elimination operation for customers serviced by the k-th vehicle in solution π. Set nk to be the number of customers visited by the k-th vehicle, then l=∑i=1kni−nk+k+1 is the position of the first customer visited by the k-th vehicle in solution π. Hence, there are as follows:(23)Ninsertπ,i=πi,j=insertπ,i,j|i=1,2,…,n−1,j=i+1,…,n,andi,j∉IJinsertinfeasible,(24)Nswapπ,i=πi,j=swapπ,i,j|i=1,2,…,n−1,j=i+1,…,n,andi,j∉IJswapinfeasible,(25)Nelim_intersectionπ,k=πi,j=elim_intersectionπ,i,j|i=l,j=i+1,…,i+nk−1.It can be seen from (23)–(25) that Ninsertπ,i is a feasible neighborhood set that any customer i inserts before or after any possible customer j within the same vehicle or between different vehicles. Nswapπ,i is a feasible neighborhood set that any customer i swaps with any possible customer j within the same vehicle or between different vehicles. Nelim_intersectionπ,k is a feasible neighborhood set that customers involved in the cross-line can be interchanged within the k-th vehicle.For the MCVRP, if customerj is too far away from customer i, the total route distance after the insert and swap operation will be greater than before, which makes the two operations meaningless. Hence, in order to avoid the search of invalid areas, the distances between all possible customers j and customer i are arranged from small to large, and just, the first S points are taken out to insert and swap. Thus, the more compact neighborhood sets Ninsertnarrow1π,i and Nswapnarrow1π,i can be obtained:(26)Ninsertnarrow1π,i=πi,j=insertπ,i,j|i=1,2,…,n−1,j∈S,andi,j∉IJinsertinfeasible,Nswapnarrow1π,i=πi,j=swapπ,i,j|i=1,2,…,n−1,j∈S,andi,j∉IJswapinfeasible.For the sake of speeding up the searching process and making the neighborhood search focus on high-quality areas, the nature of the problem can be utilized to further decrease invalid operations and reduce the feasible neighborhood set. LetLπ be the value of the objective function which is the total mileage of the route corresponding to solution π and dπi,πj be the distance between element πi which is located in position i of π and element πj which is located in position j of π.ForNinsertnarrow1π,i, if customeri is inserted before customer j, Lπi,j−Lπ=dπi,j−dπ=dπi−1,πi+1+dπj−1,πi+dπi,πj−dπi−1,πi+dπi,πi+1+dπj−1,πj; if customer i is inserted after customer j, Lπi,j−Lπ=dπi,j−dπ=dπi−1,πi+1+dπj,πi+dπi,πj+1−dπi−1,πi+dπiπi+1+dπj,πj+1. For any i,j, if dπi,j−dπ>0, then Lπi,j>Lπ. i,j in this case can be directly excluded without calculating Lπi,j. If Lπi,j<Lπ, then Lπi,j can be obtained directly from Lπi,j=Lπ+dπi,j−dπ instead of starting from the first customer of π to calculate. Therefore, a more compact and effective neighborhood set Ninsertnarrow2π,i can be obtained on the basis of Ninsertnarrow1π,i:(27)Ninsertnarrow2π,i=πi,j=insertπ,i,j|i=1,2,…,n−1,j∈S,i,j∉IJinsertinfeasible,anddπi,j−dπ<0.ForNswapnarrow1π,i, Lπi,j−Lπ=dπi,j−dπ=dπi−1,πj+dπj,πi+1+dπj−1,πi+dπi,πj+1−dπi−1,πi+dπi,πi+1+dπj−1,πj+dπj,πj+1. And for any i,j, if dπi,j−dπ>0, then Lπi,j>Lπ. i,j in this case can be directly excluded without calculating Lπi,j. If Lπi,j<Lπ, then Lπi,j can be obtained directly from Lπi,j=Lπ+dπi,j−dπ instead of starting from the first customer of π to calculate. Therefore, a more compact and effective neighborhood set Nswapnarrow2π,i can be obtained on the basis of Nswapnarrow1π,i:(28)Nswapnarrow2π,i=πi,j=swapπ,i,j|i=1,2,…,n−1,j∈S,i,j∉IJswapinfeasible,anddπi,j−dπ<0.ForNelim_intersectionπ,k, we can use the method described in Section 3.5.3 to determine whether two paths intersect. If so, according to Theorem 1, we know Lπi,j−Lπ=dπi,j−dπ=dπi−1,πj+dπi,πj+1−dπi−1,πi+dπj,πj+1<0. Then, Lπi,j can be obtained directly from Lπi,j=Lπ+dπi,j−dπinstead of starting from the first customer of π to calculate.LetSP_FindBestNinsertnarrow2π,i, SP_FindBestNswapnarrow2π,i, and SP_FindBestNelim_intersectionπ,k be the optimal neighborhoods in searching and outputting Ninsertnarrow2π,i, Nswapnarrow2π,i, and Nelim_intersectionπ,k, respectively. These three optimal neighborhoods utilize the above speed-up search method in the local search of IHACO. Meanwhile, the pipe VND described in Section 3.5.1 is also applied. The corresponding procedures are given as follows:Step 1: seti=1 and πbestprevious=πbestgen.Step 2:πlocal_besti=SP_FindBestNinsertnarrow2πbestgen,i.Step 3: ifLπlocal_besti<Lπbestgen, then πbestgen=πlocal_besti.Step 4: ifi<n−1, then i=i+1 and go to step 2.Step 5: ifi=n−1, then πbestcurrent=πbestgen.Ifπbestcurrent<πbestprevious, go to step 1.Step 6: seti=1, and πbestprevious=πbestgen.Step 7:πlocal_besti=SP_FindBestNswapnarrow2πbestgen,i.Step 8: ifLπlocal_besti<Lπbestgen, then πbestgen=πlocal_besti.Step 9: ifi<n−1, then i=i+1 and go to step 7.Step 10: ifi=n−1, then πbestcurrent=πbestgen.Ifπbestcurrent<πbestprevious, go to step 6.Step 11: setk=1 and πbestprevious=πbestgen.Step 12:πlocal_bestk=SP_FindBestNelim_intersectionπbestgen,k.Step 13: ifLπlocal_bestk<Lπbestgen, then πbestgen=πlocal_bestk.Step 14: ifk<K, then k=k+1 and go to step 12.Step 15: ifk=K, then πbestcurrent=πbestgen.Ifπbestcurrent<πbestprevious, go to step 11.Step 16: outputπbestgen.


## 3.5.5. First Move Strategy

In order to reach more different regions, a first move strategy is adopted inSP_FindBestNinsertnarrow2π,i, SP_FindBestNswapnarrow2π,i, and SP_FindBestNelim_intersectionπ,k. Here is an example to illustrate the strategy. While the first neighbor πfirsti in Ninsertnarrow2π,i which can improve πbestgen is gained, SP_FindBestNinsertnarrow2π,i terminates and outputs πfirsti. Furthermore, the sequence π in SP_FindBestNinsertnarrow2π,i is replaced by πbestgen for the purpose of implementing the search instantly from the promising areas discovered by SP_FindBestNinsertnarrow2π,i. Obviously, the algorithm will search in the vicinity of varied suboptimal solutions, but not limited to the vicinity of the optimal solution, thus expanding the scope of the search area. The changed procedures with the above strategy, named SP_FindFirstMoveNinsertnarrow2π,i, SP_FindFirstMoveNswapnarrow2π,i, and SP_FindFirstMoveNelim_intersectionπ,k, are only different from SP_FindBestNinsertnarrow2π,i, SP_FindBestNswapnarrow2π,i, and SP_FindBestNelim_intersectionπ,k in step 2, step 7, and step 12. The pseudocodes of the changed steps are as follows:Step 2:πlocal_besti=SP_FindFirstMoveNinsertnarrow2πbestgen,i.Step 7:πlocal_besti=SP_FindFirstMoveNswapnarrow2πbestgen,i.Step 12:πlocal_bestk=SP_FindFirstMoveNelim_intersectionπbestgen,k.


## 3.5.6. Procedure of Local Search

LetSPVND_FindFirstMoveNinsertnarrow2π, SPVND_FindFirstMoveNswapnarrow2π, and SPVND_FindFirstMoveNelim_intersectionπ denote SP_FindFirstMoveNinsertnarrow2π, SP_FindFirstMoveNswapnarrow2π, and SP_FindFirstMoveNelim_intersectionπ with the pipe VND method. Moreover, the basic VND described in Section 3.5.1 is also adopted. Then, the procedure of the IHACO’s local search is given in the following:Step 1. Perturbation phase.Setπ0=πbestgen;Forkk=1 to KPRandomly selecti and ji≠j;π1=insertπ0,i,j;π0=π1;End.Step 2: exploitation phase.Setloop=0;Repeatπ2=SPVND_FindFirstMoveNinsertnarrow2π1;π3=SPVND_FindFirstMoveNswapnarrow2π2;Ifπ3<π2 then π1=π3 andgo toπ2=SPVND_FindFirstMoveNinsertnarrow2π1;elseπ4=SPVND_FindFirstMoveNelim_intersectionπ3;Ifπ4<π3 then π1=π4 andgo toπ2=SPVND_FindFirstMoveNinsertnarrow2π1;elseloop++;end;end;untilloop=1.Step 3: ifLπ4<Lπbestgen, then πbestgen=π4.Step 4: outputπbestgen.In the above procedure, step 1 is a perturbation phase, which is helpful to prevent the cycling search and overcome the local optima. The special devised step 2 performs exploitation from the region obtained by step 1.


## 3.6. IHACO

Denotegen as a generation and maxgen as the maximum generation. Based on the above subsections, the proposed IHACO is given in the following:Step 1: parameter initialization. Set the crucial parameterm, maxgen, α, β, ρ, ρ0, S, and KP.Step 2: probabilistic model initialization. Setgen=1//Section 3.2.2Step 3: ant path construction. Generate the population at generationgen by sampling σkj and then evaluate its individuals to obtain the global best individual πbestgen//Section 3.3Step 4: local exploitation. Execute three types of local search toπbestgen and update πbestgen//Section 3.5Step 5: update mechanism of probability matrices. Calculate the objective function valueLπbestgen corresponding to the best individual πbestgen. Update πij and skj based on Lπbestgen and πbestgen//Section 3.4Step 6: setgen=gen+1. If gen≤maxgen, go to step 3.Step 7: output the best solution and its objective value found so far.As it can be seen from the above steps, not only does IHACO takes advantage of the new probabilistic model-based global search to execute exploration to obtain promising regions within solution space but also it makes use of the speed-up VND-based local search with the first move strategy to perform exploitation from these promising regions. On account of both global and local search are well balanced, the IHACO is expected to acquire good results.


## 4. Computational Experiments and Statistical Analysis

In this section, three sets of numerical simulations are conducted with the aim of checking the effectiveness and efficiency of IHACO. The experiments implemented have been coded in Matlab 2019a and tested on an Intel 1.8 GHz PC with 12 GB RAM.


### 4.1. Experimental Setup

The generation method of instances is derived from [13]. The test set in this paper contains 14 benchmark instances with the capacity constraint, including vrpcn1-vrpcn5 and vrpcn11-vrpcn12 of sets a and b. These scales of problem range from 50 to 200 customers.The main parameters of IHACO are the same as those in [22]: the number of ants m=20, two exponential factors α=1 and β=2, the pheromone concentration ρ=0.9, and the constant value q0=0.9. There are other parameters needed to set: the upper bound of range in local search S=n/5 and the perturbation variable KP=n/10. For the sake of having a fair comparison, all the compared algorithms terminate when either of the following two criteria is reached: the first one is the same runtime limit, and the second one is less than or equal to the optimal value found in the literature. It should be noted that when setting the termination time of the proposed algorithm, two steps are adopted. We first find the running time tbest and the running environment (that is, the CPU frequencyfbest) corresponding to the best solution of the four comprised algorithms in each instance. Suppose that the CPU frequency of this paper is f, then the termination time of the proposed algorithm t=tbest×fbest/f. In addition, to make the results clearer, in Tables 2 and 3, the best, the second-best, and the third-best values in each row are marked out with the bold, the bold and underlined, and the italic and underlined fonts.


### 4.2. Comparison of IHACO and Other Algorithms

To test the performance of the algorithm in this paper, various simulations are carried out to compare the proposed IHACO with a classical ACO (ACS) [21], an effective hybridized ACO (HAC) [22], and two new state-of-the-art approaches named HAVNS and HABC [24] for the MCVRP. In solving MCVRP, ACS is the earliest used ant colony optimization algorithm and is proved with high performance for two-compartment test problems. HAC that is presented through the inspiration of ACS combines the ant colony optimization algorithm with local search schemes. After testing new generated benchmark problem instances, the results manifest that HAC takes less computational time to get higher quality solutions for all given problems except one problem (vrpnc12a) compared to ACS. HAVNS and HABC are also heuristic-based. According to the experiment results in [24], for the 14 examples to be tested in this paper, HAVNS obtains 12 best solutions and updates 6 of them on the basis of HAC. Meanwhile, HABC gains 12 best solutions and updates 8 of them based on HAC.The experimental results of the five algorithms are shown in Table2. The first column in Table 2 represents the type of the instance, and the second column is the quantity of customers contained in every instance. The meaning of Obj. and T(s) in this table is the objective function value and runtime, respectively.Table 2
Comparison results of IHACO and four algorithms.
Inst.nACSHACHAVNSHABCIHACOObj.t(s)Obj.T(s)Obj.T(s)Obj.t(s)Obj.T(s)vrpnc1a50569.56416550.705550.74550.425550.706vrpnc1b50569.11817551.945506.184502.836548.777vrpnc2a75957.52536890.6815881.8414890.6817873.5713vrpnc2b75954.85635918.9614918.961491417888.3815vrpnc3a100964.132122874.0740878.8321875.2427870.1135vrpnc3b100959.327122895.2644898.6420895.3925871.2216vrpnc4a1501253.863451126.121461126.121031126.121101110.0633vrpnc4b1501254.513361159.481511159.481201159.481421134.90116vrpnc5a1991587.026881444.292571408.862301385.713451423.11414vrpnc5b1991640.596761525.872361515.252181492.022121456.16125vrpnc11a1201133.882811110.45751110.45701110.45821110.4575vrpnc11b1201247.492801221.73871221.73771220.43951218.6671vrpnc12a100911.861105912.6415901.2213901.1535907.7342vrpnc12b100970.833100950.7930935.6222923.2543961.7451Average1069.612225.641009.5801000.9966.43996.2382.93994.6872.79As reported from Table2, the proposed IHACO finds the best solution for 9 instances of the 14 problem instances and updates 8 of them based on HAVNS and HABC. However, on the average value of all examples, IHACO performs best. Moreover, while the quantity of customers is more than 100, IHACO can hit 5 best solutions in 6 instances. Hence, as the problem size increases, IHACO is more effective. In addition, HAVNS and HABC were implemented on Intel Core™ i5-5200U @ 2.20 GHz 6,00 GO RAM, which was different from the IHACO operating environment. Considering the speed of CPU, IHACO is much faster than HAVNS and HABC through the conversion of 72.79/2.2/1.8=59.56 on average. Furthermore, as it can be seen from Table 2, the time consumption corresponding to the 9 good solutions obtained in this study is less than that of the best solution in three comparative literatures. Considering the difference of computing environment, if the conversion is carried out, the time consumption of these 9 solutions will be shorter. So, IHACO is also more efficient.


### 4.3. Comparison of IHACO with Its Six Variants

For the purpose of evaluating the crucial operations of IHACO, IHACO is also compared with its six variants, whose abbreviations are described in the following:(1)
IHACO_V1 : in IHACO’s global search,skj is not adopted, and only τij is used to realize the generation of new individuals.(2)
IHACO_V2 : in IHACO’s global search, when the probabilistic model initialization phase is carried out,τij is initialized to the same value 1/Lneig, and skj is initialized just by the individual πneig (the individual corresponding to the objective function value Lneig).(3)
IHACO_V3 : in IHACO’s local search,Ninsertnarrow1π,i and Nswapnarrow1π,i are not employed. In other words, all possible values for customer j (j∈S is replaced with j=i+1,…,n) will be assessed.(4)
IHACO_V4 : in IHACO’s local search, the first move strategy is not employed. This means that, in the operations of insert and swap, after all possible positions have been tried, and the customer chooses the location that has the highest improvement to insert or swap.(5)
IHACO_V5 : in IHACO’s local search, the pipe VND is not employed.(6)
IHACO_V6 : in IHACO’s local search, the basic VND is not employed.The statistical results are listed in Table3. Besides, the graphical results for the compared algorithms are exhibited in Figure 2.Table 3
Comparison results of IHACO and its variants.
Inst.nIHACO_V1IHACO_V2IHACO_V3IHACO_V4IHACO_V5IHACO_V6IHACOvrpnc1a50550.70550.70553.44552.19550.70550.70550.70vrpnc1b50553.48551.94557.08548.77548.77548.77548.77vrpnc2a75888.63896.90887.90896.01888.69895.42873.57vrpnc2b75891.71905.46920.60916.60900.64908.29888.38vrpnc3a100872.29876.28885.64894.72882.26875.06870.11vrpnc3b100879.28874.96893.21897.09881.33882.38871.22vrpnc4a1501115.581111.041143.921157.331125.191127.241110.06vrpnc4b1501151.641144.281167.451174.291151.071148.681134.90vrpnc5a1991444.351428.191494.801486.061424.341460.311423.11vrpnc5b1991462.211494.861531.331526.431472.541505.951456.16vrpnc11a1201128.121128.941172.761130.861123.601133.391110.45vrpnc11b1201226.901238.301243.821237.961231.041239.541218.66vrpnc12a100920.09922.78932.85939.67919.65919.81907.73vrpnc12b100975.66974.97978.48993.87967.61979.08961.74Average1004.331007.121025.951025.131004.821012.47994.68Figure 2
Comparison results of IHACO and its variants on different-scaled instances: (a) the customer scale is no more than 100; (b) the customer scale is greater than 100.
(a)(b)Compared with IHACO’s variants, it is clear from Table3 and Figure 2 that IHACO performs best for all instances. The results indicate the importance of adopting all the crucial operations of IHACO simultaneously. It can also be found from Figure 2 that, as the problem scale increases, the difference of the objective values obtained by these comparison algorithms gradually increases, and then, the superiority of IHACO becomes more obvious. This once again verified a conclusion of the previous experiment; that is, IHACO performs better when the number of customers is large.The above results can be further analysed. IHACO_V1 does not adopt the information sequence matrix, which may lead to the low accuracy of the information left by ants on the path and reduces the similarity between the sampled individuals and the optimal individuals. IHACO_V2 initializes the two matrixes that make up the probabilistic model to a unified value. As a result, the initial individuals have no information of excellent individuals to utilize, and then, the quality of the initial population is lowered. IHACO_V3 expands the range of invalid regions in operations of insert and swap. As a result, it could take a longer time to find high-quality solutions than IHACO. IHACO_V4 must find the optimal location to execute the insert and swap operations. This causes that the search areas are relatively concentrated, which may miss some promising regions. Both IHACO_V5 and IHACO_V6 reduce the depth of exploitation, so it is not easy for them to achieve high-quality solutions. Thus, the main operations removed from these six variants are all adopted by IHACO.


### 4.4. Superiority Verification of Using Multi-Compartment Vehicles

Since some different products have different characteristic properties, they cannot be mixed. That is to say, the products cannot be transported simultaneously by vehicles with single compartment. However, they can be transported by vehicles with single compartment for several times or by vehicles with multiple compartments at one time.In this section, for the two kinds of products involved in the problem instances adopted in this paper, single-compartment vehicles and two-compartment vehicles are both utilized to gather the needs of customers within the service area. When using single-compartment vehicles to collect goods from customers, the problem is decomposed into two subproblems. One is to seek the solution of VRP for collecting product 1, and another is for collecting product 2. The sum of the total distance for the two subproblems represents the final solution listed in column 3 of Table4. The results obtained by using two-compartment vehicles are the same as reported in column 5 of Table 2. The increase percentage of total length listed in column 5 of Table 4 can be obtained from the following equation:(29)increase%=SC−TCSC×100,where SC and TC are the total distance travelled by single-compartment vehicles and two-compartment vehicles, respectively. Besides, in order to make the comparison more intuitive, plot with 95% Tukey’s honest significant difference (HSD) confidence intervals for employing these two types of compartments under two categories of data sets is shown in Figure 3.Table 4
Comparison results of using different compartment vehicles.
Inst.nTotal length of using single compartmentTotal length of using two compartmentsTotal length increase percentage (%)vrpnc1a50932.13550.7069.26vrpnc1b50934.13548.7770.22vrpnc2a751317.96873.5750.87vrpnc2b751327.43888.3849.42vrpnc3a1001423.90870.1163.65vrpnc3b1001415.06871.2262.42vrpnc4a1501709.251110.0653.98vrpnc4b1501690.641134.9048.97vrpnc5a1992043.761423.1143.61vrpnc5b1992049.471456.1640.74vrpnc11a1201583.431110.4542.59vrpnc11b1201556.181218.6627.70vrpnc12a1001281.69907.7341.20vrpnc12b1001305.50961.7435.74Average1469.32994.6847.72Figure 3
Comparison results of vehicles with different compartments.It can be obviously seen from Table4 and Figure 3 that using single-compartment vehicles brings the notable increase in the total length of each route. The reason is that although the total route cost is reduced due to the increase of vehicle capacity while transporting the same product, each customer must be served twice in order to ensure that different products are all delivered. More specifically, it can be noticed from Table 4 that the total distance raises 47.72% on average when the vehicles with the single compartment are employed to visit the customers. Moreover, it is worth mentioning that the total distance could be reduced by up to 70.22%, while two-compartment vehicles are selected. The statistical results clearly display the profit of using multi-compartment vehicles.


## 4.1. Experimental Setup

The generation method of instances is derived from [13]. The test set in this paper contains 14 benchmark instances with the capacity constraint, including vrpcn1-vrpcn5 and vrpcn11-vrpcn12 of sets a and b. These scales of problem range from 50 to 200 customers.The main parameters of IHACO are the same as those in [22]: the number of ants m=20, two exponential factors α=1 and β=2, the pheromone concentration ρ=0.9, and the constant value q0=0.9. There are other parameters needed to set: the upper bound of range in local search S=n/5 and the perturbation variable KP=n/10. For the sake of having a fair comparison, all the compared algorithms terminate when either of the following two criteria is reached: the first one is the same runtime limit, and the second one is less than or equal to the optimal value found in the literature. It should be noted that when setting the termination time of the proposed algorithm, two steps are adopted. We first find the running time tbest and the running environment (that is, the CPU frequencyfbest) corresponding to the best solution of the four comprised algorithms in each instance. Suppose that the CPU frequency of this paper is f, then the termination time of the proposed algorithm t=tbest×fbest/f. In addition, to make the results clearer, in Tables 2 and 3, the best, the second-best, and the third-best values in each row are marked out with the bold, the bold and underlined, and the italic and underlined fonts.


## 4.2. Comparison of IHACO and Other Algorithms

To test the performance of the algorithm in this paper, various simulations are carried out to compare the proposed IHACO with a classical ACO (ACS) [21], an effective hybridized ACO (HAC) [22], and two new state-of-the-art approaches named HAVNS and HABC [24] for the MCVRP. In solving MCVRP, ACS is the earliest used ant colony optimization algorithm and is proved with high performance for two-compartment test problems. HAC that is presented through the inspiration of ACS combines the ant colony optimization algorithm with local search schemes. After testing new generated benchmark problem instances, the results manifest that HAC takes less computational time to get higher quality solutions for all given problems except one problem (vrpnc12a) compared to ACS. HAVNS and HABC are also heuristic-based. According to the experiment results in [24], for the 14 examples to be tested in this paper, HAVNS obtains 12 best solutions and updates 6 of them on the basis of HAC. Meanwhile, HABC gains 12 best solutions and updates 8 of them based on HAC.The experimental results of the five algorithms are shown in Table2. The first column in Table 2 represents the type of the instance, and the second column is the quantity of customers contained in every instance. The meaning of Obj. and T(s) in this table is the objective function value and runtime, respectively.Table 2
Comparison results of IHACO and four algorithms.
Inst.nACSHACHAVNSHABCIHACOObj.t(s)Obj.T(s)Obj.T(s)Obj.t(s)Obj.T(s)vrpnc1a50569.56416550.705550.74550.425550.706vrpnc1b50569.11817551.945506.184502.836548.777vrpnc2a75957.52536890.6815881.8414890.6817873.5713vrpnc2b75954.85635918.9614918.961491417888.3815vrpnc3a100964.132122874.0740878.8321875.2427870.1135vrpnc3b100959.327122895.2644898.6420895.3925871.2216vrpnc4a1501253.863451126.121461126.121031126.121101110.0633vrpnc4b1501254.513361159.481511159.481201159.481421134.90116vrpnc5a1991587.026881444.292571408.862301385.713451423.11414vrpnc5b1991640.596761525.872361515.252181492.022121456.16125vrpnc11a1201133.882811110.45751110.45701110.45821110.4575vrpnc11b1201247.492801221.73871221.73771220.43951218.6671vrpnc12a100911.861105912.6415901.2213901.1535907.7342vrpnc12b100970.833100950.7930935.6222923.2543961.7451Average1069.612225.641009.5801000.9966.43996.2382.93994.6872.79As reported from Table2, the proposed IHACO finds the best solution for 9 instances of the 14 problem instances and updates 8 of them based on HAVNS and HABC. However, on the average value of all examples, IHACO performs best. Moreover, while the quantity of customers is more than 100, IHACO can hit 5 best solutions in 6 instances. Hence, as the problem size increases, IHACO is more effective. In addition, HAVNS and HABC were implemented on Intel Core™ i5-5200U @ 2.20 GHz 6,00 GO RAM, which was different from the IHACO operating environment. Considering the speed of CPU, IHACO is much faster than HAVNS and HABC through the conversion of 72.79/2.2/1.8=59.56 on average. Furthermore, as it can be seen from Table 2, the time consumption corresponding to the 9 good solutions obtained in this study is less than that of the best solution in three comparative literatures. Considering the difference of computing environment, if the conversion is carried out, the time consumption of these 9 solutions will be shorter. So, IHACO is also more efficient.


## 4.3. Comparison of IHACO with Its Six Variants

For the purpose of evaluating the crucial operations of IHACO, IHACO is also compared with its six variants, whose abbreviations are described in the following:(1)
IHACO_V1 : in IHACO’s global search,skj is not adopted, and only τij is used to realize the generation of new individuals.(2)
IHACO_V2 : in IHACO’s global search, when the probabilistic model initialization phase is carried out,τij is initialized to the same value 1/Lneig, and skj is initialized just by the individual πneig (the individual corresponding to the objective function value Lneig).(3)
IHACO_V3 : in IHACO’s local search,Ninsertnarrow1π,i and Nswapnarrow1π,i are not employed. In other words, all possible values for customer j (j∈S is replaced with j=i+1,…,n) will be assessed.(4)
IHACO_V4 : in IHACO’s local search, the first move strategy is not employed. This means that, in the operations of insert and swap, after all possible positions have been tried, and the customer chooses the location that has the highest improvement to insert or swap.(5)
IHACO_V5 : in IHACO’s local search, the pipe VND is not employed.(6)
IHACO_V6 : in IHACO’s local search, the basic VND is not employed.The statistical results are listed in Table3. Besides, the graphical results for the compared algorithms are exhibited in Figure 2.Table 3
Comparison results of IHACO and its variants.
Inst.nIHACO_V1IHACO_V2IHACO_V3IHACO_V4IHACO_V5IHACO_V6IHACOvrpnc1a50550.70550.70553.44552.19550.70550.70550.70vrpnc1b50553.48551.94557.08548.77548.77548.77548.77vrpnc2a75888.63896.90887.90896.01888.69895.42873.57vrpnc2b75891.71905.46920.60916.60900.64908.29888.38vrpnc3a100872.29876.28885.64894.72882.26875.06870.11vrpnc3b100879.28874.96893.21897.09881.33882.38871.22vrpnc4a1501115.581111.041143.921157.331125.191127.241110.06vrpnc4b1501151.641144.281167.451174.291151.071148.681134.90vrpnc5a1991444.351428.191494.801486.061424.341460.311423.11vrpnc5b1991462.211494.861531.331526.431472.541505.951456.16vrpnc11a1201128.121128.941172.761130.861123.601133.391110.45vrpnc11b1201226.901238.301243.821237.961231.041239.541218.66vrpnc12a100920.09922.78932.85939.67919.65919.81907.73vrpnc12b100975.66974.97978.48993.87967.61979.08961.74Average1004.331007.121025.951025.131004.821012.47994.68Figure 2
Comparison results of IHACO and its variants on different-scaled instances: (a) the customer scale is no more than 100; (b) the customer scale is greater than 100.
(a)(b)Compared with IHACO’s variants, it is clear from Table3 and Figure 2 that IHACO performs best for all instances. The results indicate the importance of adopting all the crucial operations of IHACO simultaneously. It can also be found from Figure 2 that, as the problem scale increases, the difference of the objective values obtained by these comparison algorithms gradually increases, and then, the superiority of IHACO becomes more obvious. This once again verified a conclusion of the previous experiment; that is, IHACO performs better when the number of customers is large.The above results can be further analysed. IHACO_V1 does not adopt the information sequence matrix, which may lead to the low accuracy of the information left by ants on the path and reduces the similarity between the sampled individuals and the optimal individuals. IHACO_V2 initializes the two matrixes that make up the probabilistic model to a unified value. As a result, the initial individuals have no information of excellent individuals to utilize, and then, the quality of the initial population is lowered. IHACO_V3 expands the range of invalid regions in operations of insert and swap. As a result, it could take a longer time to find high-quality solutions than IHACO. IHACO_V4 must find the optimal location to execute the insert and swap operations. This causes that the search areas are relatively concentrated, which may miss some promising regions. Both IHACO_V5 and IHACO_V6 reduce the depth of exploitation, so it is not easy for them to achieve high-quality solutions. Thus, the main operations removed from these six variants are all adopted by IHACO.


## 4.4. Superiority Verification of Using Multi-Compartment Vehicles

Since some different products have different characteristic properties, they cannot be mixed. That is to say, the products cannot be transported simultaneously by vehicles with single compartment. However, they can be transported by vehicles with single compartment for several times or by vehicles with multiple compartments at one time.In this section, for the two kinds of products involved in the problem instances adopted in this paper, single-compartment vehicles and two-compartment vehicles are both utilized to gather the needs of customers within the service area. When using single-compartment vehicles to collect goods from customers, the problem is decomposed into two subproblems. One is to seek the solution of VRP for collecting product 1, and another is for collecting product 2. The sum of the total distance for the two subproblems represents the final solution listed in column 3 of Table4. The results obtained by using two-compartment vehicles are the same as reported in column 5 of Table 2. The increase percentage of total length listed in column 5 of Table 4 can be obtained from the following equation:(29)increase%=SC−TCSC×100,where SC and TC are the total distance travelled by single-compartment vehicles and two-compartment vehicles, respectively. Besides, in order to make the comparison more intuitive, plot with 95% Tukey’s honest significant difference (HSD) confidence intervals for employing these two types of compartments under two categories of data sets is shown in Figure 3.Table 4
Comparison results of using different compartment vehicles.
Inst.nTotal length of using single compartmentTotal length of using two compartmentsTotal length increase percentage (%)vrpnc1a50932.13550.7069.26vrpnc1b50934.13548.7770.22vrpnc2a751317.96873.5750.87vrpnc2b751327.43888.3849.42vrpnc3a1001423.90870.1163.65vrpnc3b1001415.06871.2262.42vrpnc4a1501709.251110.0653.98vrpnc4b1501690.641134.9048.97vrpnc5a1992043.761423.1143.61vrpnc5b1992049.471456.1640.74vrpnc11a1201583.431110.4542.59vrpnc11b1201556.181218.6627.70vrpnc12a1001281.69907.7341.20vrpnc12b1001305.50961.7435.74Average1469.32994.6847.72Figure 3
Comparison results of vehicles with different compartments.It can be obviously seen from Table4 and Figure 3 that using single-compartment vehicles brings the notable increase in the total length of each route. The reason is that although the total route cost is reduced due to the increase of vehicle capacity while transporting the same product, each customer must be served twice in order to ensure that different products are all delivered. More specifically, it can be noticed from Table 4 that the total distance raises 47.72% on average when the vehicles with the single compartment are employed to visit the customers. Moreover, it is worth mentioning that the total distance could be reduced by up to 70.22%, while two-compartment vehicles are selected. The statistical results clearly display the profit of using multi-compartment vehicles.


## 5. Conclusion

The goal of this paper is to propose an alternative solution for the NP-hard problem, which is dedicated to minimize the total travelled mileage for the multi-compartment vehicle routing problem (MCVRP). A novel improved hybrid ant colony algorithm (IHACO) is proposed, aiming at promoting the global search to reach different regions of high-quality solutions. A new probabilistic model which utilizes two probability matrices to learn different useful information of high-quality solutions is developed. Then, a heuristic rule based on the nearest neighbor is used to initialize this probabilistic model so that the algorithm’s speed of reaching the promising regions is enhanced. Meanwhile, a new local search is constructed to execute the effective exploitation. The efficiency of this method is reflected in the fact that it does not need to judge the capacity constraint and evaluate the objective function. It only needs to judge whether the two routes cross or not, and then, it can directly carry out the exchange operation between two customers. Moreover, two types of VND combined with the speed-up strategy and the first move strategy are employed to further perform the deep, narrow but rich search in the good solution space that has been found. Finally, numerical experiment results based on a benchmark dataset display the effectiveness and efficiency of IHACO. Future work is to develop an effective ACO-based algorithm to deal with the multiobjective MCVRP and the fuzzy travel time-based MCVRP.



---

*Source: 8839526-2020-10-21.xml*