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

,


Introduction
e 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. ere are many industrial fields in which the multi-compartment vehicles are employed. e 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][2][3][4][5][6][7]. Another example is the transportation of food, where different degrees of refrigeration goods are stored in different compartments in one vehicle [8][9][10][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][13][14]. Despite the practical applications, the solution of MCVRP has been rarely studied.
us, 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][16][17][18][19][20]. However, with the increasing problem size, these algorithms are usually of limited use owing to their excessive timeconsuming nature. en, a variety of meta-heuristic algorithms were proposed to obtain desirable solutions within a reasonable calculation time [21][22][23][24]. ere are many intelligent algorithms for solving nonlinear problems [25][26][27]. e ant colony optimization algorithm (ACO) is one of the highly efficient meta-heuristic or intelligent algorithms. e 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. e ants exchange information with other ants through the pheromone. e 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. e existing results usually use only pheromone concentration matrix (hereinafter referred to as pheromone matrix) to construct a probabilistic model to generate new individuals. e 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/C nn ), where (C nn ) 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. e pheromone matrix was adopted merely and initialized to (1/nL 0 ), where n is the number of customers and L 0 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 (ph min ). 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/nL neig ), where n is the customers' number and (L neig ) 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. e 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. is 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. e 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. en, 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. e remainder of this paper is organized as follows: Section 2 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.

Problem Description and Model Formulation
e MCVRP is defined on an undirected graph (G � (N, A)) with a set of vertices N � 0, 1, . . . , n { } and a set of Complexity from i to j. Node 0 in N represents the depot, and (N c � (N/ 0 { })) are nodes of customers. e depot has an available homogeneous fleet. ese 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 ∈ N c ) has the quantity q ip to be carried for the product p, and the compartment capacity reserved for this product does not exceed Q p . 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 Table 1.
en, the integer linear program of the MCVRP can be described as follows: subject to i∈N c 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.

IHACO for the MCVR
is 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.

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 (π ant(i) ) � [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.

Probabilistic Model and Its Initialization.
e 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.
is initialization method makes the algorithm converge to the high-quality solutions faster.

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 ).
is 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. en, the probabilistic model in this paper contains a pheromone matrix and a sequence information matrix.
In a subpath visited by one vehicle, let (s kj ) 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 where N k 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 π ant(1) � [0, 1, 2, 3, 0, 4, 0], π ant(2) � [0, 2, 3, 4, 0, 1, 0], π ant(3) � [0, 3, 4, 0, 1, 2, 0], and π ant(4) � [0, 3, 2, 1, 4, 0], and then, s kj is given as follows: Since the first location is fixed as the depot, the corresponding probability values of μ kj for the second position of each customer in en, 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:

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 top t � round(n/(10)) customers into set c first . Let c first � [c first 1 , c first 2 , . . . , c first t ] denote the first customer to be visited by the t-th ant, R � [R 1 , R 2 , . . . , R t ] the complete path passed by the t-th ant (also the t-th initial solution), and L � [L 1 , L 2 , . . . , L t ] the length of the route taken by the t-th ant. e procedure of the probabilistic model initialization is given as follows: Step 1: set t � 1 and (τ 0 ij � zeros(n × n)) and (s 0 kj � zeros(n × n)).
Step 2: the ant starts from the depot 0 and visits the first customer c first t . en, the ant selects the remaining customers for access based on the nearest neighbor principle. is process forms the t-th initial solution R t .
Step 3: calculate the route length L t corresponding to R t .
Step 7: output s t kj and τ t ij .

Ant Path Construction.
In this study, each ant produces a complete solution for the MCVRP. e 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. Let m 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. en, 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). e probability of choosing customer j as the next one in the t-th iteration is given by formula (13): 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: Set q 0 be a constant value and q a random variable which is uniformly distributed in [0, 1]. If q > q 0 , select the next customer in accordance with the formulation of p t ij ; otherwise, select the customer who has the maximum value of δ ij [22].

Symbol
Description of symbol Distance between nodes i and j q ip Quantity of product p to be carried from node i Q p Capacity of the vehicle compartment dedicated to product p x ij Equals 1 if the vehicle travels from node i to j and is zero otherwise Total quantity of product p while the carried vehicle departing from node i 4 Complexity

Update Mechanism of Probability Matrices.
e probabilistic model consists of two probability matrices s kj and τ ij . ey 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: In equation (14), L t best is the total length corresponding to the optimal ant path in the t-th iteration: In equation (15), s kj (π t best ) is the probability matrix s kj obtained by the t-th optimal individual π t best . We refer to Section 3.2.1 for more details.

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. e algorithm in this section performs the problem-dependent local search on the global optimal solution π gen best , for the purpose of realizing the detailed search of the area near π gen best . 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.

Basic and Pipe VND.
e VND is a variant of the variable neighborhood search. e 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 first 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.

Operation of Insert and Swap. Assume an ant's tour is represented as
. e insert operator selects one customer like i and inserts it before or after another customer like j. en, the tour becomes π new_insert ant . e swap operator selects two customers like i and j on this tour and swaps them.

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 Figure 1, there is one vehicle route denoted by π old ant_vehicle � [0, 1, 3, 2, 4, 0]. For any two given intersecting line segments, such as lines L 13 and L 24 which intersect at point u, the end points of the two segments that can form the opposite side are connected like L 12 and L 34 . Next, the new path without intersection is obtained and denoted by π new ant_vehicle � [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), (L 1u + L 2u > L 12 ); in Δ3u4, (L u3 + L u4 > L 34 ). Let L old total be the total length of vehicle path with intersecting lines, then L old total � L 01 + L 1u + L u3 + L 32 + L 2u + L u4 + L 40 . Let L new total be the total length after eliminating intersecting lines, then L new total � L 01 + L 12 + L 23 + L 34 + L 40 . As a result, From π old ant_vehicle and π new ant_vehicle 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 Figure 1 are (x 1, y 1 ), (x 2, y 2 ), (x 3, y 3 ), and (x 4, y 4 ), and then, the linear equations of line segments L 13 and L 24 are shown in the following equations: In formulas (16) and (17), λ is the slope of the line where L 13 is located, and ψ is the slope of the line where L 24 is located.
If L 13 and L 24 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: □ 3.5.4. Speed-Up Search Strategy. Let insert(π, 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 IJ infeasible insert 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 IJ infeasible swap 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 N insert (π, i) denote the feasible neighborhood set based on insert operation for i customers in solution π, N swap (π, i) the feasible neighborhood set based on swap operation for i customers in solution π, and N elim_intersection (π, k) the feasible neighborhood set based on intersection elimination operation for customers serviced by the k-th vehicle in solution π. Set n k to be the number of customers visited by the k-th vehicle, then (l � k i�1 n i − n k + k + 1) is the position of the first customer visited by the k-th vehicle in solution π. Hence, there are as follows: It can be seen from (23)- (25) that N insert (π, 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. N swap (π, i) is a feasible neighborhood set that any customer i swaps with any possible customer j within the same vehicle or between different vehicles. N elim_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 customer j 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. us, the more compact neighborhood sets N narrow1 insert (π, i) and N narrow1 swap (π, i) can be obtained: 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. Let L(π) 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 π.  Complexity 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. erefore, a more compact and effective neighborhood set N narrow2 insert (π, i) can be obtained on the basis of N narrow1 insert (π, i): ForN narrow1 swap (π, i), π[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. erefore, a more compact and effective neighborhood set N narrow2 swap (π, i) can be obtained on the basis of N narrow1 swap (π, i): For N elim_intersection (π, k), we can use the method described in Section 3.5.3 to determine whether two paths intersect. If so, according to eorem 1, we know L( en, 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.
Let SP FindBestN narrow2 insert (π, i), SP FindBestN narrow2 swap (π, i), and SP FindBestN elim_intersection (π, k) be the optimal neighborhoods in searching and outputting N narrow2 insert (π, i), N narrow2 swap (π, i), and N elim_intersection (π, k), respectively. ese 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. e corresponding procedures are given as follows: Step 1: set i � 1 and π previous best � π gen best .
Step 4: if i < n − 1, then i � i + 1 and go to step 2.
Step 5: if i � n − 1, then π current best � π gen best . If π current best < π previous best , go to step 1.
Step 9: if i < n − 1, then i � i + 1 and go to step 7.
Step 10: if i � n − 1, then π current best � π gen best . If π current best < π previous best , go to step 6.
Step 14: if k < K, then k � k + 1 and go to step 12.
Step 15: if k � K, then π current best � π gen best . If π current best < π previous best , go to step 11.

First Move Strategy.
In order to reach more different regions, a first move strategy is adopted in SP FindBestN narrow2 insert (π, i), SP_FindBestN narrow2 swap (π, i), and SP_FindBestN elim_intersection (π, k). Here is an example to illustrate the strategy. While the first neighbor π i first in N narrow2 insert (π, i) which can improve π gen best is gained, SP_FindBestN narrow2 insert (π, i) terminates and outputs π i first . Furthermore, the sequence π in SP_FindBestN narrow2 insert (π, i) is replaced by π gen best for the purpose of implementing the search instantly from the promising areas discovered by SP_FindBestN narrow2 insert (π, 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. e changed procedures with the above strategy, named SP_FindFirstMoveN narrow2 insert (π, i), SP_FindFirstMoveN narrow2 swap (π, i), and SP_FindFirstMoveN elim_intersection (π, k), are only different from SP_FindBestN narrow2 insert (π, i), SP_FindBestN narrow2 swap (π, i), and SP_FindBestN elim_intersection (π, k) in step 2, step 7, and step 12.

Procedure of Local Search. Let SPVND_FindFirstMove
Step 1. Perturbation phase.
Step 4: output π gen best . In the above procedure, step 1 is a perturbation phase, which is helpful to prevent the cycling search and overcome the local optima. e special devised step 2 performs exploitation from the region obtained by step 1.

IHACO.
Denote gen 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 parameter m, maxgen, α, β, ρ, ρ 0 , S, and KP.
Step 2: probabilistic model initialization. Set gen � 1//Section 3.2.2 Step 3: ant path construction. Generate the population at generation gen by sampling σ kj and then evaluate its individuals to obtain the global best individual π gen best //Section 3.3 Step 4: local exploitation. Execute three types of local search to π gen best and update π gen best //Section 3.5 Step 5: update mechanism of probability matrices. Calculate the objective function value L(π gen best ) corresponding to the best individual π gen best . Update π ij and s kj based on L(π gen best ) and π gen best //Section 3.4 Step 6: set gen � 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 modelbased 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.

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. e experiments implemented have been coded in Matlab 2019a and tested on an Intel 1.8 GHz PC with 12 GB RAM.

Experimental Setup.
e generation method of instances is derived from [13]. e test set in this paper contains 14 benchmark instances with the capacity constraint, including vrpcn1-vrpcn5 and vrpcn11-vrpcn12 of sets a and b. ese scales of problem range from 50 to 200 customers. e 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 q 0 � 0.9. ere 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 (t best ) and the running environment (that is, the CPU frequency(f best )) 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 8 Complexity the proposed algorithm t � t best × (f best /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.

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 twocompartment 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. e experimental results of the five algorithms are shown in Table 2. e first column in Table 2 represents the type of the instance, and the second column is the quantity of customers contained in every instance. e meaning of Obj. and T(s) in this table is the objective function value and runtime, respectively.
As reported from Table 2, 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   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.

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, s kj 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/L neig ), and s kj is initialized just by the individual π neig (the individual corresponding to the objective function value L neig ). (3) IHACO_V3 : in IHACO's local search, N narrow1 insert (π, i) and N narrow1 swap (π, 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. is 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.
e statistical results are listed in Table 3. Besides, the graphical results for the compared algorithms are exhibited in Figure 2.
Compared with IHACO's variants, it is clear from Table 3 and Figure 2 that IHACO performs best for all instances. e 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.
is once again verified a conclusion of the previous experiment; that is, IHACO performs better when the number of customers is large.
e 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 highquality solutions than IHACO. IHACO_V4 must find the optimal location to execute the insert and swap operations.
is 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. us, the main operations removed from these six variants are all adopted by IHACO.

Superiority Verification of Using Multi-Compartment
Vehicles. Since some different products have different characteristic properties, they cannot be mixed. at 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. e sum of the total distance for the two subproblems represents the final solution listed in column 3 of Table 4. e results obtained by using twocompartment vehicles are the same as reported in column 5 of Table 2. e increase percentage of total length listed in column 5 of Table 4 can be obtained from the following equation: where SC and TC are the total distance travelled by singlecompartment 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. It can be obviously seen from Table 4 and Figure 3 that using single-compartment vehicles brings the notable increase in the total length of each route. e 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. e statistical results clearly display the profit of using multicompartment vehicles.

Conclusion
e 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. en, 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. e 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.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.