# A Two-Phase Heuristic Algorithm for the Common Frequency Routing Problem with Vehicle Type Choice in the Milk Run

Authors: Yu Lin; Tianyi Xu; Zheyong Bian
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

High frequency and small lot size are characteristics of milk runs and are often used to implement the just-in-time (JIT) strategy in logistical systems. The common frequency problem, which simultaneously involves planning of the route and frequency, has been extensively researched in milk run systems. In addition, vehicle type choice in the milk run system also has a significant influence on the operating cost. Therefore, in this paper, we simultaneously consider vehicle routing planning, frequency planning, and vehicle type choice in order to optimize the sum of the cost of transportation, inventory, and dispatch. To this end, we develop a mathematical model to describe the common frequency problem with vehicle type choice. Since the problem is NP hard, we develop a two-phase heuristic algorithm to solve the model. More specifically, an initial satisfactory solution is first generated through a greedy heuristic algorithm to maximize the ratio of the superior arc frequency to the inferior arc frequency. Following this, a tabu search (TS) with limited search scope is used to improve the initial satisfactory solution. Numerical examples with different sizes establish the efficacy of our model and our proposed algorithm.


---

## Body

## 1. Introduction

A just-in-time (JIT) supply system managed parts transportation between suppliers and a manufacturer operating under the JIT discipline [1]. By simulating the JIT process from different perspectives, researchers showed that the JIT strategy could significantly improve efficiency and reduce cost [2–4]. With progress in research, the study of JIT has become more specialized. For example, in the implementation of JIT production for the manufacturer, a minimum inventory of raw materials is required to meet production needs, which in turn requires that the manufacturer supply parts in small and multifrequency batches according to operational parts consumption (speed). For inbound logistics, the popular milk run is well suited to manufacturers’ need for JIT production because of its characteristics of high frequency and small lot size, which enable it to help reduce the cost of inventory and transportation. Therefore, many manufacturers use the milk run as the main mode of transportation for inbound logistics.The milk run originated from the traditional system of milk distribution and sales in Western culture. In this system, a milkman simultaneously supplied customers with full bottles of milk and picked up the empty ones according to a predefined route. Over time, the high frequency and small lot sizes involved in this procedure made it attractive for use in manufacturing worldwide, since it was conducive to JIT production. The method has since developed into a popular one for collecting and delivering goods for multiple suppliers and manufactures using freight cars [5]. With respect to the milk run mode, researchers [6–9] currently focused on the vehicle routing problem (VRP). Dantzig and Ramser [10] first introduced the idea of the VRP. Since then, additional scholars have conducted research in this field. With subsequent research on the problem addressing practical applications according to varying constraints, the VRP now has several formulations. For instance, the vehicle routing problem with time windows (VRPTW) adds the constraint of a hard or soft time window based on the VRP, which has encouraged various solutions [11–14]. In addition to the time window, measuring the cost of inventory is a crucial factor for decision makers. Chien et al. [15] first used the cost of inventory as a factor in the vehicle routing problem and claimed that inventory allocation and the VRP were significant logistical decisions. Based on this premise, Chuah [16] discovered that frequency was affected by the inventory in the VRP and proposed a common frequency routing (CFR) problem. Based on the traditional VRP, the CFR problem simultaneously considers the relationship between frequency and inventory. Moreover, Chuah and Yingling [17] considered the amount in the inventory required to balance the relationship between inventory and frequency because low inventories increase the frequency of milk runs whereas high inventories have the opposite effect. Chuah subsequently [18] undertook a comprehensive study of the CFR problem, where he discussed the effects of various factors on the problem and proposed a gradual change in kanban levels to attain optimal cost. Further research led to the discovery that fixed frequency decisions could change inventory cost in the CFR problem and that frequency became a decision variable in JIT systems [1]. The multiple vehicle routing problem (MVRP) has been a focus of VRP research in addition to the CFR problem. Chan et al. [19] formulated a multiple depot, multiple vehicle, location routing problem with a robust location routing strategy to solve the MVRP. Gintner et al. [20] considered the MVRP with multiple depots, an issue that arose in public transport bus routes, and proposed a two-phase method to assign buses to cover a given set of trips to solve an optimal scheduling problem.Although the traditional CFR problem considers planning of the route and frequency, it does not consider the problem of multiple types of vehicles in the MVRP and simply uses vehicle load as a constraint. However, different vehicle types have different vehicle load capacities, and vehicle load can influence pickup frequency which in turn has a significant impact on inventory cost and transportation cost. Therefore, the decision regarding the choice of vehicle is significant. Some scholars have addressed this problem. Blanton and Wainwright [21] used a genetic algorithm to research the problem of scheduling vehicles of multiple types. Ahn and Rakha [22] investigated the effects of the choice of route on different types of vehicles using microscopic and macroscopic emission estimation tools. The results showed that, from a perspective that considers the environments as well as energy consumption, the shortest route is not always optimal. Cavalcante and Roorda [23] considered the choice of vehicle as a discrete variable to solve the discrete model problem. However, their work only considered transportation cost influenced by choice of vehicle without an analysis of the effect of the frequency plan on the inventory cost and transportation cost.The abovementioned method shows that the VRP in milk runs is now being considered in the context of JIT supply systems with ever-increasing constraints and practical orientation. However, no comprehensive study has yet been conducted to simultaneously consider route decision, frequency plan, and vehicle type choice. Therefore, this paper proposes the common frequency routing problem with vehicle type choice (CFR-VTC). We consider the dispatched vehicles, the pickup frequency, and routing as the objective of minimizing the cost of transportation, inventory, and dispatch. We propose a two-phase heuristic algorithm called two-phase tabu search (TS) with limited search scope (TP-TSLSS). The effectiveness of this algorithm is verified via numerical examples. In comparison with the simulated annealing algorithm (SA) and TS, our TP-TSLSS can significantly improve the efficiency of the search process and obtain more stable and accurate solutions in a relatively short time by generating an initial satisfactory solution and limiting the search scope. Moreover, we confirm that the multiple vehicle type mode incurs lower total cost than the same vehicle type mode.The remainder of this paper is organized as follows: in Section2, we describe the CFR-VTC problem and establish the mathematical model for it. The two-phase tabu search algorithm with limited search scope (TP-TSLSS) to solve this model is proposed in Section 3. Fifty-five numerical examples are employed in the experiment to demonstrate the effectiveness of the proposed algorithm and four transportation modes are compared to demonstrate efficacy of the CFR-VTC model in Section 4. Finally, we draw conclusions and suggest directions for future work in Section 5.


## 2. Formulation

### 2.1. Problem Analysis

A logistics network system is composed of a manufacturer and multiple suppliers. To ensure JIT production, the manufacturer uses multiple vehicle types for high-frequency pickups and small lot sizes in the milk run. With a production line that consumes parts linearly, vehicle type arrangement, route, frequency, and corresponding vehicle type planning are required to minimize total transportation, inventory, and dispatch costs.Figure1 shows the relationship between the parts inventory and time. In the same period, the frequency of inventory 2 in the figure is higher than that of inventory 1. However, the average value of inventory 2 is less than that of inventory 1. These results show that higher frequencies incur lower inventories, but an increase in frequency increases transportation cost. Additionally, different vehicle type arrangements will produce different dispatch costs for a stable supply; frequency is determined according to the carrying capacity of each vehicle type. Therefore, a trade-off point must exist for the arrangement of vehicle type, pickup frequency, and transportation route that can minimize the cost of transportation, inventory, and dispatch. Based on the above analysis, vehicle type, pickup frequency, and transportation route are selected as decision variables in our mathematic model.Figure 1
Relationship between parts inventory and time.We make the following assumptions.Assumption 1.
Each vehicle type possesses a different carrying capacity. Each route requires only one vehicle type to perform the relevant transportation task in the milk run on each day.Assumption 2.
There is a sufficient number of vehicles to complete the transportation task. The demand of the manufacturer by supplieri can only be transported by one truck in the completion of one task.


### 2.2. Notation

(1) Parameters
V
=
i
∣
i
=
0,1
,
2
,
…
,
n is set of suppliers and the manufacturer.
0 represents the manufacturer;1,…,n represent the suppliers.
U
=
u
∣
u
=
0,1
,
…
,
a is set of vehicle types.
η is cost of unit transportation distance.
β
i is cost of unit parts inventory of supplier i.
γ
u is the dispatch cost of vehicle type u.
c
i
j is the distance between supplier i and j.
q
u is carrying capacity of vehicle type u.
d
i is the pickup quantity for supplier i.
t
i is unloading time at supplier i.
t
i
j is travel time from supplier i to supplier j.
T is the maximum allowable cycle time of each route for a given frequency.(2) Decision Variables
m is the number of routes.
K
=
k
∣
k
=
0,1
,
2
,
…
,
m is set of routes:(1)xijku=1,if  typeuvehicle  passes  from  supplierito  supplierjin  routek0,otherwisei,j∈V,k∈K,u∈Uyiku=1,if  typeuvehicle  picks  up  at  supplieriin  routek0,otherwisei∈V,k∈K,u∈Ubku=1,if  typeuvehicle  is  dispatched  for  routek0,otherwisek∈K,u∈U.(3) Variables That Are Dependent on the Decision Variables
f
k: the frequency of route k depends on yiku. The definition is provided in (3).
p
i
k: the pickup quantity of supplier i for a given frequency depends on yiku. The definition is provided in (4).


### 2.3. The Proposed Mathematical Model for CFR-VTC

Consider(2)min⁡z=∑i∈V∑j∈V∑k∈K∑u∈Uη·fk·cij·xijku+∑i∈V∑k∈Kβi·pik+∑k∈K∑u∈Uγu·bku,where (3)fk=∑u∑iyiku·diqu,i∈V,k∈K,u∈U. “” means rounding up(4)pik=∑uyiku·difk,i∈V,k∈K,u∈Usubject to(5)∑u∑ixijku=∑uyjku,i,j∈V,k∈K,u∈U(6)∑u∑jxijku=∑uyiku,i,j∈V,k∈K,u∈U(7)∑u∑kyiku=1,i∈V,k∈K,u∈U(8)∑u∑i∈W∑j∈Wxijku≤W-1∀W⊆V,0∉W,W≥2,k∈K,u∈U(9)∑u∑iti·yiku+∑u∑i∑jtij·xijku≤T,i,j∈V,k∈K,u∈U(10)xijku=0or1,i,j∈V;k∈K;u∈U(11)yiku=0or1,i∈V;k∈K;u∈U(12)bku=0or1,k∈K;u∈U(13)sgn⁡∑iyiku=bku,i∈V,k∈K,u∈U.The first part of the objective function above addresses transportation cost, the second part addresses inventory cost, and the third computes the cost of dispatched vehicles. The transportation cost depends on the length of a given route as well as the frequency of the route. The inventory cost depends on the pickup quantity for the suppliers in each route and the route frequency. Each vehicle type has a corresponding dispatch cost; therefore, the third part of the function depends on the vehicle type choice.The following is a detailed description of constraints (5) to (13). Equations (5) to (7) ensure that each supplier distributes in only one route, and the path of each vehicle dispatched forms a loop. Constraint (8) ensures that the subloop is avoided. Equation (9) ensures that the cycle time of each route for one frequency is less than the maximum allowable cycle time T. The value of T depends on the real situation. Equations (10) to (12) indicate the scope of variables. Equation (13) ensures that each route has only one vehicle dispatched.


## 2.1. Problem Analysis

A logistics network system is composed of a manufacturer and multiple suppliers. To ensure JIT production, the manufacturer uses multiple vehicle types for high-frequency pickups and small lot sizes in the milk run. With a production line that consumes parts linearly, vehicle type arrangement, route, frequency, and corresponding vehicle type planning are required to minimize total transportation, inventory, and dispatch costs.Figure1 shows the relationship between the parts inventory and time. In the same period, the frequency of inventory 2 in the figure is higher than that of inventory 1. However, the average value of inventory 2 is less than that of inventory 1. These results show that higher frequencies incur lower inventories, but an increase in frequency increases transportation cost. Additionally, different vehicle type arrangements will produce different dispatch costs for a stable supply; frequency is determined according to the carrying capacity of each vehicle type. Therefore, a trade-off point must exist for the arrangement of vehicle type, pickup frequency, and transportation route that can minimize the cost of transportation, inventory, and dispatch. Based on the above analysis, vehicle type, pickup frequency, and transportation route are selected as decision variables in our mathematic model.Figure 1
Relationship between parts inventory and time.We make the following assumptions.Assumption 1.
Each vehicle type possesses a different carrying capacity. Each route requires only one vehicle type to perform the relevant transportation task in the milk run on each day.Assumption 2.
There is a sufficient number of vehicles to complete the transportation task. The demand of the manufacturer by supplieri can only be transported by one truck in the completion of one task.


## 2.2. Notation

(1) Parameters
V
=
i
∣
i
=
0,1
,
2
,
…
,
n is set of suppliers and the manufacturer.
0 represents the manufacturer;1,…,n represent the suppliers.
U
=
u
∣
u
=
0,1
,
…
,
a is set of vehicle types.
η is cost of unit transportation distance.
β
i is cost of unit parts inventory of supplier i.
γ
u is the dispatch cost of vehicle type u.
c
i
j is the distance between supplier i and j.
q
u is carrying capacity of vehicle type u.
d
i is the pickup quantity for supplier i.
t
i is unloading time at supplier i.
t
i
j is travel time from supplier i to supplier j.
T is the maximum allowable cycle time of each route for a given frequency.(2) Decision Variables
m is the number of routes.
K
=
k
∣
k
=
0,1
,
2
,
…
,
m is set of routes:(1)xijku=1,if  typeuvehicle  passes  from  supplierito  supplierjin  routek0,otherwisei,j∈V,k∈K,u∈Uyiku=1,if  typeuvehicle  picks  up  at  supplieriin  routek0,otherwisei∈V,k∈K,u∈Ubku=1,if  typeuvehicle  is  dispatched  for  routek0,otherwisek∈K,u∈U.(3) Variables That Are Dependent on the Decision Variables
f
k: the frequency of route k depends on yiku. The definition is provided in (3).
p
i
k: the pickup quantity of supplier i for a given frequency depends on yiku. The definition is provided in (4).


## 2.3. The Proposed Mathematical Model for CFR-VTC

Consider(2)min⁡z=∑i∈V∑j∈V∑k∈K∑u∈Uη·fk·cij·xijku+∑i∈V∑k∈Kβi·pik+∑k∈K∑u∈Uγu·bku,where (3)fk=∑u∑iyiku·diqu,i∈V,k∈K,u∈U. “” means rounding up(4)pik=∑uyiku·difk,i∈V,k∈K,u∈Usubject to(5)∑u∑ixijku=∑uyjku,i,j∈V,k∈K,u∈U(6)∑u∑jxijku=∑uyiku,i,j∈V,k∈K,u∈U(7)∑u∑kyiku=1,i∈V,k∈K,u∈U(8)∑u∑i∈W∑j∈Wxijku≤W-1∀W⊆V,0∉W,W≥2,k∈K,u∈U(9)∑u∑iti·yiku+∑u∑i∑jtij·xijku≤T,i,j∈V,k∈K,u∈U(10)xijku=0or1,i,j∈V;k∈K;u∈U(11)yiku=0or1,i∈V;k∈K;u∈U(12)bku=0or1,k∈K;u∈U(13)sgn⁡∑iyiku=bku,i∈V,k∈K,u∈U.The first part of the objective function above addresses transportation cost, the second part addresses inventory cost, and the third computes the cost of dispatched vehicles. The transportation cost depends on the length of a given route as well as the frequency of the route. The inventory cost depends on the pickup quantity for the suppliers in each route and the route frequency. Each vehicle type has a corresponding dispatch cost; therefore, the third part of the function depends on the vehicle type choice.The following is a detailed description of constraints (5) to (13). Equations (5) to (7) ensure that each supplier distributes in only one route, and the path of each vehicle dispatched forms a loop. Constraint (8) ensures that the subloop is avoided. Equation (9) ensures that the cycle time of each route for one frequency is less than the maximum allowable cycle time T. The value of T depends on the real situation. Equations (10) to (12) indicate the scope of variables. Equation (13) ensures that each route has only one vehicle dispatched.


## 3. The Design of the Two-Phase Heuristic Algorithm

CFR-VTC, as an extension of VRP, is NP hard. Therefore, we propose a two-phase heuristic algorithm called two-phase TS with limited search scope (TP-TSLSS) to solve the problem. In the first phase, an initial satisfactory solution is generated through a greedy heuristic algorithm to maximize the ratioz2 of the superior arc frequency to the inferior arc frequency. The following steps are used to generate the initial satisfactory solution: (1) the arcs are divided into the superior arc and the inferior arc according to the 80-20 rule (the superior arc and the inferior are described in Section 3.3); (2) the greedy heuristic algorithm is then employed to maximize the ratio z2 of the frequency of the superior arc to that of the inferior arc through the four neighborhood selection approaches described in Section 3.2; (3) the final solution of the first phase with the optimal z2max in the number of iterations is exported. In the second phase, the solution from the first phase is improved using the TS. To improve search efficiency, we limit the scope of the search in the algorithm to render the z2 of the candidate solution greater than the product of the z2max of the initial satisfactory solution and a coefficient α. Calculating the z2 value is simpler than calculating the objective function value z. Consequently, the TP-TSLSS greatly enhances the search efficiency by limiting the scope of z2 rather than directly using the objective function value z for the search.


### 3.1. Representation and Evaluation of Solution

We usex and v to represent the route and vehicle type, respectively, and these form an intact solution. We set the manufacturer code to 0 and represent the supplier numbers with i(i=1,2,…,n). We randomly select the suppliers and the manufacturer as a feasible route. We first randomly generate m routes (m∈[2,n-1], where m is an integer) and randomly distribute n suppliers into these routes. For example, solution x includes routes p10-2-3-5-0,  p20-7-1-9-8-0,…, until all suppliers have finished selecting. These routes constitute solution x. Second, we calculate the sum of ti and tij of each route p1,p2,…, of solution x, and we examine whether the total time of each route meets the time constraint of (9). If not, we randomly split the route whose total time exceeds T. For example, if the total time of route p2 exceeds T, the route p20-7-1-9-8-0 is randomly split as p210-7-1-0 and p220-9-8-0. Then, recalculate the total time of each new route, respectively, and repeat the same step until all the routes meet the time constraint of (9). In that way, the new solution adjusted constitutes the feasible solution. Finally, we randomly select the vehicle type uj(j=1,2,…,m) corresponding to each route of feasible solution and set [u1,u2,…,um] as solution v for each vehicle type. We use the objective function value to evaluate the quality of solutions.


### 3.2. Neighborhood Solution Selection Method

There are four methods to select a neighborhood solution:(1)
Randomly select a route and two suppliers on the route. Reverse the traversal order of the two suppliers. For instance, randomly select routep10-7-1-9-8-0. Randomly reverse the order of 7-1-9-8; p1 becomes 0-8-9-1-7-0.(2)
Randomly select two routes and a supplier from one of these two routes. Insert the chosen supplier randomly into the other route. If there is no supplier remaining in the route from where the supplier was taken, set the route to null. For instance, randomly select routesp10-2-3-5-0 and p20-7-1-9-8-0. Randomly select supplier 2 from p1 to be inserted into p2. p1becomes 0-3-5-0, and p2 becomes 0-7-1-2-9-8-0.(3)
Randomly select a route and a supplier from it. Add the supplier to another route. For instance, say there are five routesp1,p2,p3,p4, and p5. Randomly select route p10-2-3-5-0 and supplier 2, and move the latter into p6. p6 becomes 0-2-0.(4)
Randomly select a vehicle type of solutionv, allowing the vehicle type to change within the allowable vehicle type range. For instance, v=[2,2,3,4,5] becomes [2,2,3,5,5].


### 3.3. The First-Phase Algorithm

80-20 rule is originally proposed by the Italian economist Pareto, so it is also called Pareto’s law. Commonly it is stated that 20% of all causes bring about 80% of all effects [24–26]. In this paper, we first define an arc set A that records the distances between any two nodes, which contain the suppliers and the manufacturer. Then, according to the 80-20 rule, we define 20% of arcs in set A with the shortest length as the superior arc set A1 and the other arcs are defined as the inferior arc set A2. The lengths of all arcs in A1 are smaller than those of the arcs in A2; that is, a1<a2,a1∈A1,a2∈A2. The number of arcs in A2 is four times that of the number in A1. For a specific solution (x,v), we define fs as the sum of frequencies of all superior arcs and fi as the sum of frequencies of all inferior arcs. z2(x,v) is defined as fs/fi. The objective of the first-phase algorithm is to maximize z2. We set the number of iterations to M.The greedy heuristic algorithm is designed as shown in Algorithm1.Algorithm 1
Randomly generate an initial solution of (x0,v0), which meets the time constraint in  (9). Set the current solution xnow=x0, vnow=v0;
then, calculatez2(xnow,vnow). Set it = 1 (it records the current number of iterations).
Do while  it <M
If  rand < p1
Generate a neighborhood solution (x1, v1) through the first neighborhood method
Do while the time of the new solution (x1, v1) is not less than T
Regenerate a neighborhood solution (x1, v1) through the first neighborhood method
End  do
Ifz2(x1,v1)>z2(xnow,vnow)
xnow=x1
vnow=v1
End if
End if
Execute the other three neighborhood methods using a similar process.
it = it + 1
End  do
Output the initial satisfactory solution(x0′,v0′)=(xnow,vnow), and define z2max⁡=z2(x0′,v0′).


### 3.4. The Second-Phase Algorithm

Osman [27], Taillard et al. [28], and Alonso et al. [29] applied TS to solve the VRP and verified its effectiveness. In this paper, we limit the scope of search to ensure that the z2 value of every candidate solution is greater than α·z2max (where α is a coefficient slightly smaller in value than 1. We find that α∈[0.7,0.9] is an appropriate range). We call the second-phase algorithm TS with limited search scope (TSLSS). We set the number of iterations as NI.The TSLSS is shown in Algorithm2.Algorithm 2
Set the total number of iterations performed by the TS (NI). Set tabulist =⌀.  Set the initial solution to (x0′, v0′), which is generated in
the first phase.  Let the optimal solution bexbest=x0′, vbest=v0′ and the current solution xnow=x0′, vnow=v0′.
Letn=0 (record the number of current iterations). Set the number of candidate solutions (h) in the TS.
Do while  n < NI
Employ the four neighborhood methods to generatehcandidate solutions {(x1,v1),(x2,v2),…,(xh,vh)}  satisfying the time
constraint in  (9) with z2 values greater than α⋅z2max⁡. Calculate the objective function values of these solutions
{z(x1,v1),z(x2,v2),…,z(xh,vh)}.
Obtain (xopt, vopt), whose objective function value is optimal.
If  min{z(x1,v1),z(x2,v2),…,z(xh,vh)}<z(xbest,vbest)
(xnow,vnow)=(xopt,vopt)
Update the tabulist.
Else
Do while (xopt, vopt) is in the tabulist
Select the suboptimal solution from the candidate solutions and record it as (xopt, vopt).
End  do
(xnow,vnow)=(xopt,vopt)
Update the tabulist
End if
n=n+1
End do


### 3.5. The Process of Two-Phase TS with Limited Search Scope

We combine the first-phase algorithm in Section3.3 and the second-phase algorithm in Section 3.4 as the two-phase TS with limited search scope (TP-TSLSS).Figure2 shows the specific process.Figure 2
The TP-TSLSS process.


## 3.1. Representation and Evaluation of Solution

We usex and v to represent the route and vehicle type, respectively, and these form an intact solution. We set the manufacturer code to 0 and represent the supplier numbers with i(i=1,2,…,n). We randomly select the suppliers and the manufacturer as a feasible route. We first randomly generate m routes (m∈[2,n-1], where m is an integer) and randomly distribute n suppliers into these routes. For example, solution x includes routes p10-2-3-5-0,  p20-7-1-9-8-0,…, until all suppliers have finished selecting. These routes constitute solution x. Second, we calculate the sum of ti and tij of each route p1,p2,…, of solution x, and we examine whether the total time of each route meets the time constraint of (9). If not, we randomly split the route whose total time exceeds T. For example, if the total time of route p2 exceeds T, the route p20-7-1-9-8-0 is randomly split as p210-7-1-0 and p220-9-8-0. Then, recalculate the total time of each new route, respectively, and repeat the same step until all the routes meet the time constraint of (9). In that way, the new solution adjusted constitutes the feasible solution. Finally, we randomly select the vehicle type uj(j=1,2,…,m) corresponding to each route of feasible solution and set [u1,u2,…,um] as solution v for each vehicle type. We use the objective function value to evaluate the quality of solutions.


## 3.2. Neighborhood Solution Selection Method

There are four methods to select a neighborhood solution:(1)
Randomly select a route and two suppliers on the route. Reverse the traversal order of the two suppliers. For instance, randomly select routep10-7-1-9-8-0. Randomly reverse the order of 7-1-9-8; p1 becomes 0-8-9-1-7-0.(2)
Randomly select two routes and a supplier from one of these two routes. Insert the chosen supplier randomly into the other route. If there is no supplier remaining in the route from where the supplier was taken, set the route to null. For instance, randomly select routesp10-2-3-5-0 and p20-7-1-9-8-0. Randomly select supplier 2 from p1 to be inserted into p2. p1becomes 0-3-5-0, and p2 becomes 0-7-1-2-9-8-0.(3)
Randomly select a route and a supplier from it. Add the supplier to another route. For instance, say there are five routesp1,p2,p3,p4, and p5. Randomly select route p10-2-3-5-0 and supplier 2, and move the latter into p6. p6 becomes 0-2-0.(4)
Randomly select a vehicle type of solutionv, allowing the vehicle type to change within the allowable vehicle type range. For instance, v=[2,2,3,4,5] becomes [2,2,3,5,5].


## 3.3. The First-Phase Algorithm

80-20 rule is originally proposed by the Italian economist Pareto, so it is also called Pareto’s law. Commonly it is stated that 20% of all causes bring about 80% of all effects [24–26]. In this paper, we first define an arc set A that records the distances between any two nodes, which contain the suppliers and the manufacturer. Then, according to the 80-20 rule, we define 20% of arcs in set A with the shortest length as the superior arc set A1 and the other arcs are defined as the inferior arc set A2. The lengths of all arcs in A1 are smaller than those of the arcs in A2; that is, a1<a2,a1∈A1,a2∈A2. The number of arcs in A2 is four times that of the number in A1. For a specific solution (x,v), we define fs as the sum of frequencies of all superior arcs and fi as the sum of frequencies of all inferior arcs. z2(x,v) is defined as fs/fi. The objective of the first-phase algorithm is to maximize z2. We set the number of iterations to M.The greedy heuristic algorithm is designed as shown in Algorithm1.Algorithm 1
Randomly generate an initial solution of (x0,v0), which meets the time constraint in  (9). Set the current solution xnow=x0, vnow=v0;
then, calculatez2(xnow,vnow). Set it = 1 (it records the current number of iterations).
Do while  it <M
If  rand < p1
Generate a neighborhood solution (x1, v1) through the first neighborhood method
Do while the time of the new solution (x1, v1) is not less than T
Regenerate a neighborhood solution (x1, v1) through the first neighborhood method
End  do
Ifz2(x1,v1)>z2(xnow,vnow)
xnow=x1
vnow=v1
End if
End if
Execute the other three neighborhood methods using a similar process.
it = it + 1
End  do
Output the initial satisfactory solution(x0′,v0′)=(xnow,vnow), and define z2max⁡=z2(x0′,v0′).


## 3.4. The Second-Phase Algorithm

Osman [27], Taillard et al. [28], and Alonso et al. [29] applied TS to solve the VRP and verified its effectiveness. In this paper, we limit the scope of search to ensure that the z2 value of every candidate solution is greater than α·z2max (where α is a coefficient slightly smaller in value than 1. We find that α∈[0.7,0.9] is an appropriate range). We call the second-phase algorithm TS with limited search scope (TSLSS). We set the number of iterations as NI.The TSLSS is shown in Algorithm2.Algorithm 2
Set the total number of iterations performed by the TS (NI). Set tabulist =⌀.  Set the initial solution to (x0′, v0′), which is generated in
the first phase.  Let the optimal solution bexbest=x0′, vbest=v0′ and the current solution xnow=x0′, vnow=v0′.
Letn=0 (record the number of current iterations). Set the number of candidate solutions (h) in the TS.
Do while  n < NI
Employ the four neighborhood methods to generatehcandidate solutions {(x1,v1),(x2,v2),…,(xh,vh)}  satisfying the time
constraint in  (9) with z2 values greater than α⋅z2max⁡. Calculate the objective function values of these solutions
{z(x1,v1),z(x2,v2),…,z(xh,vh)}.
Obtain (xopt, vopt), whose objective function value is optimal.
If  min{z(x1,v1),z(x2,v2),…,z(xh,vh)}<z(xbest,vbest)
(xnow,vnow)=(xopt,vopt)
Update the tabulist.
Else
Do while (xopt, vopt) is in the tabulist
Select the suboptimal solution from the candidate solutions and record it as (xopt, vopt).
End  do
(xnow,vnow)=(xopt,vopt)
Update the tabulist
End if
n=n+1
End do


## 3.5. The Process of Two-Phase TS with Limited Search Scope

We combine the first-phase algorithm in Section3.3 and the second-phase algorithm in Section 3.4 as the two-phase TS with limited search scope (TP-TSLSS).Figure2 shows the specific process.Figure 2
The TP-TSLSS process.


## 4. Computational Experiments

In this section, our proposed two-phase heuristic algorithm is programmed in MATLAB R2011b to solve 55 small, medium, and large numerical problems. These computational experiments are run on an Acer 4820TG computer with an Intel i5 CPU (2.4 GHz) and 4 GB (3.47 available) of memory.


### 4.1. Experiments with Varying Supplier Size

We use a number of examples involving varying numbers of suppliers to validate our model and algorithm. When the number of suppliers is below 100, both the horizontal and vertical coordinates of suppliers are drawn from the uniform distribution in the range[0,100] using the randi function in MATLAB. Moreover, the coordinates of the suppliers are drawn from the uniform distribution in the range [0,200], when the number of suppliers is between 100 and 200. Finally, when the number of suppliers is between 200 and 300, the coordinates are drawn from the uniform distribution in the range [0,300]. The horizontal and vertical coordinates are represented by X and Y, respectively. The pickup quantity is randomly generated from the uniform distribution in the range [80,150], whereas the inventory cost of each supplier’s parts is subject to uniform distribution in the range [0,1]. The dispatch cost for each vehicle type is uniformly distributed from 100 to 400. We set the coordinates of the manufacture to (50,50). The algorithms described in Section 3 are implemented to test these numerical instances.


### 4.2. Algorithm Comparison

All of the 55 examples are generated using the method described in Section4.1. Similar problems have been solved in some of the literature [27–29] using TS algorithms, whereas other studies [30–32] have employed simulated annealing algorithms (SA) to solve them. Thus, to demonstrate the effectiveness of our model and the TP-TSLSS algorithm, each example in this study is run five times using TP-TSLSS, TS, and the SA. We then analyze and compare the experiment results obtained by the three algorithms. Stop criteria of SA and TS are defined as follows: SA stops when the temperature reaches a specific value according to the scale of each instance, and TS stops after a specific number of iterations are performed by TS. The number of iterations performed by TS is twice as TP-TSLSS.Table1 shows the results of the first phase of the two-phase algorithm. SN represents the number of suppliers. Column “z(x0,v0)¯” represents the average objective function value of the initial solution. “z(x′,v′)¯” represents the average objective function value of solutions obtained by the first-phase algorithm, and “t′¯” represents the time taken by the first-phase algorithm (the unit of time used is second). Following the execution of the first-phase algorithm, we find that the z2 value increases significantly within a short time. Consequently, the objective function value z can be improved significantly through the optimization of the first-phase algorithm within a short time.Table 1
Results obtained for the first-phase algorithm.
Instance
SN
z
(
x
′
,
v
′
)
¯
z
(
x
0
,
v
0
)
¯
z
2
(
x
′
,
v
′
)
¯
z
2
(
x
0
,
v
0
)
¯
t
′
¯ (s)
1
10
1839.1
3069.8
2
0.1667
0.32
2
10
2112.8
2582
2.4444
0.2308
0.31
3
10
2340.7
3816.8
2.3333
0.25
0.34
4
10
2175.7
2699.8
3.0435
0.5455
0.29
5
10
2490
2935.1
0.7059
0.2444
0.3
6
20
3158
5711.4
4.8182
0.2535
0.34
7
20
5064.5
6504.5
13.2381
0.3488
0.33
8
20
3600
5894.8
6.8667
0.1227
0.31
9
20
4349.6
5709.5
10.3529
0.3793
0.3
10
20
3227.8
5273.2
8.8571
0.1333
0.32
11
30
4292.2
7362.3
14.9
0.2202
0.36
12
30
6830.9
7640.1
7.1695
0.2818
0.35
13
30
6574
7296.3
55
0.0896
0.36
14
30
6402.4
6578.9
38.0147
0.23
0.37
15
30
4881.5
9252.2
24.5
0.4107
0.36
16
50
12523
14449
39.7647
0.244
0.36
17
50
11190
14649
18.5238
0.3557
0.35
18
50
11945
14672
129.7
0.3254
0.37
19
50
11919
12149
91
0.25
0.38
20
50
10930
17825
21.1429
0.2551
0.36
21
70
11737
17540
47.9091
0.2878
0.34
22
70
15508
18462
54.5625
0.1973
0.39
23
70
12484
20227
100.3333
0.2653
0.38
24
70
15890
20555
29
0.3419
0.37
25
70
14998
19819
17.2857
0.2903
0.35
26
90
18475
20596
45.0556
0.3529
0.38
27
90
19643
20328
47.2356
0.3427
0.4
28
90
18792
19873
46.729
0.2957
0.39
29
90
16268
20921
87.46
0.3268
0.38
30
90
20846
28943
36.67
0.3832
0.41
31
110
25081
46116
20.68
0.4035
0.83
32
110
30848
45420
18.57
0.4689
0.78
33
110
27515
42895
19.39
0.5073
0.79
34
110
28371
37652
18.5
0.4971
0.73
35
110
23627
39712
30.46
0.4837
0.78
36
150
38053
59964
18.0213
0.2545
0.92
37
150
36892
56279
32.63
0.2384
0.89
38
150
39604
62874
20.481
0.1953
0.93
39
150
40285
58164
16.306
0.2491
0.92
40
150
37681
63692
26.17
0.1902
0.9
41
200
86334
111050
5.5
0.2611
21.02
42
200
88274
113130
3.5806
0.2194
20.87
43
200
85145
113670
4.1758
0.2874
20.74
44
200
89732
112823
5.28
0.2391
21.23
45
200
87692
135713
4.825
0.261
20.84
46
250
111760
145270
3.9955
0.2601
40.29
47
250
111830
146220
4.316
0.2206
41.81
48
250
118713
149254
3.8421
0.2059
40.17
49
250
116915
153873
3.617
0.1973
39.82
50
250
120317
148148
3.291
0.2391
40.59
51
300
131980
169670
4.4194
0.2328
59.72
52
300
129468
163200
4.841
0.304
62.92
53
300
139642
161726
4.163
0.2845
58.17
54
300
130613
159212
4.0157
0.2174
61.83
55
300
129862
163944
4.9153
0.316
60.84Table2 shows the comparison results of the optimal objective function values of solutions obtained by the three algorithms. We choose examples 23, 31, 36, and 41 to show the convergence of TS and TP-TSLSS in Figures 3, 4, 5, and 6. The horizontal axis represents the number of iterations and the vertical axis represents the average objective function value. The column “average” represents the average objective function value of the optimal solution in each algorithm, “t-” represents the CPU time for each algorithm, and “σ” represents the standard deviation.Table 2
Results obtained for the three algorithms (SA, TS, and TP-TSLSS).
Instance
SN
SA
TS
TP-TSLSS
Average
t
SA
¯ (s)
σ
SA
Average
t
TS
¯ (s)
σ
TS
Average
t
TP-TSLSS
¯ (s)
σ
TP-TSLSS
1
10
1619.1
16.54
87.31
1386.2
32.47
0
1386.2
14.63
0
2
10
1574.3
15.26
76.31
1388
32.58
0
1388
14.19
0
3
10
1491
15.94
40.20
1421.9
32.82
0
1421.9
13.46
0
4
10
1582.5
16.31
72.18
1477.8
33.54
0
1477.8
13.22
0
5
10
1604.9
16.82
71.94
1487.5
31.96
40.31
1438.2
14.07
0
6
20
2568.9
17.27
82.56
2333.3
46.64
41.07
2354.9
16.86
32.18
7
20
2602.1
17.14
84.1
2669.8
45.61
56.2
2564.9
18.11
35.91
8
20
2492.8
17.93
86.46
2510.6
46.28
47.9
2421.7
17.65
34.5
9
20
2618.3
17.37
87.29
2765.6
46.19
60.1
2571.4
18.79
33.14
10
20
2310.5
17.49
84.39
2275.8
45.94
50.4
2242.3
18.14
36.2
11
30
3421.9
20.15
118.3
3527.8
54.92
76.60
3317.3
18.25
42.71
12
30
3679.8
21.25
110.38
3421.4
51.35
71.35
3308.2
18.19
41.6
13
30
3419.2
20.62
109.49
3507.3
52.84
70.58
3301.9
17.46
40.84
14
30
3519.3
20.48
111.64
3646.5
51.75
75.63
3396
18.22
40.23
15
30
3827.6
21.16
126.45
4214.5
54.28
79.37
3736.1
17.07
42.5
16
50
6425
24.85
410.72
5970.9
65.14
120.47
5405.1
21.16
55.57
17
50
6015.8
24.72
382.64
6110.1
66.28
110.62
5998
23.92
58.35
18
50
5702.8
24.93
341.83
5741.3
66.04
84.59
5627.6
22.16
51.94
19
50
5602.6
24.19
339.63
5870.8
65.92
92.46
5218.5
21.45
57.49
20
50
6092.4
25.62
329.74
6368.6
65.85
87.94
5686
22.81
53.85
21
70
8510.6
40.57
370.16
9079.7
78.49.
185.84
7615.3
21.36
117.03
22
70
8217.9
39.47
382.94
8820
76.32
172.94
7846.5
20.73
107.52
23
70
8246.2
39.29
411.74
8459.4
79.40
168.48
7386.5
24.91
104.73
24
70
8817.1
40.19
382.65
8990.5
78.35
171.82
8170.2
26.15
101.49
25
70
8109.6
40.95
347.93
8246.8
75.29
174.9
8051.6
24.91
113.52
26
90
10492
64.18
403.27
10737
106.73
272.32
9275.6
42.15
177.28
27
90
11600
64.29
410.46
12446
104.28
234.72
9079.4
38.81
171.36
28
90
10843.3
63.16
406.84
13647.5
105.82
258.19
8961.5
39.74
173.85
29
90
11692.3
62.85
409.58
11436.7
104.18
241.74
9235.9
40.39
174.82
30
90
10239.5
61.93
407.15
10869.1
108.37
229.64
9572.5
41.39
169.83
31
110
20487
89.72
779.67
21116.8
132.48
541.39
18016.1
57.97
338.29
32
110
18732.1
80.17
789.24
20429.2
130.82
527.64
11742.3
58.26
310.34
33
110
19035.5
85.26
835.92
19816.3
135.69
510.74
18535.4
54.94
340.93
34
110
19024.1
84.38
931.64
19252.5
131.94
509.86
17939.4
56.18
358.5
35
110
19348.2
85.25
764.92
21624.1
130.74
511.96
17310.6
57.69
379.29
36
150
30813
180.38
1348.4
27982
359.22
973.46
24094
150.06
835.35
37
150
29813
181.57
1246.7
30583
355.42
1004.8
24695
146.78
821.46
38
150
30847
182.62
1329.5
31764
356.86
986.83
28174
149.03
782.48
39
150
29139
180.37
1482.9
30269
357.42
904.74
27108
142.86
795.83
40
150
28742
179.84
1279.3
29264
352.91
981.34
26495
149.82
849.73
41
200
50954
240.44
3231.6
56340
412.49
2874.3
48246
191.51
1576.3
42
200
58192
238.59
3593.5
58861
414.65
2763.1
52407
200.18
1639.4
43
200
56102
241.68
3819.3
57928
416.72
2492.6
52887
194.27
1582.6
44
200
57149
240.27
3719.4
59218
418.39
2742.8
50961
187.37
1782.4
45
200
56291
239.39
3691.3
58316
419.62
2972.6
51620
194.72
1629.7
46
250
67194
380.54
4691.1
75935
630.28
3687.7
62171
320.06
2428.7
47
250
62719
379.37
4851.2
77501
635.42
3582.5
59954
310.78
2592.6
48
250
61461
378.46
4923.6
78519
638.87
3601.4
58917
314.03
2174.5
49
250
70192
377.35
4392.6
76110
640.41
3582.6
63625
317.86
2852.8
50
250
64192
373.42
4729.5
71926
639.93
3571.4
62128
326.82
2394.8
51
300
82155
520.69
7397.9
92167
912.49
5821.2
78166
463.51
4232.6
52
300
84124
510.53
7492.4
91135
894.65
5614.9
77640
447.18
3872.6
53
300
88307
517.49
7148.6
90114
906.72
5937.6
75783
482.27
3741.9
54
300
83491
512.83
7395.1
87910
918.39
5286.4
78629
419.37
4173.6
55
300
76529
509.59
7505.6
88412
909.62
5825.1
73210
426.72
4394.5Figure 3
Convergence graph of 70 suppliers.Figure 4
Convergence graph of 110 suppliers.Figure 5
Convergence graph of 150 suppliers.Figure 6
Convergence graph of 200 suppliers.We draw the following conclusions from a comparison of SA, TS, and TP-TSLSS, shown in Table2 and the convergence graphs (Figures 3, 4, 5, and 6):(1)
When the scales of the instances are small, the objective function values of the solutions obtained by SA are greater than those obtained by the other two algorithms. Therefore, the performance of SA is inferior to that of TS and TP-TSLSS in solving small-scale problems. By contrast, SA outperforms TS but is worse than TP-TSLSS when the number of suppliers exceeds 30. A comparison of TP-TSLSS and TS shows that the former yields solutions with smaller objective function values than the latter, except for a few small-scale instances. Moreover, Figures3, 4, 5, and 6 show that the number of iterations of TS is twice that of TP-TSLSS, but TS is inferior to TP-TSLSS. The performance of TP-TSLSS is better than that of the other two algorithms for most instances. The advantage of TP-TSLSS becomes more obvious as the scale of problems increases.(2)
The computational time of TP-TSLSS is the shortest, followed by SA and TS.(3)
Solutions obtained by TP-TSLSS have the smallest standard deviation. The comparison of standard deviation values proves that TP-TSLSS is the most robust, followed by TS with moderate robustness and SA with the least. Therefore, the two-phase algorithm (TP-TSLSS) can obtain better and more stable solutions in a relatively short time.To verify the effectiveness of the two-phase heuristic algorithm, we select an optimal solution for TP-TSLSS from an example involving 30 suppliers. Figure7 shows the path diagram.Figure 7
The path diagram of TP-TSLSS.In addition, we take thet-test of statistical analysis to compare difference in objective function values of the three algorithms. T test employs t distribution theory to infer the occurring probability of differences and to examine whether the differences are significant in two averages [33, 34]. This paper takes t-test to compare the differences in average objective function values of the three algorithms. First, we run 10 times for all 55 examples. We employ the SPSS 18 to take normality test for the objective function values of the 55 examples. We find that except the objective function values of the five numerical examples with the scale of 10 suppliers the rest of the numerical examples’ objective function values obtained by all the three algorithms have a good normality. Then we take two hypothesis tests at 0.05 confidence level conditions. The first hypothesis test is designed for SA and TP-TSLSS with the null hypothesis H: μ1≤μ0. The second hypothesis test is designed for TS and TP-TSLSS with the null hypothesis H: μ2≤μ0. Finally, we take independent sample t-test for the two hypotheses, respectively. The p values obtained are listed in Table 3.Table 3
p value of t-test.
Instance
SN
p value (SA and TP-TSLSS)
p value (TS and TP-TSLSS)
1
10
—
—
2
10
—
—
3
10
—
—
4
10
—
—
5
10
—
—
6
20
0.00031
0.16
7
20
0.0000028
0.18
8
20
0.000026
0.21
9
20
0.0043
0.12
10
20
0.000030
0.28
11
30
0.0051
0.0085
12
30
0.0032
0.027
13
30
0.0017
0.010
14
30
0.0072
0.0036
15
30
0.015
0.0029
16
50
0.0000022
0.0000025
17
50
0.0000063
0.00041
18
50
0.012
0.0061
19
50
0.0000027
0.0000035
20
50
0.0000012
0.0000018
21
70
0.00023
0.000010
22
70
0.00014
0.0000029
23
70
0.000092
0.0000012
24
70
0.000011
0.00000063
25
70
0.00047
0.000021
26
90
0.000016
0.000047
27
90
0.0000024
0.00000018
28
90
0.0000051
0.00000010
29
90
0.000030
0.000028
30
90
0.000029
0.000011
31
110
0.000089
0.000062
32
110
0.0000033
0.00000010
33
110
0.00032
0.00029
34
110
0.000092
0.000059
35
110
0.000075
0.0000033
36
150
0.0000013
0.0000014
37
150
0.0000014
0.0000020
38
150
0.000067
0.000051
39
150
0.00018
0.000053
40
150
0.00051
0.000047
41
200
0.000069
0.00000015
42
200
0.00000098
0.00000083
43
200
0.000074
0.0000019
44
200
0.0000023
0.00000012
45
200
0.000082
0.0000038
46
250
0.00000080
0.000000037
47
250
0.00000015
0.000000013
48
250
0.00000011
0.0000000098
49
250
0.0000064
0.00000013
50
250
0.000094
0.00000027
51
300
0.0000011
0.00000060
52
300
0.00000081
0.00000056
53
300
0.00000042
0.00000017
54
300
0.0000098
0.00000082
55
300
0.000020
0.000000041For the comparison of SA and TP-TSLSS, allp values of the examples with 20 or more suppliers are less than 0.05; therefore the null hypotheses of these examples are rejected, demonstrating that TP-TSLSS is better than SA.For the comparison of TS and TP-TSLSS, apart fromp values of examples of 20 suppliers scale (TS and TP-TSLSS) being greater than 0.05, p values of the remaining examples are less than 0.05; therefore the null hypotheses of these examples are rejected, demonstrating TP-TSLSS is better than TS when the suppliers are more than 30.


### 4.3. Results Analysis of the Different Modes

General studies of the routing problems did not consider the effect of vehicle type choice on transportation frequency and inventory cost, while we take the vehicle type choice into consideration. To examine this effect, we make a comparison of the general transportation mode (Mode1) with the CFR-VTC mode (Mode 2) for different scales of instances.Mode 1.
Only one vehicle type is allowed for the transportation task.Mode 2.
All vehicle types are available.Figure8 shows a comparison of the results for the two modes (Mode 1 includes three cases: only vehicle type 1, only vehicle type 2, and only vehicle type 3, in which the vehicle types are not required for decision; Mode 2 signifies that all the three vehicle types are involved in the decision). The load capacity of vehicle type 1 is the smallest whose value is 100, and that of vehicle type 2 is moderate whose value is 300, whereas the capacity of vehicle type 3 is the largest whose value is 600. The horizontal axis in Figure 8 represents the number of suppliers and the vertical axis shows the four types of cost—transportation cost, inventory cost, dispatch cost, and total cost. The four cost trends are shown in Figure 8.Figure 8
A comparison of the mode results.The transportation cost graph in Figure8 shows that the case involving only vehicle type 3 records the lowest transportation cost because a vehicle with a large load capacity requires less frequency to complete the pickup task. However, the inventory cost and dispatch cost in this case are the highest when the scales of instances are large because of low frequency and high load capacity. By contrast, the case involving only vehicle type 1 has the advantage of low inventory and dispatch costs, while the transportation cost in this case is the highest. However, the total cost of Mode 2 is lower than all the three cases in Mode 1. Therefore, we find that we can balance the cost of transportation, inventory, and dispatch by a reasonable choice of vehicle to acquire the optimal total cost. Mode 2 has greater cost advantages than Mode 1.


## 4.1. Experiments with Varying Supplier Size

We use a number of examples involving varying numbers of suppliers to validate our model and algorithm. When the number of suppliers is below 100, both the horizontal and vertical coordinates of suppliers are drawn from the uniform distribution in the range[0,100] using the randi function in MATLAB. Moreover, the coordinates of the suppliers are drawn from the uniform distribution in the range [0,200], when the number of suppliers is between 100 and 200. Finally, when the number of suppliers is between 200 and 300, the coordinates are drawn from the uniform distribution in the range [0,300]. The horizontal and vertical coordinates are represented by X and Y, respectively. The pickup quantity is randomly generated from the uniform distribution in the range [80,150], whereas the inventory cost of each supplier’s parts is subject to uniform distribution in the range [0,1]. The dispatch cost for each vehicle type is uniformly distributed from 100 to 400. We set the coordinates of the manufacture to (50,50). The algorithms described in Section 3 are implemented to test these numerical instances.


## 4.2. Algorithm Comparison

All of the 55 examples are generated using the method described in Section4.1. Similar problems have been solved in some of the literature [27–29] using TS algorithms, whereas other studies [30–32] have employed simulated annealing algorithms (SA) to solve them. Thus, to demonstrate the effectiveness of our model and the TP-TSLSS algorithm, each example in this study is run five times using TP-TSLSS, TS, and the SA. We then analyze and compare the experiment results obtained by the three algorithms. Stop criteria of SA and TS are defined as follows: SA stops when the temperature reaches a specific value according to the scale of each instance, and TS stops after a specific number of iterations are performed by TS. The number of iterations performed by TS is twice as TP-TSLSS.Table1 shows the results of the first phase of the two-phase algorithm. SN represents the number of suppliers. Column “z(x0,v0)¯” represents the average objective function value of the initial solution. “z(x′,v′)¯” represents the average objective function value of solutions obtained by the first-phase algorithm, and “t′¯” represents the time taken by the first-phase algorithm (the unit of time used is second). Following the execution of the first-phase algorithm, we find that the z2 value increases significantly within a short time. Consequently, the objective function value z can be improved significantly through the optimization of the first-phase algorithm within a short time.Table 1
Results obtained for the first-phase algorithm.
Instance
SN
z
(
x
′
,
v
′
)
¯
z
(
x
0
,
v
0
)
¯
z
2
(
x
′
,
v
′
)
¯
z
2
(
x
0
,
v
0
)
¯
t
′
¯ (s)
1
10
1839.1
3069.8
2
0.1667
0.32
2
10
2112.8
2582
2.4444
0.2308
0.31
3
10
2340.7
3816.8
2.3333
0.25
0.34
4
10
2175.7
2699.8
3.0435
0.5455
0.29
5
10
2490
2935.1
0.7059
0.2444
0.3
6
20
3158
5711.4
4.8182
0.2535
0.34
7
20
5064.5
6504.5
13.2381
0.3488
0.33
8
20
3600
5894.8
6.8667
0.1227
0.31
9
20
4349.6
5709.5
10.3529
0.3793
0.3
10
20
3227.8
5273.2
8.8571
0.1333
0.32
11
30
4292.2
7362.3
14.9
0.2202
0.36
12
30
6830.9
7640.1
7.1695
0.2818
0.35
13
30
6574
7296.3
55
0.0896
0.36
14
30
6402.4
6578.9
38.0147
0.23
0.37
15
30
4881.5
9252.2
24.5
0.4107
0.36
16
50
12523
14449
39.7647
0.244
0.36
17
50
11190
14649
18.5238
0.3557
0.35
18
50
11945
14672
129.7
0.3254
0.37
19
50
11919
12149
91
0.25
0.38
20
50
10930
17825
21.1429
0.2551
0.36
21
70
11737
17540
47.9091
0.2878
0.34
22
70
15508
18462
54.5625
0.1973
0.39
23
70
12484
20227
100.3333
0.2653
0.38
24
70
15890
20555
29
0.3419
0.37
25
70
14998
19819
17.2857
0.2903
0.35
26
90
18475
20596
45.0556
0.3529
0.38
27
90
19643
20328
47.2356
0.3427
0.4
28
90
18792
19873
46.729
0.2957
0.39
29
90
16268
20921
87.46
0.3268
0.38
30
90
20846
28943
36.67
0.3832
0.41
31
110
25081
46116
20.68
0.4035
0.83
32
110
30848
45420
18.57
0.4689
0.78
33
110
27515
42895
19.39
0.5073
0.79
34
110
28371
37652
18.5
0.4971
0.73
35
110
23627
39712
30.46
0.4837
0.78
36
150
38053
59964
18.0213
0.2545
0.92
37
150
36892
56279
32.63
0.2384
0.89
38
150
39604
62874
20.481
0.1953
0.93
39
150
40285
58164
16.306
0.2491
0.92
40
150
37681
63692
26.17
0.1902
0.9
41
200
86334
111050
5.5
0.2611
21.02
42
200
88274
113130
3.5806
0.2194
20.87
43
200
85145
113670
4.1758
0.2874
20.74
44
200
89732
112823
5.28
0.2391
21.23
45
200
87692
135713
4.825
0.261
20.84
46
250
111760
145270
3.9955
0.2601
40.29
47
250
111830
146220
4.316
0.2206
41.81
48
250
118713
149254
3.8421
0.2059
40.17
49
250
116915
153873
3.617
0.1973
39.82
50
250
120317
148148
3.291
0.2391
40.59
51
300
131980
169670
4.4194
0.2328
59.72
52
300
129468
163200
4.841
0.304
62.92
53
300
139642
161726
4.163
0.2845
58.17
54
300
130613
159212
4.0157
0.2174
61.83
55
300
129862
163944
4.9153
0.316
60.84Table2 shows the comparison results of the optimal objective function values of solutions obtained by the three algorithms. We choose examples 23, 31, 36, and 41 to show the convergence of TS and TP-TSLSS in Figures 3, 4, 5, and 6. The horizontal axis represents the number of iterations and the vertical axis represents the average objective function value. The column “average” represents the average objective function value of the optimal solution in each algorithm, “t-” represents the CPU time for each algorithm, and “σ” represents the standard deviation.Table 2
Results obtained for the three algorithms (SA, TS, and TP-TSLSS).
Instance
SN
SA
TS
TP-TSLSS
Average
t
SA
¯ (s)
σ
SA
Average
t
TS
¯ (s)
σ
TS
Average
t
TP-TSLSS
¯ (s)
σ
TP-TSLSS
1
10
1619.1
16.54
87.31
1386.2
32.47
0
1386.2
14.63
0
2
10
1574.3
15.26
76.31
1388
32.58
0
1388
14.19
0
3
10
1491
15.94
40.20
1421.9
32.82
0
1421.9
13.46
0
4
10
1582.5
16.31
72.18
1477.8
33.54
0
1477.8
13.22
0
5
10
1604.9
16.82
71.94
1487.5
31.96
40.31
1438.2
14.07
0
6
20
2568.9
17.27
82.56
2333.3
46.64
41.07
2354.9
16.86
32.18
7
20
2602.1
17.14
84.1
2669.8
45.61
56.2
2564.9
18.11
35.91
8
20
2492.8
17.93
86.46
2510.6
46.28
47.9
2421.7
17.65
34.5
9
20
2618.3
17.37
87.29
2765.6
46.19
60.1
2571.4
18.79
33.14
10
20
2310.5
17.49
84.39
2275.8
45.94
50.4
2242.3
18.14
36.2
11
30
3421.9
20.15
118.3
3527.8
54.92
76.60
3317.3
18.25
42.71
12
30
3679.8
21.25
110.38
3421.4
51.35
71.35
3308.2
18.19
41.6
13
30
3419.2
20.62
109.49
3507.3
52.84
70.58
3301.9
17.46
40.84
14
30
3519.3
20.48
111.64
3646.5
51.75
75.63
3396
18.22
40.23
15
30
3827.6
21.16
126.45
4214.5
54.28
79.37
3736.1
17.07
42.5
16
50
6425
24.85
410.72
5970.9
65.14
120.47
5405.1
21.16
55.57
17
50
6015.8
24.72
382.64
6110.1
66.28
110.62
5998
23.92
58.35
18
50
5702.8
24.93
341.83
5741.3
66.04
84.59
5627.6
22.16
51.94
19
50
5602.6
24.19
339.63
5870.8
65.92
92.46
5218.5
21.45
57.49
20
50
6092.4
25.62
329.74
6368.6
65.85
87.94
5686
22.81
53.85
21
70
8510.6
40.57
370.16
9079.7
78.49.
185.84
7615.3
21.36
117.03
22
70
8217.9
39.47
382.94
8820
76.32
172.94
7846.5
20.73
107.52
23
70
8246.2
39.29
411.74
8459.4
79.40
168.48
7386.5
24.91
104.73
24
70
8817.1
40.19
382.65
8990.5
78.35
171.82
8170.2
26.15
101.49
25
70
8109.6
40.95
347.93
8246.8
75.29
174.9
8051.6
24.91
113.52
26
90
10492
64.18
403.27
10737
106.73
272.32
9275.6
42.15
177.28
27
90
11600
64.29
410.46
12446
104.28
234.72
9079.4
38.81
171.36
28
90
10843.3
63.16
406.84
13647.5
105.82
258.19
8961.5
39.74
173.85
29
90
11692.3
62.85
409.58
11436.7
104.18
241.74
9235.9
40.39
174.82
30
90
10239.5
61.93
407.15
10869.1
108.37
229.64
9572.5
41.39
169.83
31
110
20487
89.72
779.67
21116.8
132.48
541.39
18016.1
57.97
338.29
32
110
18732.1
80.17
789.24
20429.2
130.82
527.64
11742.3
58.26
310.34
33
110
19035.5
85.26
835.92
19816.3
135.69
510.74
18535.4
54.94
340.93
34
110
19024.1
84.38
931.64
19252.5
131.94
509.86
17939.4
56.18
358.5
35
110
19348.2
85.25
764.92
21624.1
130.74
511.96
17310.6
57.69
379.29
36
150
30813
180.38
1348.4
27982
359.22
973.46
24094
150.06
835.35
37
150
29813
181.57
1246.7
30583
355.42
1004.8
24695
146.78
821.46
38
150
30847
182.62
1329.5
31764
356.86
986.83
28174
149.03
782.48
39
150
29139
180.37
1482.9
30269
357.42
904.74
27108
142.86
795.83
40
150
28742
179.84
1279.3
29264
352.91
981.34
26495
149.82
849.73
41
200
50954
240.44
3231.6
56340
412.49
2874.3
48246
191.51
1576.3
42
200
58192
238.59
3593.5
58861
414.65
2763.1
52407
200.18
1639.4
43
200
56102
241.68
3819.3
57928
416.72
2492.6
52887
194.27
1582.6
44
200
57149
240.27
3719.4
59218
418.39
2742.8
50961
187.37
1782.4
45
200
56291
239.39
3691.3
58316
419.62
2972.6
51620
194.72
1629.7
46
250
67194
380.54
4691.1
75935
630.28
3687.7
62171
320.06
2428.7
47
250
62719
379.37
4851.2
77501
635.42
3582.5
59954
310.78
2592.6
48
250
61461
378.46
4923.6
78519
638.87
3601.4
58917
314.03
2174.5
49
250
70192
377.35
4392.6
76110
640.41
3582.6
63625
317.86
2852.8
50
250
64192
373.42
4729.5
71926
639.93
3571.4
62128
326.82
2394.8
51
300
82155
520.69
7397.9
92167
912.49
5821.2
78166
463.51
4232.6
52
300
84124
510.53
7492.4
91135
894.65
5614.9
77640
447.18
3872.6
53
300
88307
517.49
7148.6
90114
906.72
5937.6
75783
482.27
3741.9
54
300
83491
512.83
7395.1
87910
918.39
5286.4
78629
419.37
4173.6
55
300
76529
509.59
7505.6
88412
909.62
5825.1
73210
426.72
4394.5Figure 3
Convergence graph of 70 suppliers.Figure 4
Convergence graph of 110 suppliers.Figure 5
Convergence graph of 150 suppliers.Figure 6
Convergence graph of 200 suppliers.We draw the following conclusions from a comparison of SA, TS, and TP-TSLSS, shown in Table2 and the convergence graphs (Figures 3, 4, 5, and 6):(1)
When the scales of the instances are small, the objective function values of the solutions obtained by SA are greater than those obtained by the other two algorithms. Therefore, the performance of SA is inferior to that of TS and TP-TSLSS in solving small-scale problems. By contrast, SA outperforms TS but is worse than TP-TSLSS when the number of suppliers exceeds 30. A comparison of TP-TSLSS and TS shows that the former yields solutions with smaller objective function values than the latter, except for a few small-scale instances. Moreover, Figures3, 4, 5, and 6 show that the number of iterations of TS is twice that of TP-TSLSS, but TS is inferior to TP-TSLSS. The performance of TP-TSLSS is better than that of the other two algorithms for most instances. The advantage of TP-TSLSS becomes more obvious as the scale of problems increases.(2)
The computational time of TP-TSLSS is the shortest, followed by SA and TS.(3)
Solutions obtained by TP-TSLSS have the smallest standard deviation. The comparison of standard deviation values proves that TP-TSLSS is the most robust, followed by TS with moderate robustness and SA with the least. Therefore, the two-phase algorithm (TP-TSLSS) can obtain better and more stable solutions in a relatively short time.To verify the effectiveness of the two-phase heuristic algorithm, we select an optimal solution for TP-TSLSS from an example involving 30 suppliers. Figure7 shows the path diagram.Figure 7
The path diagram of TP-TSLSS.In addition, we take thet-test of statistical analysis to compare difference in objective function values of the three algorithms. T test employs t distribution theory to infer the occurring probability of differences and to examine whether the differences are significant in two averages [33, 34]. This paper takes t-test to compare the differences in average objective function values of the three algorithms. First, we run 10 times for all 55 examples. We employ the SPSS 18 to take normality test for the objective function values of the 55 examples. We find that except the objective function values of the five numerical examples with the scale of 10 suppliers the rest of the numerical examples’ objective function values obtained by all the three algorithms have a good normality. Then we take two hypothesis tests at 0.05 confidence level conditions. The first hypothesis test is designed for SA and TP-TSLSS with the null hypothesis H: μ1≤μ0. The second hypothesis test is designed for TS and TP-TSLSS with the null hypothesis H: μ2≤μ0. Finally, we take independent sample t-test for the two hypotheses, respectively. The p values obtained are listed in Table 3.Table 3
p value of t-test.
Instance
SN
p value (SA and TP-TSLSS)
p value (TS and TP-TSLSS)
1
10
—
—
2
10
—
—
3
10
—
—
4
10
—
—
5
10
—
—
6
20
0.00031
0.16
7
20
0.0000028
0.18
8
20
0.000026
0.21
9
20
0.0043
0.12
10
20
0.000030
0.28
11
30
0.0051
0.0085
12
30
0.0032
0.027
13
30
0.0017
0.010
14
30
0.0072
0.0036
15
30
0.015
0.0029
16
50
0.0000022
0.0000025
17
50
0.0000063
0.00041
18
50
0.012
0.0061
19
50
0.0000027
0.0000035
20
50
0.0000012
0.0000018
21
70
0.00023
0.000010
22
70
0.00014
0.0000029
23
70
0.000092
0.0000012
24
70
0.000011
0.00000063
25
70
0.00047
0.000021
26
90
0.000016
0.000047
27
90
0.0000024
0.00000018
28
90
0.0000051
0.00000010
29
90
0.000030
0.000028
30
90
0.000029
0.000011
31
110
0.000089
0.000062
32
110
0.0000033
0.00000010
33
110
0.00032
0.00029
34
110
0.000092
0.000059
35
110
0.000075
0.0000033
36
150
0.0000013
0.0000014
37
150
0.0000014
0.0000020
38
150
0.000067
0.000051
39
150
0.00018
0.000053
40
150
0.00051
0.000047
41
200
0.000069
0.00000015
42
200
0.00000098
0.00000083
43
200
0.000074
0.0000019
44
200
0.0000023
0.00000012
45
200
0.000082
0.0000038
46
250
0.00000080
0.000000037
47
250
0.00000015
0.000000013
48
250
0.00000011
0.0000000098
49
250
0.0000064
0.00000013
50
250
0.000094
0.00000027
51
300
0.0000011
0.00000060
52
300
0.00000081
0.00000056
53
300
0.00000042
0.00000017
54
300
0.0000098
0.00000082
55
300
0.000020
0.000000041For the comparison of SA and TP-TSLSS, allp values of the examples with 20 or more suppliers are less than 0.05; therefore the null hypotheses of these examples are rejected, demonstrating that TP-TSLSS is better than SA.For the comparison of TS and TP-TSLSS, apart fromp values of examples of 20 suppliers scale (TS and TP-TSLSS) being greater than 0.05, p values of the remaining examples are less than 0.05; therefore the null hypotheses of these examples are rejected, demonstrating TP-TSLSS is better than TS when the suppliers are more than 30.


## 4.3. Results Analysis of the Different Modes

General studies of the routing problems did not consider the effect of vehicle type choice on transportation frequency and inventory cost, while we take the vehicle type choice into consideration. To examine this effect, we make a comparison of the general transportation mode (Mode1) with the CFR-VTC mode (Mode 2) for different scales of instances.Mode 1.
Only one vehicle type is allowed for the transportation task.Mode 2.
All vehicle types are available.Figure8 shows a comparison of the results for the two modes (Mode 1 includes three cases: only vehicle type 1, only vehicle type 2, and only vehicle type 3, in which the vehicle types are not required for decision; Mode 2 signifies that all the three vehicle types are involved in the decision). The load capacity of vehicle type 1 is the smallest whose value is 100, and that of vehicle type 2 is moderate whose value is 300, whereas the capacity of vehicle type 3 is the largest whose value is 600. The horizontal axis in Figure 8 represents the number of suppliers and the vertical axis shows the four types of cost—transportation cost, inventory cost, dispatch cost, and total cost. The four cost trends are shown in Figure 8.Figure 8
A comparison of the mode results.The transportation cost graph in Figure8 shows that the case involving only vehicle type 3 records the lowest transportation cost because a vehicle with a large load capacity requires less frequency to complete the pickup task. However, the inventory cost and dispatch cost in this case are the highest when the scales of instances are large because of low frequency and high load capacity. By contrast, the case involving only vehicle type 1 has the advantage of low inventory and dispatch costs, while the transportation cost in this case is the highest. However, the total cost of Mode 2 is lower than all the three cases in Mode 1. Therefore, we find that we can balance the cost of transportation, inventory, and dispatch by a reasonable choice of vehicle to acquire the optimal total cost. Mode 2 has greater cost advantages than Mode 1.


## 5. Conclusions

In this paper, we considered the common frequency routing problem with vehicle type choice in milk runs in logistical systems. We developed a mathematical model to describe the vehicle type choice, frequency planning, and vehicle route planning in the milk run system. To solve this model, we developed a two-phase TS algorithm with limited search scope (TP-TSLSS) that increased the search efficiency. The proposed TP-TSLSS algorithm was tested on 55 numerical examples with varying scales for verification and was compared with the TS and the SA methods. The results showed that our TP-TSLSS obtained better solutions in a shorter time and a more stable manner. By comparing different transportation modes, we concluded that consideration of the vehicle type choice could help save on cost of transportation, inventory, and dispatch.This was the first study to comprehensively consider frequency planning, vehicle type choice, and routing planning. Although the consideration of a variety of vehicle types significantly increased the complexity of the original CFR problem, the TP-TSLSS algorithm could reduce search time and improve accuracy through limiting search scope. In the future, we intended to extend this model to cases involving multiple manufacturers. The problem will be more complex, and the decision might depend on additional factors. Moreover, multiple depots could be considered from a practical perspective. Finally, we planned to develop even more efficient and accurate algorithms to solve these problems.



---

*Source: 404868-2015-10-18.xml*