Mutation Ant Colony Algorithm of Milk-Run Vehicle Routing Problem with Fastest Completion Time Based on Dynamic Optimization

Abstract

The objective of vehicle routing problem is usually to minimize the total traveling distance or cost. But in practice, there are a lot of problems needed to minimize the fastest completion time. The milk-run vehicle routing problem (MRVRP) is widely used in milk-run distribution. The mutation ACO is given to solve MRVRP with fastest completion time in this paper. The milk-run VRP with fastest completion time is introduced first, and then the customer division method based on dynamic optimization and split algorithm is given to transform this problem into finding the optimal customer order. At last the mutation ACO is given and the numerical examples verify the effectiveness of the algorithm.

1. Introduction

The vehicle routing problem (VRP) was firstly brought forward by Dantzig and Ramser in 1960 [1]. With the development of VRP, there exist several variations and specializations. Mostly VRP aims to minimize the total travel distance (or travel time) and total cost. But for distribution of fast foods, express delivery, and emergency supplies, these objectives are not suitable and the completion time is more important.

Though the vehicle routing problem with time windows takes the service time into consideration [2], it cannot solve the problem when all customers want to be served as early as possible. Nikolakopoulou et al. solved a vehicle routing problem by balancing the vehicles time utilization [3]. But it maybe takes a long time for some vehicles to complete the distribution tasks when the vehicles time utilizations are balanced. The VRP with fastest completion time was studied to minimize the fastest completion time [4–6].

Traditional vehicle routing problem assumes that the total distribution can be completed by one vehicle in a round trip. In practice, the number of vehicles is limited and there are many customers to be served. So a round trip by one vehicle is impossible and more round trips and vehicles are needed, such as the tobacco distribution with thousands of customers. As JIT production, small-scaled and multiple-batch distributions are more popular, in which one vehicle is needed to collect goods on multiple round trips. This kind of vehicle routing problem is called milk-run vehicle routing problem (MRVRP).

For problem aiming to minimize total travel distance or total cost, MRVRP can be transferred into the traditional VPR by increasing the routes number and allowing a vehicle serving customers in different routes. But for problem to minimize the fastest completion time, the completion time in the first period affects that in the second period. It is multiperiod optimization problem, so the algorithm of the MRVRP with fastest completion time will be studied based on dynamic optimization.

VRP is an NP-hard problem [7]; heuristics and evolutionary algorithms are used to solve VRP. In this paper, mutation ant colony algorithm is used to solve MRVRP with fastest completion time. In the next section, we will give the description of the MRVRP with fastest completion time. In Section 3, optimal division methods of customer orders based on dynamic optimization are given. In Section 4, mutation ant colony algorithm is given to solve the problem. In the last section, a numerical example is given.

2. MRVRP with Fastest Completion Time

Suppose that there is one depot serving n customers, which have 𝑚 vehicles with capacity 𝑤. The demand of customer 𝑖 is 𝑞_𝑖(𝑖=1,…,𝑛) and satisfies (𝑘−1)𝑤≤∑^𝑛_𝑖=1𝑞_𝑖<𝑘𝑤, which means that for every vehicle, it needs 𝑘 times to finish the distribution tasks on average.

Because the distances between the depot and the customers and the distances between the customers are different, the travel time between any two nodes is different. Supposing that the traveling time between node 𝑖 and node 𝑗 is 𝑡_𝑖𝑗,𝑖=0,1,…,𝑛,𝑗=0,1,…,𝑛, node 0 is depot and the times satisfy the triangle inequality.

The vehicle routes satisfy the following.(i)Each customer is served by a certain vehicle.(ii)One vehicle can serve many customers. When the demands of the customers exceed the vehicle capacity, the vehicle returns to the depot to unload and goes back to serve the next customer.(iii)The demands of the customers served by a vehicle cannot exceed 𝑘 times vehicle capacity.(iv)All vehicles should depart from the depot and return to the depot.

When the demands of the customers exceed the vehicle capacity, the customers should be divided into different groups. The customers in the first group are served first, then the second group, until all the customers are served, which is shown in Figure 1.

Figure 1 

The milk-run vehicle routes of two vehicles.

For the MRVRP aiming to minimize total travel distance or total cost, because the traveling distance of a vehicle is the sum of the distances of each round and the total distance is the sum of all vehicle traveling distances, the MRVRP can be transferred as the VRP with km vehicles. The customers are divided into km routes and each vehicle serves 𝑘 routes. The total traveling distance of 𝑚 vehicles equals the total distance of km routes. Then if the total distance of km routes is minimized, the total traveling distance of 𝑚 vehicles is also minimized.

But for the MRVRP with fastest completion time, the completion time of the whole distribution task is the completion time of the last vehicle. Suppose that the completion time of vehicle 𝑙 is 𝑇_𝑙,𝑙=1,2,…,𝑚, and the completion time of the whole distribution task is 𝑇; then we have, 𝑇=max_{𝑙=1,2,…𝑚}𝑇_𝑙. (1)

For each vehicle, the distribution task needs to complete several routes and the completion time is determined by the time when last customer served. Because there is order of each route, the completion time of vehicle 𝑖 is the sum of every route completion time. In the last route, the time from the customer to the depot is not considered.

If this problem is transferred to the problem with km routes in which each vehicle is in charge of 𝑘 routes, the completion time is the sum of the travel times of previous 𝑘−1 routes and the completion time of the last routes. Supposing that the 𝑙th vehicle is in charge of the distribution task in 𝑘 routes, the traveling time of the 𝑗th route is 𝑎_𝑗,𝑗=1,2,…,𝑘, and the traveling time from the last served customer to the depot is 𝑡, the completion time of this vehicle is 𝑇_𝑙=^𝑘_𝑗=1𝑎_𝑗−𝑡. (2)

The difficulty of solving VRP lies in too many arrays of customer service order. It is hard to solve VRP by dynamic programming since there are too many states. For a given array of customers, the problem is transferred to how to divide the customers into groups and the customers in the same group are served by one vehicle. The feasible groups’ division scheme is much less and it is easy to be solved by dynamic programming. In this paper, we will fix the serving array first and then give the optimal division scheme and calculate the fastest completion time. Taking this completion time as the objective, we will determine the final customer serving order by ACO.

3. Customer Division Method Based on Dynamic Programming

The problem of dividing customers into groups is how to divide a given serving array 𝑖₁,𝑖₂,…,𝑖_𝑛 into m groups to ensure that the completion time of all the distribution tasks is minimized and the total demands of the same group do not exceed kw. The customers belonging to the same group are served by one vehicle on the array order.

Usually there are a lot of division schemes and their completion times are different. Nikolakopoulou gave the split method to minimize the total traveling distances by transferring the dividing problem into the shortest-path problem [3]. We gave the improved split method to minimize the fastest completion time by transferring the division problem into the longest-edge shortest problem [5]. The improved split method spends a long calculation time, as it repeats to find the shortest paths. Then the customer division method based on dynamic programming is given which spends less time [6].

For the MRVRP, there are two divisions: the first is to divide the customers into 𝑚 groups and the second is to divide the customers in the same groups into 𝑘 routes. The aim of the first division is to minimize the longest completion time and is solved with dynamic programming. The aim of the second division is to minimize the total traveling time and is solved with split method.

3.1. The First Division of Customers

The problem of dividing customers into 𝑚 groups is to determine 𝑚−1 cut points and can be looked as at 𝑚-period decision problem. In a period, the number of customers served by a vehicle is determined. So the state variable is the number of customers who have not assigned, denoted as 𝑠, and the number of vehicles which have not assigned, denoted as 𝑥. Then the state variable is (𝑠,𝑥), and the initial state is (𝑛,𝑚). The decision variable is the number 𝑧 of customers served by the next vehicle.

In the next period, the numbers of unassigned customers and vehicles are 𝑠−𝑧 and 𝑥−1. The state variable is (𝑠−𝑧,𝑥−1). The state variable and decision variable are both confined by the vehicle capacity. For a given (𝑠,𝑥), the demands of assigned customer cannot exceed 𝑘 times vehicle capacity. Supposing that 𝑑_𝑖(𝑖=1,…,𝑛) is the total demand from the first customer to the 𝑖th customer, we have 𝑑_𝑛−𝑠≤(𝑚−𝑥)𝑘𝑤. (3)

When the decision variable in the present period is 𝑧 and the unassigned customers need 𝑥−1 vehicles to distribute, the total demand of the 𝑧 customers cannot exceed 𝑘𝑤; that is, 𝑑_𝑛−𝑑_{𝑛−𝑠+𝑧}≤()𝑑𝑥−1𝑘𝑤,_{𝑛−𝑠+𝑧}−𝑑_𝑛−𝑠≤𝑘𝑤. (4) Denote decision variable set as Ω; then Ω=𝑧∈Ζ∣𝑑_𝑛−𝑑_{𝑛−𝑠+𝑧}≤(𝑥−1)𝑘𝑤,𝑑_{𝑛−𝑠+𝑧}−𝑑_𝑛−𝑠≤𝑘𝑤. (5)

Supposing that 𝑓(𝑠,𝑥) is the fastest completion time of serving the 𝑠 customers with 𝑥 vehicles and ℎ(𝑠,𝑧) is the completion time if the former 𝑧 customers are served by a vehicle, we have 𝑓(𝑠,𝑥)=min_𝑧∈Ωmaxℎ(𝑠,𝑧),𝑓(𝑠−𝑧,𝑥−1). (6)

When there is a vehicle left, if the total demands of 𝑠 customers do not exceed kw, the 𝑠 customers are served by this vehicle and the completion time is ℎ(𝑠,𝑠); otherwise the completion time is infinite; that is, 𝑓()=ℎ()𝑠,1𝑠,𝑠,if𝑞_{𝑛−𝑠+1}+⋯𝑞_𝑛≤𝑤,+∞,if𝑞_{𝑛−𝑠+1}+⋯𝑞_𝑛>𝑤. (7)

The state variable and decision variable are discrete. It can be solved by enumeration method. In a period, the vehicle number decreases one unit and the unassigned vehicle is determined by the periods. Then the number of states is determined by the number of unassigned customers, the number of states of each period does not exceed 𝑛. Since the upper bound of vehicle capacity is 𝑘𝑤, there is also an upper bound of the number of customers to be served, denoted as 𝑔. The maximum iteration time of the method is 𝑚𝑛𝑔. In practice, since the state variable should satisfy the inequality (3) and the decision variable should satisfy the inequality (4), the actual iteration time is much less than 𝑚𝑛𝑔.

In each period, the completion time of every vehicle is needed to calculate and the second division of customers is needed.

3.2. The Second Division of Customers

Suppose that the customers served by a given vehicle are 𝑖₁,𝑖₂,…,𝑖_𝑧. Because the total demands exceed the vehicle capacity 𝑤, the route should be cut into smaller groups to ensure that the customers in the same group can be served by the same vehicle. This is also a division problem which is different from that in Section 3.1 in the scale. The scale is smaller and the objective is not to minimize the maximal completion but to minimize the total milk run completion time.

To get the division schemes, a directed graph is constructed. The vertex set has the depot and the customers served by the same vehicle. For nodes 𝑗₁ and 𝑗₂, supposing that 𝑖_𝑗1 is in front of 𝑖_𝑗2, the total demand of all customers between 𝑖_𝑗1and 𝑖_𝑗2 (Including customer 𝑖_𝑗1but not 𝑖_𝑗2) is calculated. If the total demand is less than 𝑤, an arc from𝑗₁ to 𝑗₂ is drawn. The weight is the traveling time from the depot in turn reach all customers between 𝑖_𝑗1and 𝑖_𝑗2 (including customer 𝑖_𝑗1but not 𝑖_𝑗2), and then return. The weighted direct graph is called split graph which is shown in Figure 2.

Figure 2 

Split graph.

The shortest distance from 1 to 𝑧+1 is the shortest travel time of all distribution tasks and the fastest completion time is this time minus the return time from customer 𝑖_𝑧 to depot.

Suppose the capacity of a vehicle is 15 and the order of customers to be served by the vehicle is 1-2-5-3-4-7; the traveling times between the nodes are shown in Figure 3, where the numbers in the brackets are the demands.

Figure 3 

The distribution of customers.

The total demand of customers 1, 2, and 5 is no more than 15; then an edge is drawn between node 1 and 3 with the weight 90. The other weights are got similarly. After calculation, the directed weighted graph is got in Figure 4.

Figure 4 

Split graph.

The shortest path is from node 1 to node 3 plus from node 3 to node 0. The optimal distance is 180. Customers 1, 2, and 5 construct the first route. Customers 3, 4, and 7 construct the second route. The completion time is 155(=180−25).

The algorithm to solve the shortest path problem such as Dijkstra algorithm is polynomial. Its complexity is 𝑂(𝑛²), where 𝑛 is the number of nodes. For this problem, the maximum number of nodes is 𝑔 and the complexity is 𝑂(𝑔²). Because for every iteration the shortest path problem should be solved to calculate the completion time of a vehicle, the complexity of the whole algorithm is 𝑂(𝑚𝑛𝑔³).

4. Mutation ACO Algorithm

In this section, we will give an improved ACO to solve the optimal customers array.

4.1. ACO Algorithm

The main idea of ant colony algorithms is to mimic the pheromone trail used by real ants as a medium for communication and feedback among ants. Basically, the ACO algorithm is a population-based, cooperative search procedure that is derived from the behavior of real ants. ACO algorithms make use of simple agents called ants that iteratively construct solutions to combinatorial optimization problems. The key problem to solve VRP with ACO is how an individual ant constructs a complete solution by starting with a null solution and iteratively adding solution components until a complete solution is constructed. The key problem of ACO is to determine the pheromone matrix. The pheromone matrix is (𝑛+1)×𝑛 when there are 𝑛 customers, where the last row stands for the information from the depot to the customer and the 𝑖th row stands for the information from customer 𝑖 to other customers. Initially, since for a given customer there is the same possibility following other customers, the pheromone matrix starts with equal probable matrix. When 𝑡=0, the pheromone matrix is 𝐵_𝑖𝑗(0)=01⋯1𝑛−11𝑛−1⋯⋯⋯⋯1𝑛−11𝑛−1⋯0𝑛1𝑛⋯1𝑛_{(𝑛+1)×𝑛}, (8) where 𝑏_𝑖𝑗,𝑖=1,2,…,𝑛+1,𝑗=1,2,…,𝑛 is the pheromone from 𝑖 to j.

In the 𝑡th iteration, there are 𝑚 ants and 𝑚 customer arrays can be got. Supposing that the fastest completion time of each array is 𝐿_𝑘(𝑡);𝑘=1,2,…,𝑚, the increased value of pheromone in the 𝑡th iteration is Δ𝑏_𝑖𝑗(𝑡)=∑^𝑚_𝑘=1Δ𝑏^𝑘_𝑖𝑗(𝑡), where Δ𝑏^𝑘_𝑖𝑗(𝑡)=𝑄𝐿_𝑘(𝑡);if𝑗follows𝑖closely,where𝑄isaconstant,0;otherwise. (9) Local update is performed during the ant constructive procedure in the following way: 𝑏_𝑖𝑗(𝑡+1)=1−𝜌𝑏_𝑖𝑗(𝑡)+Δ𝑏_𝑖𝑗(𝑡), (10) where 𝜌 is the evaporation coefficient.

The heuristic information is in an (𝑛+1)×𝑛 matrix, with the last row 𝜂_𝑛+1,𝑗=𝐴/𝑡_𝑗 and the other rows 𝜂_𝑖𝑗=𝐴/𝑡_𝑖𝑗, where 𝑡_𝑖𝑗 is the traveling time between two nodes, 𝑡_𝑗 is the average traveling time from customer 𝑗 to the depot, and 𝐴 is a constant depending on the situations.

In a given service order, the first served customer is determined by probability. The next served customer is selected from the allowing set allowed_𝑖, which is the customers set that can be served after 𝑖. The customer 𝑗 is chosen with the probability 𝑝_𝑖𝑗, where 𝑝_𝑖𝑗(𝑡)=𝑏_𝑖𝑗(𝑡)^𝛼𝜂_𝑖𝑗^𝛽∑_{𝑠∈allowed𝑖}𝑏_𝑖𝑠(𝑡)^𝛼𝜂_𝑖𝑠^𝛽;𝑗∈allowed_𝑖,0;𝑗∉allowed_𝑖. (11)𝛼 is the importance of pheromone information, and 𝛽 is the importance of heuristics information. After calculation of the probability matrix, the node is selected by the probability 𝑝_𝑖𝑗 and the node is deleted from the allowing set. Then a customer service order is got with the pheromone.

4.2. Mutation Operator

ACO algorithm possibly runs into prematurity just as other evolutionary algorithms. The main reason is the concentrations of pheromone which makes the same solutions be got. In reality, the pheromone may be changed by rain and other factors and this change may help the ants find a new route. A mutation operator will be introduced into the pheromone of ACO to escape from local optima and strengthen its global search ability.

The initial pheromone matrix’s factors are the same. As time goes on, some elements become large but other elements become small in the same row. The pheromone is concentrated which obstructs the ability to find more optimal solutions. The concentration of a row is defined as the ratio of the maximum factor and the sum of the factors in the row. The 𝑖th row concentration 𝜇_𝑖(𝑡) is 𝜇_𝑖(𝑡)=max_{𝑗=1,2,…,𝑛}𝜏_𝑖𝑗∑^𝑛_𝑗=1𝜏_𝑖𝑗,𝑖=1,2,…,𝑛+1. (12)

The concentrations of different row are not the same. Whether to mutate a row should depend on the concentration of the row. When the concentrations of a row exceed the threshold value, the factors in the row should be mutated.

For the row mutation, the maximum factor decreases randomly and the decreased value is assigned to other positions. Suppose 𝑗₀th factor is the maximum factor in the 𝑖th row, and given a random number 𝑟(0≤𝑟≤1) and a random vector 𝑅_1×𝑛. A new row is got with the 𝑗₀th factor multiplied by 𝑟, and the decreased value is assigned to other positions; that is, 𝜏^′_𝑖=𝜏^′_𝑖+𝜏_𝑖𝑗0(1−𝑟)∑^𝑛_𝑗=1𝑅_𝑗𝑅. (13)

Row mutation is a local mutation and only affects the choice of routes partly. When the matrix has high concentration, local mutation cannot ensure escaping the local optima and matrix mutation is necessary. The minimum concentration of the row is the matrix concentration, which is 𝜈(𝑡): 𝜈(𝑡)=min_{𝑖=1,2,…,𝑛+1}𝜇_𝑖(𝑡). (14)

Whether to mutate a matrix should also depend on the concentration of the matrix. When the concentrations of a matrix exceed the threshold value, the matrix should be mutated. For the matrix mutation, every element in the matrix decreases randomly and the decreased value is again randomly assigned. A random number 𝑟(0≤𝑟≤1) and a random vector 𝑅_𝑛+1×𝑛 are got. New values are got with every factor multiplied by 𝑟. Then the decreased value is assigned randomly’ that is, 𝜏^′=𝜏𝑟∗𝜏,_′=𝜏^′+∑^𝑛+1_𝑖=1∑^𝑛_𝑗=1𝜏_𝑖𝑗−𝜏^′_𝑖𝑗∑^𝑛+1_𝑖=1∑^𝑛_𝑗=1𝑅_𝑖𝑗𝑅. (15)

Then a new pheromone matrix is got and the mutation can ensure that the total pheromones do not change.

The mutation ACO algorithm is as follows.

Step 1. Initiation: Determine the parameters 𝑚,𝑇,𝛼,𝛽,𝜀,𝜙, and 𝑝. Input the initial pheromone 𝜏_𝑖𝑗(0) and heuristics matrix 𝜂_𝑖𝑗 and get the initial ants. Give 𝑚 customers arrays randomly as 𝐴(0)=𝐴₁(0),𝐴₂(0),…,𝐴_𝑚(0). (16)

Step 2. For a given array, assign the customers to the vehicle and get the corresponding completion time by dynamic programming and split algorithm. Record the present best route besttrip and the fastest completion time finishtime.

Step 3. The local pheromone update is performed by all the ants after each construction step on the formulas (9) and (10).

Step 4. Row Mutation: Inspect whether the pheromone concentration of the row in the pheromone matrix is more than 𝜀; if yes, a random number 𝑟 is got. If 𝑟≤𝑝, a new row is got by (13).

Step 5. Matrix Mutation: Inspect whether the pheromone matrix concentration is more than 𝜙; if yes, a random number 𝑟 is got. If 𝑟≤𝑝, a new pheromone matrix is got by the formula (15).

Step 6. According to the new pheromone matrix 𝐵_𝑖𝑗(𝑡) and the heuristics information matrix 𝜂_𝑖𝑗, Give 𝑚 customers arrays randomly as 𝐴(𝑡+1)=𝐴₁(𝑡+1),𝐴₂(𝑡+1),…,𝐴_𝑚(𝑡+1). (17)

Step 7. Let 𝑡=𝑡+1. ACO procedure stops if 𝑡=𝑇 and output the finishtime and besttrip; otherwise return to Step 2.

5. Numerical Examples

In this section, we will consider the emergency supplies distribution after the earthquake in Wenchuan. After the earthquake, the roads on earth are not fluent. The emergency supplies are mainly transported by helicopter. As the number of helicopters is limited comparED to the broad place, the helicopter 𝑠 is needed for distribution many times. The problem is a typical MRVRP with fastest completion time.

Suppose there are 3 helicopters in charge of the distribution of emergency supplies to 20 settlements. The location of material distribution center is (30, 40), the demands and location of the settlements are given in Table 1.

Table 1 

Locations and demands of settlement.

The capacity of a helicopter is 12t. The total demand is 92t. Every plane is needed to fly 3 times on average and distribute at most 5 settlements every time.

The algorithm is realized by Scilab. The parameters are as follows: the ants number is 10 and 𝛼=1.2,𝛽=0.5,𝑃=0.15,𝜀=0.8,𝜙=0.75. The iteration time is 400. Run the program 20 times for the same problem, and get the fastest completion which times are as shown in Figure 5.

Figure 5 

20 times calculation results.

Average fastest completion time is 244.9978, their maximum gap is 7.291, and the gap is 3.1014% of the best scheme. It shows that the algorithm has good convergence.

Best distribution route is as shown in Figure 6.

Figure 6 

Optimal scheme of MRVRP with fastest completion time.

Its fastest completion time is 235.08428 and total travel time is 865.66435. A new scheme to minimize total travel time is as shown in Figure 7.

Figure 7 

Optimal scheme.

Its shortest total travel time is 845.89676, the fastest completion time of this customer array is 273.379, and it is 38.295 minutes more than the first scheme. And the total flight time is reduced by 19.76759 minutes. In emergency management, effectively shortening the completion time is very necessary, and therefore the first scheme is more reasonable.

6. Conclusions

In this paper, the MRVRP with fastest completion time is proposed, which has many applications in fast foods distribution, express delivery, and emergency supplies. Solving the problem is more difficult than the general VRP. The key problem to solve MRVRP with fastest completion time is to give the division method for customer array. The customer division method based on dynamic programming and split method is given in this paper, which can transfer MRVRP with fastest completion time into the problem of finding the optimal customer service order. Then the problem is solved with mutation ACO.