# A Two-Stage Simulated Annealing Algorithm for the Many-to-Many Milk-Run Routing Problem with Pipeline Inventory Cost

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

---

## Abstract

In recent years, logistics systems with multiple suppliers and plants in neighboring regions have been flourishing worldwide. However, high logistics costs remain a problem for such systems due to lack of information sharing and cooperation. This paper proposes an extended mathematical model that minimizes transportation and pipeline inventory costs via the many-to-many Milk-run routing mode. Because the problem is NP hard, a two-stage heuristic algorithm is developed by comprehensively considering its characteristics. More specifically, an initial satisfactory solution is generated in the first stage through a greedy heuristic algorithm to minimize the total number of vehicle service nodes and the best insertion heuristic algorithm to determine each vehicle’s route. Then, a simulated annealing algorithm (SA) with limited search scope is used to improve the initial satisfactory solution. Thirty numerical examples are employed to test the proposed algorithms. The experiment results demonstrate the effectiveness of this algorithm. Further, the superiority of the many-to-many transportation mode over other modes is demonstrated via two case studies.


---

## Body

## 1. Introduction

Concomitant with the accelerated pace of globalization, economic competition has intensified. For better development, companies around the world have committed themselves to reducing costs, of which one of the most important constituents is logistics cost. A sign of this is the flourishing economic and technological development zones in China, where numerous enterprises are collocated in the same area. Other examples include the four major motor cities of the world (Toyota Motor City, Detroit Motor City, Turin Motor City, and Stuttgart Motor City). Accordingly, a multisupplier and multiplant logistics network system based on a factory’s spatial location or functional interaction is formed. Nevertheless, many systems of this kind have not fully developed into a real “cohesive cluster,” especially China’s automotive industry parks. Most automotive parks in China are simply “clustered” based on spatial location instead of the functional interaction among them. Therefore, inbound logistics transportation for this many-to-many (n-n) system remains a problem. Without doubt, inbound logistics transportation problems of this kind are more complex. However, proper planning can generate enormous benefits: it can not only strengthen cooperation among enterprises and contribute to their further development, but also reduce the cost of the entire logistics system and mitigate the pressure on logistics and circulation costs for society.Inbound logistics transportation methods can generally be classified into three types: direct shipping, Milk-run, and cross-docking. Berman and Wang [1] describe these three modes and their respective characteristics. Cross-docking is characterized by long transportation time and large pipeline inventory and is therefore suitable for long-distance transportation. In contrast, direct shipping and Milk-run have the advantage of short transportation time and small pipeline inventory, which make them suitable for short-range transportation. Considering the large number of suppliers and plants and the short distances between them, direct shipping or Milk-run is more applicable for an industrial park. This paper focuses on Milk-run route planning with multiple suppliers and plants within a certain adjacent area.The name “Milk-run” originated from the traditional system by which milk was sold in the West, in which the milkman simultaneously delivered full bottles of milk and collected the empty ones according to a predefined route. After completing a roundtrip, he returned with the empty bottles to the starting point. This method has subsequently been extensively utilized in miscellaneous industries and the auto manufacturing companies of the world and has developed into a popular way of picking up and delivering goods for multiple suppliers and plants using the same freight car [2]. Nemoto et al. [3] presented two Milk-run methods used by Toyota’s inbound logistics. The first is a many-to-one (n-1) mode, in which the car assembly manufacturer dispatches one truck at a specified time to visit various suppliers (i.e., parts suppliers) following a predefined route to collect parts or products and deliver them to the factories. The second is a one-to-many (1-n) mode, in which the parts are consolidated at the place of departure in such a way that they are fully loaded on the truck, transported, and unloaded at the various destinations. For both of these logistics modes, managers need to design and plan for specific routes, so as to minimize the overall logistics costs. Taking this rationale a step further, when all suppliers and plants share information with each other, then the many-to-many (n-n) Milk-run transportation described in [1] is achieved; that is, trucks pick up products from one or several suppliers and deliver them to one or several plants. Compared with Toyota’s n-1 and 1-n Milk-run system, this integrated n-n method can significantly reduce the transportation costs to a larger extent for the entire logistics system.The problem this paper focuses on is based on the classical vehicle routing problem (VRP). The generic VRP and many of its practical occurrences were studied intensively in the literature (Raff [4], Laporte and Nobert [5]). Then some extensive versions of VRP, such as vehicle routing with time windows (VRPTW) [6, 7] and multidepot vehicle routing problems (MDVRP) [8, 9], have been widely studied. These achievements of VRP have been already applied in the Milk-run system. Chuah and Yingling [10] and Jiang et al. [11] studied the common frequency routing problem in the Milk-run system and minimized transportation and inventory costs using column generation combined with the tabu search algorithm. Using genetic algorithm and a robust optimization approach, Sadjadi et al. [2] and Jafari-Eskandari et al. [12, 13] solved general and specific Milk-run issues with time windows and inventory uncertainty. Zhang and Li [14] analyzed the impact of the Milk-run method on supply chains, the conditions for using this method, and the changes in cost before and after adopting the Milk-run plan for a company. They also optimized the Milk-run mode by classifying containers into various carrying capacities based on the integral slice algorithm model. Yun et al. [15] analyzed the relationship between inventory and transportation, established an integrated transportation and inventory optimization model based on Milk-run, proposed genetic algorithm to solve the vehicle scheduling problem, applied a stepwise iterated algorithm to balance inventory and transportation costs, and eventually minimized the total costs for the entire Milk-run logistics system. Du et al. [16] investigated the parameter settings of a real-time vehicle-dispatching system for consolidating Milk-runs. Ma and Sun [17] employed a mutation ant colony algorithm to solve the Milk-run vehicle routing problem with the fastest completion time based on dynamic optimization. However, all of these studies targeted the multiple suppliers to one plant (n-1) problem. Effectually, none of them can be directly applied to Milk-run systems with multiple suppliers and plants. Given the status quo, this paper proposes a solution to the Milk-run routing problem involving multiple suppliers and plants (n-n).A number of researchers have already studied many-to-many transportation problems. The Pickup and Delivery Problem (PDP) is one of the research areas inn-n transportation issues. In general PDP, there is no restriction on the choice of pickup point to provide products for each specified delivery point, as long as the demand of the delivery point is satisfied. However, in practical terms, for inbound logistics systems, a plant’s demand for each supplier is always clearly determined because the plants will place specific orders with each supplier. Thus, the general PDP is not entirely consistent with the actual n-n inbound logistics systems. Furthermore, the set of variables employed to represent the commodities in PDP [18, 19] is not suitable for problems in practical terms. Further, because, in many manufacturing industries, such as auto manufacturing, the number of parts used in an individual assembly plant is more than 3,000 [20], the problem will become more complex if the set is defined based on the types of products. In light of the discussions above, PDP is not applicable for the type of logistics systems focused on in this paper. Therefore, an effective mathematical model with determined demand (hereafter, indicating the total demand of plant p for supplier i regarding one or more products, which will be described in more detail in Section 2) should be constructed.The cross-docking system is another transportation mode in inbound logistics involving multiple suppliers and plants. Several scholars have already studied vehicle routing problems in this system. For example, Berman and Wang [1] considered the problem of selecting the appropriate distribution strategy (direct shipping and cross-docking) for delivering a family of products from a set of suppliers to a set of plants so that the total transportation, pipeline inventory, and plant inventory costs were minimized. Lee et al. [21] optimized vehicle routing schedules by planning and scheduling vehicle routes during the pickup and delivery stages of cross-docking. Musa et al. [22] employed the ant colony algorithm to solve the n-n transportation problem of the cross-docking network. Miao et al. [23] employed a hybrid genetic algorithm to solve the multiple crossdocks problem with supplier and customer time windows. Mousavi and Tavakkoli-Moghaddam [24] minimized the total cost by simultaneously taking cross-dock’s location and route planning into consideration. The cross-docking transportation mode has been used extensively in long-distance transportation. However, it is not applicable for short-range transportation [2, 25]. In the cross-docking system, a vehicle must first pick up goods from suppliers and then transport them to the cross-dock, where inbound materials are sorted, consolidated, and stored until the outbound shipment is complete and ready to ship and finally deliver the goods to each plant. This increases the distance and the time needed for each transportation task and, thus, unnecessary pipeline inventory costs. In contrast, the Milk-run system has no transfer system similar to the cross-dock because all goods are directly delivered to the corresponding plant after they are picked up from each supplier. Consequently, all processes, including loading and unloading, are completed in the vehicle within a shorter time, which makes the Milk-run mode more suitable for short-distance transportation. Therefore, it is necessary to develop an effective method to solve the Milk-run routing problem involving multiple suppliers and plants in neighboring regions.From the discussion above, this paper contributes to solving the Milk-run routing problem involving multiple suppliers and plants in neighboring regions, that is, the many-to-many Milk-run routing problem (MM-MRP). In addition, we also consider the pipeline inventory cost caused by damage to products during the transportation process. It is necessary to consider the pipeline inventory costs because in some specific situations, some products, materials, and commodities, such as parts of some precision instruments, fragile parts, refrigerated foods, and some primary products, are more likely to be damaged with longer transportation time and more frequent loading and unloading times. On the basis of the above discussion, we firstly formulate the many-to-many Milk-run routing problem with pipeline inventory cost (MM-MRP-PIC) and extend the vehicle routing problem (VRP) model to MM-MRP-PIC model.Because VRP is NP hard, MM-MRP-PIC, an extension of VRP, is surely NP hard. Exact algorithms are difficult to find the optimal solutions within a reasonable time. Many scholars have developed heuristic algorithms to solve the VRP. Alfa et al. [26], Osman [27], and Van Breedam [28] applied simulated annealing to solve VRP and verified its effectiveness. Tabu search was developed by Osman [27] and Gendreau et al. [29] to solve VRP. Furthermore, genetic algorithms [30, 31], ant colony optimization [32–35], and particle swarm optimization [36, 37] are also developed to solve VRP successfully. However, some evolutionary algorithms, such as genetic algorithm, are inefficient to solve MM-MRP-PIC because the problem is so complex that they will spend much time to tackle illegal solutions generated during iterative process. In addition, ant colony can solve the problems whose objective function is only involved with transportation cost while it is difficult to tackle the pipeline inventory cost. Particle swarm optimization is also difficult to be designed for this problem because so various mutations should be considered that the PSO formulas cannot define these mutations together. However, simulated annealing can easily tackle illegal solutions and various mutations simultaneously considering the pipeline inventory cost. Therefore, we employ simulated annealing and make some improvements to solve MM-MRP-PIC. The algorithm is two-stage simulated annealing with limited search scope (TSSALSS). The effectiveness of this algorithm is verified via thirty numerical examples. Compared with traditional SA, TSSALSS not only significantly reduces the computing time, but also improves the quality of the solutions.The remainder of this paper is organized as follows. In Section2, we describe the problem and propose the MM-MRP-PIC model. Section 3 is devoted to the solution method (TSSALSS) for this model. Thirty numerical examples are presented to demonstrate the efficacy of our algorithm in Section 4. In Section 5, we present two case studies to demonstrate the validity of the many-to-many transportation mode. Finally, we conclude this paper in Section 6.


## 2. Formulation

### 2.1. Problem Description

This paper focuses on logistics systems comprising multiple suppliers and multiple plants in certain neighboring regions. In such a system, every supplier produces one or more parts, and each of the plants places corresponding orders to suppliers. The total demand of plantp for the products of supplier i (Dpi) can be indicated by the ratio of the products’ total volume to the truck capacity. The transportation process is organized as follows. A certain number of trucks, each of which has its own flow, are first scheduled for the shipping task. In each flow, the truck sets off from a plant, picks up products at a number of suppliers, delivers those products to corresponding plants, and finally returns to the departure point upon finishing the transportation task. This paper targets route planning for this kind of logistics systems, with the objective of minimizing transportation and pipeline inventory costs for the entire system.The following three assumptions are made to this problem.Assumption 1.
Product quantities cannot be split and each vehicle’s transportation frequency is one.Assumption 2.
There are enough vehicles at each plant to complete the transportation task.Assumption 3.
All vehicles employed for transportation are of the same type and have cargo capacityQ≥Dpi.


### 2.2. Notations

(1) Constants. Consider the following.S={1,2,…,n}: set of suppliers.P={n+1,n+2,…,n+m}: set of plants.N={1,2,…,n+m}(N=S∪P): set of nodes, including all suppliers and plants.V={1,2,…,s}: set of vehicles employed to complete the transportation task, with each vehicle corresponding to a path; that is, V is also the set of flows.Cij: transportation cost for shipping one truckload of products from node i to node j,i,j∈N.Dpi: the demand of plant p for the products of supplier i.vpi: pipeline inventory cost per unit time per unit of product Dpi.tij: time spent in transportation from node i to node j.Q: capacity of each vehicle.β: loading or unloading time per unit of product.(2) Decision Variables. Consider the following:(1)xijk=1,ifvehiclekpassesfromnodeitonodej0,otherwise,i,j∈N,k∈V,ypik=1,ifvehiclektransportsp’sdemandfori’sproduct0,otherwisep∈P,i∈S,k∈V.(3) Variables That Are Dependent on the Decision Variables. Consider the following:Lijk: pipeline inventory cost for products on vehicle k per unit time from node i to node j,i,j∈N,k∈V.Tik: time spent when vehicle k loads products at supplier i,i∈S.Tpk: time spent when vehicle k unloads products at plant p,p∈P.


### 2.3. The Proposed Mathematical Model for MM-MRP-PIC

Consider the following:(2)MinZ1=∑k∈V∑i∈N∑j∈NCij·xijk+∑k∈V∑j∈S∑i∈SLijkTjk+tij+∑k∈V∑i∈S∑p∈PLipktip+∑k∈V∑p∈P∑q∈PLpqkTpk+tpq,where (3)Tik=β∑p∈PypikDpi,i∈S,k∈VTpk=β∑i∈SypikDpi,p∈P,k∈Vsubject to(4)∑p∈P∑i∈SypikDpi≤Q,k∈V(5)∑k∈Vypik=sgn⁡Dpi,p∈P,i∈S(6)∑j∈Nxjik=sgn⁡∑p∈Pypik,i∈S,k∈V(7)∑j∈Nxijk=sgn⁡∑p∈Pypik,i∈S,k∈V(8)∑j∈Nxpjk=sgn⁡∑i∈Sypik,p∈P,k∈V(9)∑j∈Nxjpk=sgn⁡∑i∈Sypik,p∈P,k∈V(10)∑p∈P∑i∈Sxpik=1,k∈V(11)∑i∈S∑p∈Pxipk=1,k∈V(12)∑i,j∈Uxijk≤U-1,2≤U≤N-2,U⊂N,k∈V(13)Lpik=0,p∈P,i∈S,k∈V(14)Lijk=xijk∑l∈NLlik+∑p∈PDpiypikvpi,i∈S,j∈N,k∈V(15)Lpjk=xpik∑l∈NLlpk-∑i∈SDpiypikvpi,p∈P,j∈N,k∈V(16)xijk∈0,1,i,j∈N,k∈V(17)ypik∈0,1,p∈P,i∈S,k∈V.The first part of the objective function is direct transportation cost; the second, third, and fourth parts are pipeline inventory costs. The pipeline inventory costs depend on the road transportation time and residence time for freight loading or unloading at each point. An immediate idea for calculating the pipeline inventory costs for productsDpi is multiplying the pipeline inventory costs for products per unit time and the corresponding pipeline transportation time. The pipeline transportation time of Dpi is the period starting from the vehicle’s departure time at supplier i and ends at the arrival time at plant p. It includes transportation time on the way, loading time at various other suppliers, and unloading time at other plants in the same flow (we ignore the loading time of Dpi at supplier i and the unloading time at plant p because the pipeline inventory cost of Dpi caused by the loading time at supplier i and the unloading time at plant p are constant and independent of the decision variables). Although this idea is easy to understand, it is difficult to describe using mathematical models. Therefore, this paper divides the above process into three stages and calculates the respective cost for each part.(1)
Pickup stage: the pipeline inventory cost when vehiclek picks up products at supplier i and goes to the next supplier j can be represented by the product of the costs per unit time for the freight (Lijk) and the sum of the transportation time on road (tij) and the pickup time at supplier j (Tjk); that is, Lijk×(tij+Tjk). Tjk is proportional to the loading volume of vehicle k at supplier j, and the scale factor is β.(2)
Stage between the final pickup point and the first delivery point: the pipeline inventory cost is determined by the product of the costs per unit time for the freight (Lipk) and the time spent on transportation (tip); that is, Lipk×tip.(3)
Delivery stage: the pipeline inventory cost when vehiclek unloads products at plant p and goes to the next plant q can be represented by the product of the costs per unit time for the freight (Lpqk) and the sum of the transportation time on road (tpq) and the unloading time at p (Tpk); that is, Lpqk×(tpq+Tpk). Tpk is proportional to the unloading volume of vehicle k at plant p, and the scale factor is also assumed to be β.The following is a detailed description of constraints (4)–(17). Equation (4) signifies that the transportation task of vehicle k cannot exceed its capacity. Equation (5) means that as long as plant p has a demand for products from supplier i, a vehicle should be arranged to facilitate the transportation task. Equations (6)–(9) ensure that the path of each vehicle forms a loop. In classical VRP, there can only be one vehicle flow into and out of each node. In contrast, in this paper, more than one truck can go into and out of each node because multiple vehicles can be employed to execute different tasks at the same node. However, every vehicle that flows into a node must also leave that node; that is, the inflow and outflow for each vehicle at each node must be the same. In each flow, a vehicle must first pick up products from suppliers, then deliver goods to corresponding plants, and finally return to the departure point. To ensure that this process is conducted in an orderly manner, (10)-(11) restrict each vehicle to one single route from the plants to the suppliers, and vice versa. Equation (12) ensures that the solution does not contain any illegal isolated subtour. Equation (13) signifies that the pipeline inventory cost from the first plant to the starting supplier is zero. This is obvious because the vehicle sets off from a plant; therefore, it is empty during this period. Equation (14) represents the pipeline inventory cost per unit time for vehicle k from node i (i∈S) to node j (j∈N) in the pickup stage, which is the sum of the pipeline inventory cost per unit time for vehicle k before node i and that of the newly loading freight. Equation (15) is similar to (14); it represents the pipeline inventory cost per unit time for vehicle k from node p (p∈P) to node j (j∈N) in the delivery stage. Equations (16)-(17) signify the range of decision variables.


## 2.1. Problem Description

This paper focuses on logistics systems comprising multiple suppliers and multiple plants in certain neighboring regions. In such a system, every supplier produces one or more parts, and each of the plants places corresponding orders to suppliers. The total demand of plantp for the products of supplier i (Dpi) can be indicated by the ratio of the products’ total volume to the truck capacity. The transportation process is organized as follows. A certain number of trucks, each of which has its own flow, are first scheduled for the shipping task. In each flow, the truck sets off from a plant, picks up products at a number of suppliers, delivers those products to corresponding plants, and finally returns to the departure point upon finishing the transportation task. This paper targets route planning for this kind of logistics systems, with the objective of minimizing transportation and pipeline inventory costs for the entire system.The following three assumptions are made to this problem.Assumption 1.
Product quantities cannot be split and each vehicle’s transportation frequency is one.Assumption 2.
There are enough vehicles at each plant to complete the transportation task.Assumption 3.
All vehicles employed for transportation are of the same type and have cargo capacityQ≥Dpi.


## 2.2. Notations

(1) Constants. Consider the following.S={1,2,…,n}: set of suppliers.P={n+1,n+2,…,n+m}: set of plants.N={1,2,…,n+m}(N=S∪P): set of nodes, including all suppliers and plants.V={1,2,…,s}: set of vehicles employed to complete the transportation task, with each vehicle corresponding to a path; that is, V is also the set of flows.Cij: transportation cost for shipping one truckload of products from node i to node j,i,j∈N.Dpi: the demand of plant p for the products of supplier i.vpi: pipeline inventory cost per unit time per unit of product Dpi.tij: time spent in transportation from node i to node j.Q: capacity of each vehicle.β: loading or unloading time per unit of product.(2) Decision Variables. Consider the following:(1)xijk=1,ifvehiclekpassesfromnodeitonodej0,otherwise,i,j∈N,k∈V,ypik=1,ifvehiclektransportsp’sdemandfori’sproduct0,otherwisep∈P,i∈S,k∈V.(3) Variables That Are Dependent on the Decision Variables. Consider the following:Lijk: pipeline inventory cost for products on vehicle k per unit time from node i to node j,i,j∈N,k∈V.Tik: time spent when vehicle k loads products at supplier i,i∈S.Tpk: time spent when vehicle k unloads products at plant p,p∈P.


## 2.3. The Proposed Mathematical Model for MM-MRP-PIC

Consider the following:(2)MinZ1=∑k∈V∑i∈N∑j∈NCij·xijk+∑k∈V∑j∈S∑i∈SLijkTjk+tij+∑k∈V∑i∈S∑p∈PLipktip+∑k∈V∑p∈P∑q∈PLpqkTpk+tpq,where (3)Tik=β∑p∈PypikDpi,i∈S,k∈VTpk=β∑i∈SypikDpi,p∈P,k∈Vsubject to(4)∑p∈P∑i∈SypikDpi≤Q,k∈V(5)∑k∈Vypik=sgn⁡Dpi,p∈P,i∈S(6)∑j∈Nxjik=sgn⁡∑p∈Pypik,i∈S,k∈V(7)∑j∈Nxijk=sgn⁡∑p∈Pypik,i∈S,k∈V(8)∑j∈Nxpjk=sgn⁡∑i∈Sypik,p∈P,k∈V(9)∑j∈Nxjpk=sgn⁡∑i∈Sypik,p∈P,k∈V(10)∑p∈P∑i∈Sxpik=1,k∈V(11)∑i∈S∑p∈Pxipk=1,k∈V(12)∑i,j∈Uxijk≤U-1,2≤U≤N-2,U⊂N,k∈V(13)Lpik=0,p∈P,i∈S,k∈V(14)Lijk=xijk∑l∈NLlik+∑p∈PDpiypikvpi,i∈S,j∈N,k∈V(15)Lpjk=xpik∑l∈NLlpk-∑i∈SDpiypikvpi,p∈P,j∈N,k∈V(16)xijk∈0,1,i,j∈N,k∈V(17)ypik∈0,1,p∈P,i∈S,k∈V.The first part of the objective function is direct transportation cost; the second, third, and fourth parts are pipeline inventory costs. The pipeline inventory costs depend on the road transportation time and residence time for freight loading or unloading at each point. An immediate idea for calculating the pipeline inventory costs for productsDpi is multiplying the pipeline inventory costs for products per unit time and the corresponding pipeline transportation time. The pipeline transportation time of Dpi is the period starting from the vehicle’s departure time at supplier i and ends at the arrival time at plant p. It includes transportation time on the way, loading time at various other suppliers, and unloading time at other plants in the same flow (we ignore the loading time of Dpi at supplier i and the unloading time at plant p because the pipeline inventory cost of Dpi caused by the loading time at supplier i and the unloading time at plant p are constant and independent of the decision variables). Although this idea is easy to understand, it is difficult to describe using mathematical models. Therefore, this paper divides the above process into three stages and calculates the respective cost for each part.(1)
Pickup stage: the pipeline inventory cost when vehiclek picks up products at supplier i and goes to the next supplier j can be represented by the product of the costs per unit time for the freight (Lijk) and the sum of the transportation time on road (tij) and the pickup time at supplier j (Tjk); that is, Lijk×(tij+Tjk). Tjk is proportional to the loading volume of vehicle k at supplier j, and the scale factor is β.(2)
Stage between the final pickup point and the first delivery point: the pipeline inventory cost is determined by the product of the costs per unit time for the freight (Lipk) and the time spent on transportation (tip); that is, Lipk×tip.(3)
Delivery stage: the pipeline inventory cost when vehiclek unloads products at plant p and goes to the next plant q can be represented by the product of the costs per unit time for the freight (Lpqk) and the sum of the transportation time on road (tpq) and the unloading time at p (Tpk); that is, Lpqk×(tpq+Tpk). Tpk is proportional to the unloading volume of vehicle k at plant p, and the scale factor is also assumed to be β.The following is a detailed description of constraints (4)–(17). Equation (4) signifies that the transportation task of vehicle k cannot exceed its capacity. Equation (5) means that as long as plant p has a demand for products from supplier i, a vehicle should be arranged to facilitate the transportation task. Equations (6)–(9) ensure that the path of each vehicle forms a loop. In classical VRP, there can only be one vehicle flow into and out of each node. In contrast, in this paper, more than one truck can go into and out of each node because multiple vehicles can be employed to execute different tasks at the same node. However, every vehicle that flows into a node must also leave that node; that is, the inflow and outflow for each vehicle at each node must be the same. In each flow, a vehicle must first pick up products from suppliers, then deliver goods to corresponding plants, and finally return to the departure point. To ensure that this process is conducted in an orderly manner, (10)-(11) restrict each vehicle to one single route from the plants to the suppliers, and vice versa. Equation (12) ensures that the solution does not contain any illegal isolated subtour. Equation (13) signifies that the pipeline inventory cost from the first plant to the starting supplier is zero. This is obvious because the vehicle sets off from a plant; therefore, it is empty during this period. Equation (14) represents the pipeline inventory cost per unit time for vehicle k from node i (i∈S) to node j (j∈N) in the pickup stage, which is the sum of the pipeline inventory cost per unit time for vehicle k before node i and that of the newly loading freight. Equation (15) is similar to (14); it represents the pipeline inventory cost per unit time for vehicle k from node p (p∈P) to node j (j∈N) in the delivery stage. Equations (16)-(17) signify the range of decision variables.


## 3. Approach

As an extension of VRP, MM-MRP-PIC is NP hard. As a result, a two-stage simulated annealing algorithm with limited search scope (TSSALSS) is developed in this paper. Geng et al. [38] employed an adaptive simulated annealing algorithm with greedy search to solve the traveling salesman problem (TSP) and verified its effectiveness. In this paper, the process of limiting the search scope in TSSALSS is similar to the greedy search in [38], but we use a much simpler method to limit the search scope instead of greedy search which requires calculating the objective function frequently. A detailed solution process for TSSALSS is given below.Varanelli and Cohoon [39] stated that conventional SA can be initiated at a low temperature in an attempt to improve the heuristic solution. Thus, a faster heuristic (the first stage) is used to replace the SA actions occurring at the highest temperatures in the cooling schedule. In the first stage, an initial satisfactory solution is constructed in the following three steps: (1) Determine an initial feasible task allocation plan. (2) Minimize the total number of vehicle service nodes Z2 (calculated using formula (19) or (20)) by adjusting the task allocation for each vehicle. (3) Employ the best insertion heuristic algorithm to determine the route for each vehicle, which is achieved when the transportation cost is relatively optimal. In the second stage, the solution obtained in the first stage is improved by theSimulated Annealing (SA) algorithm. In order to improve the search efficiency, in this paper, the scope of the search in the algorithm is limited to a specific range. That is, in the search process, all feasible solutions are subject to the constraint that the total number of vehicle service nodes be not greater than γ·Z2min, in which Z2min is obtained in the first stage and γ is a number slightly larger than one; we find 1.1<γ<1.25 to be the most appropriate. Because it is much simpler to calculate ΔZ2 (the two neighbor solutions’ difference in the total number of vehicle service nodes, Z2) than to calculate ΔZ1 (the two neighbor solutions’ difference in the objective function value, Z1), the search process will be much more efficient with the limited search scope. This method can not only reduce the unwanted time caused by searching the overall feasible solution domain, but also prevent the global optimal solution being overlooked by appropriate expansion of the search scope adjusted by coefficient γ. The structure of the algorithm is shown in Figure 1. The following is a detailed description of the algorithm.Figure 1
Flow chart for the TSSALSS.


### 3.1. The First Stage of Algorithm

This stage consists chiefly of three steps: determination of an initial feasible task allocation plan, minimization of the total number of service nodes (Z2) accessed by all vehicles via the greedy heuristic algorithm, and determination of each vehicle’s route via the best insertion heuristic algorithm. The algorithms used in each step are as follows.Algorithm 1 (formation of an initial feasible task allocation plan).
We first determine the number of vehicles (s) needed:(18)s≤∑p∈P∑i∈SDpiα·Q,where α is a decimal slightly less than one. Note that if α is too close to one, we may miss the optimal solution, but if α is too small, it will cause unnecessary search times. We find that α∈[0.85,0.9] is most appropriate.
To ensure that (4) is satisfied, that is, no vehicle’s transportation task exceeds its loading capacity, the following method is used for task allocation: Randomly assign all transportation tasks to every vehicle.Do while at least one vehicle’s loading tasks are larger than the capacity Q.
Select the vehicles with the maximum and minimum loading tasks, respectively.Do while the transportation task of the maximum-loading- task vehicle is no longer greater than its loading capacity.
Assign one task of the maximum-loading-task vehicle to the minimum-loading-task vehicle.End doEnd doA feasible task allocation plan is thus obtained; this is indicated by the setY={y1,y2,…,ys}, where yk={ypik∣∀p,i}. Determine the service nodes, represented by nodesk, for each vehicle. For vehicle k, if ∑p∈Pypik≥1(∀i∈S,k∈V), then i∈nodesk; otherwise i∉nodesk. Similarly, if ∑i∈Sypik≥1(∀p∈P,k∈V), then p∈nodesk; otherwise p∉nodesk. The collection of nodes for all the vehicles is then represented by the set flow={nodes1,nodes2,…,nodess}.Algorithm 2 (greedy heuristic algorithm).
This algorithm minimizes the total number of service nodes (Z2) accessed by all vehicles, represented by(19)minZ2=∑k∈V∑i∈Ssgn⁡∑p∈Pypik+∑k∈V∑p∈Psgn⁡∑i∈Sypik;that is,(20)min∑k∈Vnodesk.The greedy heuristic algorithm is designed as follows:Set the current vehicle task allocation plan as the one obtained in Algorithm1 and calculate its total number of service nodes (Z2). Initialize Z2min=Z2.Forf=1:M
Randomly choose a vehiclek1 (k1∈V); one of its tasks is Dp0i0 (yp0i0k1=1), and another vehicle is k2 (k2∈V,k2≠k1).Do while the sum of Dp0i0 and k2’s transportation task exceeds the loading capacity ofk2
Repeat the random selection ofk1, k2, Dp0i0End doIfnodesk1+nodesk2 decreases or remains unchanged (ΔZ2≤0) when k2 is employed to finish transportation task Dp0i0
Accept the change andupdateyk1, yk2 (yp0i0k1=1→yp0i0k1=0,yp0i0k2=0→yp0i0k2=1) Z2min=Z2min+ΔZ2, and nodesk1, nodesk2 accordinglyEnd ifEnd forOutput the task allocation plan and the total number of service nodes (Z2min) accessed by all vehicles.Algorithm 3 (best insertion heuristic algorithm).
For the logistics system in this paper, each vehicle passes by suppliers to pick up products first and then delivers goods to plants. Considering this point, we modify the best insertion heuristic algorithm in [40] to ensure that (7)-(8) must be satisfied when inserting nodes. The pickup and delivery routes for all vehicles are represented by the sets pickup_route={PR1,PR2,…,PRs} and delivery_route={DR1,DR2,…,DRs}, in which PRk and DRk indicate vehicle k’s routes. The best insertion heuristic algorithm is designed as follows: InitializePRk∈pickup_route and DRk∈delivery_route as empty.Fork=1:s //construct the kth vehicle’s path
SelectCp0i0=max⁡(Cpi),p,i∈nodesk,p∈P,i∈S. Insert i0 into PRk and p0 into DRk and put i0, p0 into tabu list tabuk.Do whilenodesk>tabuk
Randomly select noded∈nodesk∖tabuk (d∈nodesk,d∉tabuk)Ifd≤n (n is the total number of suppliers)
Insertd into PRk, The insertion position of d is the place where the least increase in transportation cost is achievedElse
Insertd into DRk.End ifPutd into tabukEnd doEnd forOutput satisfactory solutionX0.


### 3.2. The Second Stage of Algorithm

For the SA in this paper, four neighbors are first designed, and the search scope is limited to a specific range by requiring that the total number of vehicle service nodes (∑k∈Vnodesk) be not greater than γ·Z2min. The four neighbors are performed at probabilities P1, P2, P3, and P4, respectively. A detailed description of the four neighbors with limited search scope designed in this paper is given below, and instances of the four neighbors are shown in Figure 2.Figure 2
Instances of the four neighbors.Algorithm 4 (the four neighbor structures with limited search scope).
Note that only Neighbors 3 and 4 have the possibility to generate new solutions whose total number of vehicle service nodes (Z2) will be greater than γ·Z2min. Neighbors 1 and 2 will not change the total number of vehicle service nodes because Neighbors 1 and 2 only change the vehicle routes and do not change the task allocation plan. Ifrand<P1 //Neighbor 1: block reverse mutation between suppliers (i.e., 2-opt)
Randomly choosePRk∈pickup_routeIf nodes in PRk are more than two
Randomly select two positions inPRk and reverse all the nodes between the two positionsElse
Give up the neighborhood operationEnd ifEnd ifIfrand<P2 //Neighbor 2: block reverse mutation between plants
Randomly chooseDRk∈delivery_routeIf nodes in DRk are more than two
Randomly select two positions inDRk and reverse all the nodes between the two positionsElse
Give up the neighborhood operationEnd ifEnd ifIfrand<P3 //Neighbor 3: a certain transportation task for a vehicle is handed over to another one
Randomly choose a vehiclek1 (k1∈V), one of its tasks Dp0i0 (Dp0i0≠0), and another vehicle k2 (k2∈V,k2≠k1).Do while the sum of Dp0i0 and k2’s transportation task exceeds the loading capacity of k2, or ∑k∈Vnodesk exceeds γ·Z2min after the handing over operation //limit the search scope
Repeat the random selection ofk1, k2, Dp0i0End doHandDp0i0 over to k2, update ypik1=1→ypik1=0, ypik2=0→ypik2=1,If∑p∈Pypik1=0//modify vehicle k1’s pickup route
Delete nodei from PRk1End ifIf∑i∈Sypik1=0 //modify vehicle k1’s delivery route
Delete nodep from DRk1End ifIfi∉nodesk2 before the neighborhood operation //modify vehicle k2’s pickup route
Insert nodei into any position of PRk2End ififp∉nodesk2 before the neighborhood operation //modify vehicle k2’s delivery route
Insert nodep into any position of DRk2End ifEnd ifIfP4<rand //Neighbor 4: exchange of transportation tasks among vehicles
Do while after exchanging the two tasks, either vehicle’s transportation task exceeds its loading capacity, or ∑k∈Vnodesk≥γ·Z2min  //limit the search scope
Repeat the random selection ofk1, k2, Dp1i1, Dp2i2End doModify the two vehicles’ pickup and delivery routes as in Neighbor 3End ifAlgorithm 5 (simulated annealing).
The temperature decreases at an exponential rate, that is,tcurrent=α·tcurrent, where α is the cooling rate and tcurrent is the current temperature. tstart is the initial temperature. The algorithm terminates when it reaches the final temperature (tend). The epoch length (Len) is constant.Initializetcurrent=tstart, Xcurrent=X0, which is obtained by the first stage of algorithm, Xbest=X0, Ebest=Z1(X0)Do whiletcurrent≥tend
Do whileit≤Len
Use Algorithm4 to generate a neighbor solution (X′) of XcurrentIfZ1(X′)≤Z1(Xcurrent)∥exp⁡(Z1(Xcurrent)-Z1(X′)/tcurrent)<rand
AcceptX′ as the current solution XcurrentIfZ1(X′)≤Z1(Xbest)
Xbest=X′, Ebest=Z1X′End ifEnd ifit=it+1End dotcurrent=α·tcurrentEnd do


## 3.1. The First Stage of Algorithm

This stage consists chiefly of three steps: determination of an initial feasible task allocation plan, minimization of the total number of service nodes (Z2) accessed by all vehicles via the greedy heuristic algorithm, and determination of each vehicle’s route via the best insertion heuristic algorithm. The algorithms used in each step are as follows.Algorithm 1 (formation of an initial feasible task allocation plan).
We first determine the number of vehicles (s) needed:(18)s≤∑p∈P∑i∈SDpiα·Q,where α is a decimal slightly less than one. Note that if α is too close to one, we may miss the optimal solution, but if α is too small, it will cause unnecessary search times. We find that α∈[0.85,0.9] is most appropriate.
To ensure that (4) is satisfied, that is, no vehicle’s transportation task exceeds its loading capacity, the following method is used for task allocation: Randomly assign all transportation tasks to every vehicle.Do while at least one vehicle’s loading tasks are larger than the capacity Q.
Select the vehicles with the maximum and minimum loading tasks, respectively.Do while the transportation task of the maximum-loading- task vehicle is no longer greater than its loading capacity.
Assign one task of the maximum-loading-task vehicle to the minimum-loading-task vehicle.End doEnd doA feasible task allocation plan is thus obtained; this is indicated by the setY={y1,y2,…,ys}, where yk={ypik∣∀p,i}. Determine the service nodes, represented by nodesk, for each vehicle. For vehicle k, if ∑p∈Pypik≥1(∀i∈S,k∈V), then i∈nodesk; otherwise i∉nodesk. Similarly, if ∑i∈Sypik≥1(∀p∈P,k∈V), then p∈nodesk; otherwise p∉nodesk. The collection of nodes for all the vehicles is then represented by the set flow={nodes1,nodes2,…,nodess}.Algorithm 2 (greedy heuristic algorithm).
This algorithm minimizes the total number of service nodes (Z2) accessed by all vehicles, represented by(19)minZ2=∑k∈V∑i∈Ssgn⁡∑p∈Pypik+∑k∈V∑p∈Psgn⁡∑i∈Sypik;that is,(20)min∑k∈Vnodesk.The greedy heuristic algorithm is designed as follows:Set the current vehicle task allocation plan as the one obtained in Algorithm1 and calculate its total number of service nodes (Z2). Initialize Z2min=Z2.Forf=1:M
Randomly choose a vehiclek1 (k1∈V); one of its tasks is Dp0i0 (yp0i0k1=1), and another vehicle is k2 (k2∈V,k2≠k1).Do while the sum of Dp0i0 and k2’s transportation task exceeds the loading capacity ofk2
Repeat the random selection ofk1, k2, Dp0i0End doIfnodesk1+nodesk2 decreases or remains unchanged (ΔZ2≤0) when k2 is employed to finish transportation task Dp0i0
Accept the change andupdateyk1, yk2 (yp0i0k1=1→yp0i0k1=0,yp0i0k2=0→yp0i0k2=1) Z2min=Z2min+ΔZ2, and nodesk1, nodesk2 accordinglyEnd ifEnd forOutput the task allocation plan and the total number of service nodes (Z2min) accessed by all vehicles.Algorithm 3 (best insertion heuristic algorithm).
For the logistics system in this paper, each vehicle passes by suppliers to pick up products first and then delivers goods to plants. Considering this point, we modify the best insertion heuristic algorithm in [40] to ensure that (7)-(8) must be satisfied when inserting nodes. The pickup and delivery routes for all vehicles are represented by the sets pickup_route={PR1,PR2,…,PRs} and delivery_route={DR1,DR2,…,DRs}, in which PRk and DRk indicate vehicle k’s routes. The best insertion heuristic algorithm is designed as follows: InitializePRk∈pickup_route and DRk∈delivery_route as empty.Fork=1:s //construct the kth vehicle’s path
SelectCp0i0=max⁡(Cpi),p,i∈nodesk,p∈P,i∈S. Insert i0 into PRk and p0 into DRk and put i0, p0 into tabu list tabuk.Do whilenodesk>tabuk
Randomly select noded∈nodesk∖tabuk (d∈nodesk,d∉tabuk)Ifd≤n (n is the total number of suppliers)
Insertd into PRk, The insertion position of d is the place where the least increase in transportation cost is achievedElse
Insertd into DRk.End ifPutd into tabukEnd doEnd forOutput satisfactory solutionX0.


## 3.2. The Second Stage of Algorithm

For the SA in this paper, four neighbors are first designed, and the search scope is limited to a specific range by requiring that the total number of vehicle service nodes (∑k∈Vnodesk) be not greater than γ·Z2min. The four neighbors are performed at probabilities P1, P2, P3, and P4, respectively. A detailed description of the four neighbors with limited search scope designed in this paper is given below, and instances of the four neighbors are shown in Figure 2.Figure 2
Instances of the four neighbors.Algorithm 4 (the four neighbor structures with limited search scope).
Note that only Neighbors 3 and 4 have the possibility to generate new solutions whose total number of vehicle service nodes (Z2) will be greater than γ·Z2min. Neighbors 1 and 2 will not change the total number of vehicle service nodes because Neighbors 1 and 2 only change the vehicle routes and do not change the task allocation plan. Ifrand<P1 //Neighbor 1: block reverse mutation between suppliers (i.e., 2-opt)
Randomly choosePRk∈pickup_routeIf nodes in PRk are more than two
Randomly select two positions inPRk and reverse all the nodes between the two positionsElse
Give up the neighborhood operationEnd ifEnd ifIfrand<P2 //Neighbor 2: block reverse mutation between plants
Randomly chooseDRk∈delivery_routeIf nodes in DRk are more than two
Randomly select two positions inDRk and reverse all the nodes between the two positionsElse
Give up the neighborhood operationEnd ifEnd ifIfrand<P3 //Neighbor 3: a certain transportation task for a vehicle is handed over to another one
Randomly choose a vehiclek1 (k1∈V), one of its tasks Dp0i0 (Dp0i0≠0), and another vehicle k2 (k2∈V,k2≠k1).Do while the sum of Dp0i0 and k2’s transportation task exceeds the loading capacity of k2, or ∑k∈Vnodesk exceeds γ·Z2min after the handing over operation //limit the search scope
Repeat the random selection ofk1, k2, Dp0i0End doHandDp0i0 over to k2, update ypik1=1→ypik1=0, ypik2=0→ypik2=1,If∑p∈Pypik1=0//modify vehicle k1’s pickup route
Delete nodei from PRk1End ifIf∑i∈Sypik1=0 //modify vehicle k1’s delivery route
Delete nodep from DRk1End ifIfi∉nodesk2 before the neighborhood operation //modify vehicle k2’s pickup route
Insert nodei into any position of PRk2End ififp∉nodesk2 before the neighborhood operation //modify vehicle k2’s delivery route
Insert nodep into any position of DRk2End ifEnd ifIfP4<rand //Neighbor 4: exchange of transportation tasks among vehicles
Do while after exchanging the two tasks, either vehicle’s transportation task exceeds its loading capacity, or ∑k∈Vnodesk≥γ·Z2min  //limit the search scope
Repeat the random selection ofk1, k2, Dp1i1, Dp2i2End doModify the two vehicles’ pickup and delivery routes as in Neighbor 3End ifAlgorithm 5 (simulated annealing).
The temperature decreases at an exponential rate, that is,tcurrent=α·tcurrent, where α is the cooling rate and tcurrent is the current temperature. tstart is the initial temperature. The algorithm terminates when it reaches the final temperature (tend). The epoch length (Len) is constant.Initializetcurrent=tstart, Xcurrent=X0, which is obtained by the first stage of algorithm, Xbest=X0, Ebest=Z1(X0)Do whiletcurrent≥tend
Do whileit≤Len
Use Algorithm4 to generate a neighbor solution (X′) of XcurrentIfZ1(X′)≤Z1(Xcurrent)∥exp⁡(Z1(Xcurrent)-Z1(X′)/tcurrent)<rand
AcceptX′ as the current solution XcurrentIfZ1(X′)≤Z1(Xbest)
Xbest=X′, Ebest=Z1X′End ifEnd ifit=it+1End dotcurrent=α·tcurrentEnd do


## 4. Computational Experiments

In this section, thirty numerical examples are designed to test the performance of the proposed algorithms. We compared our algorithms with the standard simulated annealing (SA) and two-stage simulated annealing (TSSA) to demonstrate the effectiveness of the improved algorithm. The same initial solutions are set for the three algorithms. We apply this benchmark method in order to compare the three algorithms fairly. We can not only obtain the gaps among the three algorithms, but also the number of cases for which TSSALSS performs worse than TSSA and SA with the same initial solutions. In order to verify the robustness and validity of the proposed algorithms, we test our algorithms for three groups of instances: small-scale problems with 3 to 5 plants and 10 to 20 suppliers, medium-scale problems with 10 to 20 plants and 20 to 50 suppliers, and large-scale problems with 10 to 30 plants and 100 to 300 suppliers. Simultaneously, it is obvious to observe the trend of gaps among the three algorithms with the scales of instances increasing. Furthermore, in order to test whether the proposed algorithms can adapt to numerical examples with different characteristics, the data are generated by the following methods: the coordinates of suppliers and plants are drawn from uniform distributions and the ranges of these distributions increase with the scales of instances increasing; some of the transportation tasksDpi are set to be zeros and the nonzero tasks Dpi are drawn from the uniform distributions from 5 to 15 or from 5 to 30; vpi (pipeline inventory cost per unit time per unit of product Dpi) are all drawn from the same uniform distribution from 0 to 0.1; the loading or unloading time per unit of product (β) is also drawn from the uniform distribution from 0.015 to 0.025. The detailed setting of data and parameters in the algorithms are listed in Appendix A. The algorithms were programmed in MATLAB 2012b. All of the experiments were carried out on a Lenovo T4900v i5 3470 computer with 4 GB (3.47 available) RAM and a 3.20 GHz CPU. Computation times are in seconds.We first generated ten initial feasible task allocation plans using Algorithm1 and determined ten initial feasible solutions X0={X10,X20,…,X100} by randomly sequencing the service nodes for each vehicle. Then, ten initial satisfactory solutions X={X1,X2,…,X10} were obtained using Algorithms 2 and 3 based on the same ten initial feasible task allocation plans. The results are presented in Table 1, where Z2¯ indicates the average number of nodes that all the vehicles should serve in the ten initial feasible solutions X0={X10,X20,…,X100}, Z2min¯ indicates the average number of nodes that all the vehicles should serve in ten initial satisfactory solutions X={X1,X2,…,X10}, Z1(X)¯ is the average value of the objective function, T¯ represents the average CPU time for the first stage of algorithm, and M indicates the total number of iterations. From Table 1, it can be seen that significant improvements in the value of the objective function were obtained via the first stage of algorithm within a short time.Table 1
Results obtained for the first stage of algorithm.
Instance
Z
2
¯
Z
2
min
⁡
¯
Z
1
(
X
0
)
¯
Z
1
(
X
)
¯
T
¯
(
s
)
M
S1
37.1
21.3
650.756
415.469
0.0038
7000
S2
42.1
28.4
997.775
589.488
0.0053
7000
S3
75.2
41.3
1717.88
827.196
0.011
25000
S4
87.9
59.1
2160.15
1271.33
0.015
25000
S5
70.8
36.8
1555.01
721.974
0.0089
20000
S6
78.5
46.0
1402.56
800.643
0.012
20000
S7
130.6
70.2
2938.82
1310.18
0.023
60000
S8
162.0
90.7
3845.58
1855.90
0.035
60000
M1
304.2
142.0
6571.25
2610.43
0.085
200000
M2
347.4
184.4
9322.98
4429.06
0.11
200000
M3
459.7
213.1
14318.9
5404.58
0.21
500000
M4
531.6
273.0
14118.2
6618.70
0.27
500000
M5
1026.1
425.5
31096.5
11244.5
0.60
1500000
M6
1101.0
546.2
34617.8
15927.7
0.81
1500000
M7
867.0
374.0
27949.6
9759.69
0.39
1000000
M8
928.9
467.9
30983.4
13920.5
0.54
1000000
M9
1710.2
710.3
52170.3
18493.5
1.85
3000000
M10
1856.0
916.0
57091.2
25443.4
1.66
3000000
L1
899.1
472.1
36851.0
17388.2
0.60
1000000
L2
1795.3
913.2
73359.4
33653.0
1.73
3000000
L3
1799.8
946.8
78274.9
37801.9
1.86
3000000
L4
1606.1
1857.2
174523.0
77891.9
5.41
9000000
L5
1877.7
947.7
68970.8
31358.7
1.78
3000000
L6
3770.7
1831.3
148326.3
65545.2
5.17
9000000
L7
3768.3
1857.3
144844.0
64654.8
5.54
9000000
L8
7564.6
3719.6
287945.5
129019.1
17.13
30000000
L9
1926.6
963.2
75724.3
33979.2
1.70
3000000
L10
3823.3
1860.0
149613.3
65616.0
5.07
9000000
L11
3830.4
1927.5
208620.9
91394.0
5.31
9000000
L12
7680.7
3746.5
416486.8
181603.0
16.84
30000000Table2 presents the results of the following three algorithms: (1) two-stage simulated annealing algorithm with limited search scope (TSSALSS), (2) two-stage simulated annealing algorithm (TSSA), and (3) simulated annealing (SA) algorithm. We utilized SA to improve the initial feasible solutions X0={X10,X20,…,X100} and obtained X3∗={X13∗,X23∗,…,X103∗}. TSSALSS and TSSA were employed to generate X1∗={X11∗,X21∗,…,X101∗} and X2∗={X12∗,X22∗,…,X102∗} based on the initial satisfactory solutions X={X1,X2,…,X10}, which were obtained via the first stage of algorithm from the initial feasible task allocation plans. S(X2∗) is defined as follows: S(X2∗)={Xi2∗∣Z1Xi2∗<Z1Xi1∗,Xi1∗∈X1∗}, where Z1(X) indicates the objective function value of solution X. S(X3∗) is similarly defined. Column SX2∗(SX3∗) is the number of solutions in S(X2∗)(S(X3∗)), which, respectively, represents the number of cases for which TSSALSS performs worse than TSSA and SA. Column gap 1 (gap 2) is the gap between average Z1(Xi1∗) and average Z1(Xi2∗)(Z1(Xi3∗)). The columns “average,” “minimum,” “std,” “T,” and “NI,” respectively, represent the objective function value of the best solution among the ten trials, the average objective function value of the ten trials, the standard deviation of the ten solutions, the average time spent by the two algorithms, and the number of iterations performed by the two algorithms.Table 2
Results obtained for the three algorithms (TSSALSS, TSSA, and SA).
Instances
TSSALSS
TSSA
SA
Average
Minimum
Std
T
NI
Average
Minimum
Std
T
|
S
(
X
2
∗
)
|
Gap 1
NI
Average
Minimum
Std
T
|
S
(
X
3
∗
)
|
Gap 1
NI
S1
304.901
297.15
4.41
1.73
30988
302.645
295.69
3.65
1.96
5
−0.74%
37414
303.769
298.92
2.46
3.12
6
−0.37%
56618
S2
456.536
449.86
3.98
1.88
30988
452.370
444.21
4.83
2.04
6
−0.91%
37414
455.791
448.72
4.04
3.16
7
−0.16%
56618
S3
602.408
585.28
8.14
5.67
72306
602.332
587.65
9.84
6.83
6
−0.01%
99770
598.319
584.15
7.82
8.89
4
−0.68%
152817
S4
879.412
848.56
18.20
5.70
72306
876.568
837.68
27.13
7.29
5
−0.32%
99770
875.906
846.73
19.89
8.61
6
−0.40%
152817
S5
501.301
493.02
5.18
3.59
51647
503.956
490.13
8.55
4.11
3
0.53%
74828
502.876
493.48
5.24
6.02
3
0.31%
114613
S6
629.063
612.58
10.24
3.43
51647
630.897
620.04
7.24
4.34
4
0.29%
74828
629.803
614.61
9.64
6.14
3
0.12%
114613
S7
934.877
902.45
19.29
4.05
92965
929.226
893.19
24.96
7.62
5
−0.60%
149655
936.248
906.73
20.39
13.16
2
0.15%
235956
S8
1239.98
1192.3
27.04
4.47
92965
1244.04
1217.4
19.98
7.31
3
0.33%
149655
1247.59
1216.6
20.09
13.82
3
0.61%
235956
M1
1788.32
1748.3
33.93
9.44
206589
1790.89
1768.9
12.69
19.35
1
0.14%
374138
1789.84
1760.8
15.76
34.67
2
0.08%
602942
M2
3056.33
2961.1
46.96
10.38
206589
3099.02
3067.3
20.22
19.58
0
1.40%
374138
3100.79
3047.8
36.31
36.12
1
1.45%
619767
M3
3362.66
3316.6
34.38
23.01
413178
3423.01
3371.1
32.18
47.55
0
1.79%
872990
3403.08
3351.0
37.11
83.74
1
1.20%
1446124
M4
4304.45
4244.3
42.61
24.32
413178
4400.65
4343.8
43.99
46.28
0
2.23%
872990
4393.82
4311.3
51.87
84.20
0
2.08%
1459131
M5
6668.53
6596.3
48.92
73.86
1239535
6757.88
6669.6
55.29
134.65
0
1.34%
2244830
6803.68
6673.4
80.34
233.27
0
2.03%
3810673
M6
10201.1
10107
60.24
76.35
1239535
10415.9
10338
86.07
135.77
0
2.11%
2244830
10426.8
10302
93.99
239.09
0
2.21%
3796910
M7
5933.78
5880.7
59.37
74.02
1239535
6025.07
5928.7
60.39
127.23
1
1.54%
2244830
6034.00
5941.9
62.11
226.57
0
1.69%
3796910
M8
9039.82
8905.8
84.22
79.19
1239535
9192.46
9040.1
122.74
137.51
0
1.69%
2244830
9212.56
9086.4
100.60
256.14
0
1.91%
3860891
M9
10582.9
10396
132.05
160.04
2582364
10790.7
10561
173.01
320.22
0
1.96%
4988511
10805.6
10536
222.69
570.52
0
2.10%
8804100
M10
13712.1
13523
137.71
184.77
2788953
14000.7
13810
151.68
395.58
0
2.10%
5237937
14005.5
13834
133.30
672.36
0
2.14%
9244305
L1
9829.9
9766.5
58.41
83.20
1342829
9972.3
9813.0
87.24
151.05
0
1.45%
2494256
10033.1
9917.0
88.39
275.78
0
2.07%
4327041
L2
18127.8
17930
119.90
203.29
2995542
18621.3
18342
186.38
411.13
0
2.72%
5487363
18809.6
18477
273.08
728.70
0
3.76%
9785545
L3
25307.3
25023
170.72
358.27
3098836
26132.2
25663
283.87
722.54
0
3.26%
5612075
26134.0
25733
209.27
1103.96
0
3.27%
10100378
L4
44704.2
44276
319.65
1532.21
8712249
46402.4
45235
730.92
3212.35
0
3.80%
16941776
46353.9
45637
503.28
5932.82
0
3.69%
30370067
L5
20877.3
20651
143.83
224.84
3098836
21468.2
21039
308.38
423.63
0
2.83%
5612075
21518.3
21356
273.78
743.19
0
3.07%
10100378
L6
39325.9
38912
251.02
824.93
8712249
40496.4
40076
296.43
1564.54
0
2.98%
16941776
40566.8
40233
219.51
2688.31
0
3.16%
30370067
L7
38951.6
38327
389.55
1675.24
10759079
40706.9
39124
1266.30
3973.80
0
4.51%
19527397
40560.1
39646
495.65
7432.02
0
4.13%
34207264
L8
77932.8
77134
472.58
5708.22
29953307
>10000
56492589
>10000
96608294
L9
21374.5
21195
144.58
259.16
3615309
21933.1
21573
271.82
488.21
0
2.61%
6235639
21954.4
21373
365.09
882.13
1
2.71%
11222642
L10
39113.6
38850
210.16
867.88
9801280
40401.4
39762
430.52
1631.88
0
3.29%
18244990
40463.4
39831
354.59
2982.57
0
3.45%
32706226
L11
54512.3
54081
373.65
1787.21
11680795
56645.6
54883
1122.00
4209.09
0
3.91%
20733929
57025.3
56219
617.77
8397.48
0
4.61%
35694937
L12
110389.7
109834
765.05
6045.78
34276626
>10000
71917382
>10000
121540613From the computational results above, we can conclude the following:(1)
TSSA is much faster than SA but has no obvious advantage in terms of the quality of the solution. This is because the initial temperature of TSSA is lower, and its total number of iterations is less.(2)
Although TSSALSS has no significant advantage over TSSA in the small-scale numerical examples, it generates higher-quality solutions within shorter time frames in both the medium-scale and large-scale numerical examples. Moreover, in Figure3, it can be seen that the solution obtained using TSSALSS tends to improve in quality as the scale increases.(3)
The convergence rate of these three algorithms decreases in the order of TSSALSS, TSSA, and SA, which is shown in Figure4. Figure 4 is the computational result of L9 using the three algorithms. Note that the abscissa axis is the cooling times rather than the total number of iterations. The epoch length of TSSALSS, TSSA, and SA at each temperature is 35000, 50000, and 50000, respectively (see Appendix A, row L9).Figure 3
Trends in the gaps.Figure 4
Computational result for L9 using the three algorithms.In addition, we compared our algorithm (TSSALSS) with three common metaheuristic algorithms from the literature: simulated annealing with greedy search (SA-GS) [38], genetic algorithm (GA) [30], and tabu search (TS) [29]. The neighborhood structures of SA-GS and TS are the same as those proposed in the Section 3.2. The crossover operators and mutation operators of GA are newly designed for the proposed problem. The mutation operators are the same as the neighborhood structures in Section 3.2. The crossover operators include the task allocation crossover and the routes crossover. The task allocation crossover will generate illegal solutions; thus we legalize these illegal solutions once they occur. The routes’ crossover is similar to that of the traveling salesman problem (TSP) [41]. These three algorithms were all executed ten times in the same running condition as that of TSSALSS. The results obtained are listed in Table 3, where the columns “average,” “minimum,” and “T” indicate the objective function value of the best solution among the ten trials, average objective function value of the ten trials, and average time spent by these algorithms.Table 3
Comparison of SA-GS, TS, GA, and TSSALSS.
Instances
TSSALSS
SA-GS
TS
GA
Average
Minimum
T
Average
Minimum
T
Average
Minimum
T
Average
Minimum
T
S1
304.901
297.15
1.73
302.733
296.80
9.85
306.544
300.18
5.92
369.470
359.03
23.43
S2
456.536
449.86
1.88
453.692
447.84
10.79
457.090
446.83
6.43
550.552
539.90
32.75
S3
602.408
585.28
5.67
600.790
586.33
39.02
599.563
588.24
19.95
728.918
712.55
109.65
S4
879.412
848.56
5.70
880.433
850.10
36.73
883.270
853.16
18.28
1066.26
1032.8
123.67
S5
501.301
493.02
3.59
500.094
494.37
21.69
500.032
492.16
11.60
615.338
606.99
92.34
S6
629.063
612.58
3.43
624.186
610.55
20.38
629.761
614.59
12.61
773.413
753.07
91.45
S7
934.877
902.45
4.05
930.638
901.27
30.76
936.710
909.08
16.67
1146.19
1113.0
124.53
S8
1239.98
1192.3
4.47
1236.74
1207.8
37.09
1248.65
1200.4
15.26
1518.47
1480.3
138.05
M1
1788.32
1748.3
9.44
1789.61
1747.5
89.83
1803.16
1763.0
49.87
2232.40
2212.0
384.95
M2
3056.33
2961.1
10.38
3045.58
2964.2
83.74
3082.92
2986.6
45.02
3811.97
3737.0
429.05
M3
3362.66
3316.6
23.01
3379.41
3317.3
203.34
3388.22
3348.8
82.99
4232.26
4152.9
948.38
M4
4304.45
4244.3
24.32
4304.95
4237.4
219.50
4347.06
4280.8
95.04
5446.28
5371.0
1010.25
M5
6668.53
6596.3
73.86
6670.79
6600.3
718.94
6675.87
6662.9
282.29
8420.00
8349.1
3415.91
M6
10201.1
10107
76.35
10196.3
10116
711.86
10384.7
10226
439.45
12848.1
12751
3794.89
M7
5933.78
5880.7
74.02
5947.83
5890.1
704.68
5986.00
5940.1
297.43
7576.72
7479.4
3596.00
M8
9039.82
8905.8
79.19
9038.74
8905.3
731.35
9198.02
9065.2
359.17
11483.2
11435
4117.48
M9
10582.9
10396
160.04
10593.4
10403
1319.08
10826.3
10651
735.17
14028.2
13371
9126.51
M10
13712.1
13523
184.77
13717.8
13531
1503.93
14096.0
13926
888.05
>10000
L1
9829.93
9766.5
83.20
9837.62
9774.0
788.55
10019.6
9951.1
378.80
12693.8
12626
4488.83
L2
18127.8
17930
203.29
18139.5
17934
1737.07
18780.4
18558
897.94
>10000
L3
25307.3
25023
358.27
25316.2
25091
3198.06
26261.4
26026
1829.06
>10000
L4
44704.2
44276
1532.21
>10000
46890.2
46547
8276.54
>10000
L5
20877.3
20651
224.84
20868.0
20635
1968.49
21526.6
21310
1159.61
>10000
L6
39325.9
38912
824.93
39335.6
38934
7034.29
40950.1
40387
3759.47
>10000
L7
38951.6
38327
1675.24
>10000
40809.6
40309
8570.97
>10000
L8
77932.8
77134
5708.22
>10000
>10000
>10000
L9
21374.5
21195
259.16
21387.4
21293
2217.22
22163.2
21998
1231.80
>10000
L10
39113.6
38850
867.88
39141.1
38904
7419.42
40431.7
40435
4159.30
>10000
L11
54512.3
54081
1787.21
>10000
>10000
>10000
L12
110389.7
109834
6045.78
>10000
>10000
>10000As the table shows, TSSALSS spends the least amount of time. Although the quality of the solutions obtained by TSSALSS is not particularly high compared with that of the other two algorithms (SA-GS and TS), especially SA-GS, when the scale of the problems is relatively small, the advantage of TSSALSS becomes evident as the scale of the instances increases. GA is less efficient and less precise compared with the other three algorithms because the crossover operators generate many illegal solutions and it spends much time to tackle these illegal solutions.


## 5. Case Study

### 5.1. Case One

The first case study in this paper focuses on the inbound logistics routing problem for a feed factory in Shandong, China. This feed factory has five plants and twelve suppliers nearby. The location of every plant and supplier, tasks, and pipeline inventory costs per unit time are given in AppendixB. We optimized the following three transportation modes.(1)
Mode 1: Milk-run among all plants and suppliers, that is,n-n.(2)
Mode 2: Milk-run among all suppliers for each plant, that is,n-1.(3)
Mode 3: Milk-run among all plants for each supplier, that is,1-n.The routes and task allocation plans of the three models are given in Tables 4–7.Table 4
Vehicle routes of Mode 1.
Vehicles
Routes of each vehicle
1
O-G-B-F-N-Q-P-O
2
Q-J-Q
3
M-D-C-P-M
4
O-J-E-M-O
5
Q-H-Q
6
N-K-A-F-O-M-P-Q-N
7
P-H-E-P
8
O-H-O
9
M-J-M
10
N-A-F-Q-N
11
N-H-N
12
M-I-L-P-Q-N-O-M
13
M-H-M
14
M-G-A-M
15
P-J-P
16
Q-J-L-E-N-Q
17
P-D-C-N-Q-P
18
M-B-F-O-MTable 5
Task allocation plan of Mode 1.
Plants
Suppliers
A
B
C
D
E
F
G
H
I
J
K
L
M
14∗
18
3
0
4
18
14
13
12
9
6
12
N
10
0
0
17
16
1
1
11
12
16
6
16
O
6
18
3
3
0
6
1
8
12
4
6
0
P
6
1
17
17
7
6
1
7
12
15
6
0
Q
10
1
17
17
16
10
1
5
12
2
6
12
∗14 indicates that task DMA is performed by vehicle 14. If DMA is zero, ∗ will be zero.Table 6
Vehicle routes of Mode 2.
Plants
Routes for each plant
M
0-B-F-0-G-A-0-J-0-E-L-I-0-C-K-0-H-0
N
0-G-I-L-J-E-D-0-H-F-A-K-0
O
0-F-B-A-K-C-D-0-J-I-G-0-H-0
P
0-F-A-B-G-I-0-E-H-0-J-0-K-C-D-0
Q
0-F-A-B-G-0-E-J-0-H-I-L-0-D-C-K-0Table 7
Vehicle routes of Mode 3.
Suppliers
Routes for each supplier
A
0-M-P-0-Q-N-0-O-0
B
0-Q-P-M-0-O-0
C
0-Q-P-0-M-O-0
D
0-N-O-Q-P-0
E
0-P-N-Q-M-0
F
0-O-M-P-Q-N-0
G
0-M-P-Q-N-O-0
H
0-M-0-N-0-O-0-P-0-Q-0
I
0-M-O-N-Q-P-0
J
0-M-0-O-0-P-0-N-0-Q-0
K
0-N-Q-P-M-O-0
L
0-Q-N-M-0Table8 presents the results for the three transportation plans. For the total objective function, Mode 1 is 24.8% less than Mode 2 and 22.7% less than Mode 3. For traveling distance, Mode 1 is 25.3% less than Mode 2 and 24.4% less than Mode 3. Although the pipeline inventory cost of Mode 3 is only 18.3% lower than that of Mode 1, the number of vehicles employed in Mode 3 is 33.3% more than that of Mode 1. From the discussions above, it is obvious that Mode 1 holds significant advantages for logistics systems comprising multiple suppliers and plants.Table 8
Results for the three transportation plans.
Transportation modes
Objective function
Total distance
Pipeline inventory cost
Number of vehicles
Mode 1
894.6
841.5
53.1
18
Mode 2
1190.1
1126.5
63.6
19
Mode 3
1156.9
1113.5
43.4
24


### 5.2. Case Two

In order to demonstrate the validity of the many-to-many Milk-run mode in the case of large scale, we use the computer to generate ann-n problem randomly, in which the case has ten plants and thirty suppliers. The location of every plant and supplier, volumes of tasks, and pipeline inventory costs per unit time are given in Appendix C. We similarly optimized the three transportation modes mentioned in Section 5.1. The routes and task allocation plans of the three modes are given in Tables 9–12. Suppliers are numbered from 1 to 30, and plants are numbered from A to J.Table 9
Vehicle routes of Mode 1.
Vehicles
Routes of each vehicle
1
C-29-14-9-28-24-30-H-J-C
2
A-19-18-25-8-I-E-A
3
F-17-23-16-15-26-G-F
4
H-21-6-11-2-B-J-H
5
B-28-9-14-29-13-3-4-H-B
6
H-11-1-23-A-E-I-H
7
I-29-14-15-E-A-I
8
F-28-9-29-14-G-F
9
J-8-25-18-19-7-H-J
10
B-2-1-11-30-9-G-F-B
11
E-28-9-4-I-A-E
12
D-19-5-7-30-4-C-D
13
I-3-13-26-22-A-E-I
14
D-3-9-14-13-C-D
15
D-27-26-22-C-D
16
C-29-12-10-D-C
17
A-6-21-20-I-E-A
18
D-8-24-18-25-28-C-D
19
E-10-12-27-26-I-A-E
20
I-26-27-13-3-4-H-J-I
21
B-11-1-H-B
22
G-22-20-21-6-8-F-G
23
D-15-17-23-16-I-C-D
24
J-26-27-12-10-20-H-B-J
25
C-7-5-2-11-I-A-E-D-C
26
F-19-5-7-8-G-F
27
J-1-6-21-20-C-D-J
28
J-17-23-16-15-22-H-B-J
29
B-7-30-24-8-25-18-19-5-B
30
F-25-18-24-4-3-G-F
31
H-7-5-30-20-22-I-J-H
32
G-12-10-27-13-F-G
33
J-17-23-16-E-A-J
34
A-30-24-I-E-ATable 10
Task allocation plan of Mode 1.
Plants
Suppliers
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
A
6
25
13
11
25
17
25
2
11
19
6
19
13
7
7
33
33
2
2
17
17
13
33
34
2
13
19
11
7
34
B
21
10
5
5
29
4
29
29
5
24
21
24
5
5
28
28
28
29
29
24
4
28
28
29
29
24
24
5
5
29
C
27
25
14
12
12
27
12
18
14
16
25
16
14
14
23
23
23
18
12
27
27
15
23
18
18
15
15
1
16
12
D
27
25
14
12
12
27
12
18
14
16
25
16
14
14
23
23
23
18
12
27
27
15
23
18
18
15
15
18
16
12
E
6
25
13
11
25
17
25
2
11
19
6
19
13
7
7
33
33
2
2
17
17
13
33
34
2
13
19
11
7
34
F
10
10
30
30
26
22
26
22
8
32
10
32
32
8
3
3
3
30
26
22
22
22
3
30
30
3
32
8
8
10
G
10
10
30
30
26
22
26
26
10
32
10
32
32
8
3
3
3
30
26
22
22
22
3
30
30
3
32
8
8
10
H
21
4
20
20
31
4
9
9
1
24
21
24
20
5
28
28
28
9
9
31
4
31
6
1
9
20
24
1
5
1
I
6
25
13
11
31
17
31
2
11
19
6
19
13
7
7
23
23
2
2
17
17
31
6
34
2
19
20
11
7
34
J
27
4
20
20
31
4
9
9
1
24
4
24
20
1
28
33
33
9
9
31
4
31
33
1
9
20
20
1
1
31Table 11
Vehicle routes of Mode 2.
Plants
Routes for each plant
A
A-20-24-19-18-25-8-28-9-4-22-A-23-17-1-11-2-5-7-30-6-21-A-26-3-13-14-29-12-10-27-15-16-A
B
B-1-17-23-16-15-22-20-21-6-30-B-2-11-B-26-27-10-12-29-14-9-25-18-B-7-24-4-3-13-28-8-19-5-B
C
C-13-14-29-12-10-27-C-4-30-7-5-2-11-1-6-21-20-C-26-15-17-23-16-22-C-24-19-18-25-8-28-9-3-C
D
D-15-16-23-17-6-21-20-22-26-D-10-D-13-3-25-18-19-5-7-2-1-11-30-D-27-4-24-8-28-9-14-29-12-D
E
E-22-3-9-28-8-25-18-19-24-4-20-E-23-E-15-27-10-12-29-14-13-26-E-16-21-6-30-7-5-2-11-1-17-E
F
F-28-9-13-27-10-12-29-14-F-25-18-19-8-F-24-30-7-5-2-11-1-6-21-20-F-3-26-15-17-23-16-22-4-F
G
G-29-14-G-25-18-5-7-19-8-28-9-G-24-30-2-11-1-6-21-20-4-3-G-12-10-27-15-17-23-16-22-26-13-G
H
H-1-11-2-7-5-19-18-25-8-24-30-H-27-10-12-29-14-13-3-9-28-4-H-17-23-16-15-26-22-20-21-6-H
I
I-16-23-17-1-11-2-5-7-30-6-I-26-27-10-12-29-14-13-3-9-28-I-15-I-22-4-8-25-18-19-24-21-20-I
J
J-22-3-13-29-14-9-28-4-21-6-J-20-26-27-12-10-15-16-23-17-J-11-30-24-8-25-18-19-5-7-2-1-JTable 12
Vehicle routes of Mode 3.
Suppliers
Routes for each supplier
1
1-J-A-E-I-D-C-G-F-H-B-1
2
2-H-J-A-E-I-D-C-G-F-2-B-2
3
3-G-F-B-H-J-A-E-I-D-C-3
4
4-H-B-J-A-E-I-D-C-G-F-4
5
5-B-5-H-J-A-E-I-D-C-G-F-5
6
6-J-B-H-6-F-G-C-D-E-A-I-6
7
7-F-G-C-D-I-E-A-J-7-H-B-7
8
8-F-G-8-B-H-J-A-E-I-D-C-8
9
9-C-D-I-E-A-J-H-B-F-9-G-9
10
10-D-10-C-I-E-A-J-H-B-F-G-10
11
11-H-J-A-E-I-D-C-G-F-B-11
12
12-D-C-I-E-A-J-H-B-F-G-12
13
13-C-D-I-E-A-J-H-B-F-G-13
14
14-G-C-14-D-I-E-A-J-H-B-F-14
15
15-E-A-I-J-H-B-F-G-C-D-15
16
16-E-A-J-H-B-F-G-C-D-I-16
17
17-A-E-I-D-C-G-F-B-H-J-17
18
18-F-G-C-D-I-E-A-J-H-B-18
19
19-F-G-C-D-I-E-A-J-H-B-19
20
20-I-A-E-D-C-G-F-B-J-H-20
21
21-H-B-J-21-F-G-C-D-E-A-I-21
22
22-D-C-G-F-B-H-J-A-E-I-22
23
23-E-I-D-C-G-F-B-H-J-A-23
24
24-I-A-E-D-C-G-F-24-B-J-H-24
25
25-B-H-J-A-E-I-D-C-G-F-25
26
26-I-E-A-J-H-B-F-G-C-26-D-26
27
27-C-27-D-G-F-B-H-J-A-E-I-27
28
28-H-B-J-A-E-I-D-C-G-28-F-28
29
29-D-29-C-I-E-A-J-H-B-F-G-29
30
30-H-B-J-A-E-I-D-C-G-F-30Table13 presents the results for the three transportation plans. From the comparison between Mode 1 and Mode 2, although Mode 1 has a greater pipeline inventory cost (853.6) than Mode 2 does (803.9), the total cost and the total distance of Mode 1 (3036.2 and 2182.6, resp.) are both smaller than those of Mode 2 (3285.6 and 2481.7, resp.). Furthermore, Mode 1 needs three vehicles less than Mode 2 does. From the comparison between Mode 1 and Mode 3, all the four values of Mode 1 (objective function, total distance, pipeline inventory cost, and number of vehicle) are lower than those of Mode 3. From the analysis of Case Two, we can demonstrate the superiority of many-to-many transportation mode in the case of large-scale instances.Table 13
Results for the three transportation plans.
Transportation modes
Objective function
Total distance
Pipeline inventory cost
Number of vehicles
Mode 1
3036.2
2182.6
853.6
34
Mode 2
3285.6
2481.7
803.9
37
Mode 3
4082.7
5056.4
973.7
44


## 5.1. Case One

The first case study in this paper focuses on the inbound logistics routing problem for a feed factory in Shandong, China. This feed factory has five plants and twelve suppliers nearby. The location of every plant and supplier, tasks, and pipeline inventory costs per unit time are given in AppendixB. We optimized the following three transportation modes.(1)
Mode 1: Milk-run among all plants and suppliers, that is,n-n.(2)
Mode 2: Milk-run among all suppliers for each plant, that is,n-1.(3)
Mode 3: Milk-run among all plants for each supplier, that is,1-n.The routes and task allocation plans of the three models are given in Tables 4–7.Table 4
Vehicle routes of Mode 1.
Vehicles
Routes of each vehicle
1
O-G-B-F-N-Q-P-O
2
Q-J-Q
3
M-D-C-P-M
4
O-J-E-M-O
5
Q-H-Q
6
N-K-A-F-O-M-P-Q-N
7
P-H-E-P
8
O-H-O
9
M-J-M
10
N-A-F-Q-N
11
N-H-N
12
M-I-L-P-Q-N-O-M
13
M-H-M
14
M-G-A-M
15
P-J-P
16
Q-J-L-E-N-Q
17
P-D-C-N-Q-P
18
M-B-F-O-MTable 5
Task allocation plan of Mode 1.
Plants
Suppliers
A
B
C
D
E
F
G
H
I
J
K
L
M
14∗
18
3
0
4
18
14
13
12
9
6
12
N
10
0
0
17
16
1
1
11
12
16
6
16
O
6
18
3
3
0
6
1
8
12
4
6
0
P
6
1
17
17
7
6
1
7
12
15
6
0
Q
10
1
17
17
16
10
1
5
12
2
6
12
∗14 indicates that task DMA is performed by vehicle 14. If DMA is zero, ∗ will be zero.Table 6
Vehicle routes of Mode 2.
Plants
Routes for each plant
M
0-B-F-0-G-A-0-J-0-E-L-I-0-C-K-0-H-0
N
0-G-I-L-J-E-D-0-H-F-A-K-0
O
0-F-B-A-K-C-D-0-J-I-G-0-H-0
P
0-F-A-B-G-I-0-E-H-0-J-0-K-C-D-0
Q
0-F-A-B-G-0-E-J-0-H-I-L-0-D-C-K-0Table 7
Vehicle routes of Mode 3.
Suppliers
Routes for each supplier
A
0-M-P-0-Q-N-0-O-0
B
0-Q-P-M-0-O-0
C
0-Q-P-0-M-O-0
D
0-N-O-Q-P-0
E
0-P-N-Q-M-0
F
0-O-M-P-Q-N-0
G
0-M-P-Q-N-O-0
H
0-M-0-N-0-O-0-P-0-Q-0
I
0-M-O-N-Q-P-0
J
0-M-0-O-0-P-0-N-0-Q-0
K
0-N-Q-P-M-O-0
L
0-Q-N-M-0Table8 presents the results for the three transportation plans. For the total objective function, Mode 1 is 24.8% less than Mode 2 and 22.7% less than Mode 3. For traveling distance, Mode 1 is 25.3% less than Mode 2 and 24.4% less than Mode 3. Although the pipeline inventory cost of Mode 3 is only 18.3% lower than that of Mode 1, the number of vehicles employed in Mode 3 is 33.3% more than that of Mode 1. From the discussions above, it is obvious that Mode 1 holds significant advantages for logistics systems comprising multiple suppliers and plants.Table 8
Results for the three transportation plans.
Transportation modes
Objective function
Total distance
Pipeline inventory cost
Number of vehicles
Mode 1
894.6
841.5
53.1
18
Mode 2
1190.1
1126.5
63.6
19
Mode 3
1156.9
1113.5
43.4
24


## 5.2. Case Two

In order to demonstrate the validity of the many-to-many Milk-run mode in the case of large scale, we use the computer to generate ann-n problem randomly, in which the case has ten plants and thirty suppliers. The location of every plant and supplier, volumes of tasks, and pipeline inventory costs per unit time are given in Appendix C. We similarly optimized the three transportation modes mentioned in Section 5.1. The routes and task allocation plans of the three modes are given in Tables 9–12. Suppliers are numbered from 1 to 30, and plants are numbered from A to J.Table 9
Vehicle routes of Mode 1.
Vehicles
Routes of each vehicle
1
C-29-14-9-28-24-30-H-J-C
2
A-19-18-25-8-I-E-A
3
F-17-23-16-15-26-G-F
4
H-21-6-11-2-B-J-H
5
B-28-9-14-29-13-3-4-H-B
6
H-11-1-23-A-E-I-H
7
I-29-14-15-E-A-I
8
F-28-9-29-14-G-F
9
J-8-25-18-19-7-H-J
10
B-2-1-11-30-9-G-F-B
11
E-28-9-4-I-A-E
12
D-19-5-7-30-4-C-D
13
I-3-13-26-22-A-E-I
14
D-3-9-14-13-C-D
15
D-27-26-22-C-D
16
C-29-12-10-D-C
17
A-6-21-20-I-E-A
18
D-8-24-18-25-28-C-D
19
E-10-12-27-26-I-A-E
20
I-26-27-13-3-4-H-J-I
21
B-11-1-H-B
22
G-22-20-21-6-8-F-G
23
D-15-17-23-16-I-C-D
24
J-26-27-12-10-20-H-B-J
25
C-7-5-2-11-I-A-E-D-C
26
F-19-5-7-8-G-F
27
J-1-6-21-20-C-D-J
28
J-17-23-16-15-22-H-B-J
29
B-7-30-24-8-25-18-19-5-B
30
F-25-18-24-4-3-G-F
31
H-7-5-30-20-22-I-J-H
32
G-12-10-27-13-F-G
33
J-17-23-16-E-A-J
34
A-30-24-I-E-ATable 10
Task allocation plan of Mode 1.
Plants
Suppliers
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
A
6
25
13
11
25
17
25
2
11
19
6
19
13
7
7
33
33
2
2
17
17
13
33
34
2
13
19
11
7
34
B
21
10
5
5
29
4
29
29
5
24
21
24
5
5
28
28
28
29
29
24
4
28
28
29
29
24
24
5
5
29
C
27
25
14
12
12
27
12
18
14
16
25
16
14
14
23
23
23
18
12
27
27
15
23
18
18
15
15
1
16
12
D
27
25
14
12
12
27
12
18
14
16
25
16
14
14
23
23
23
18
12
27
27
15
23
18
18
15
15
18
16
12
E
6
25
13
11
25
17
25
2
11
19
6
19
13
7
7
33
33
2
2
17
17
13
33
34
2
13
19
11
7
34
F
10
10
30
30
26
22
26
22
8
32
10
32
32
8
3
3
3
30
26
22
22
22
3
30
30
3
32
8
8
10
G
10
10
30
30
26
22
26
26
10
32
10
32
32
8
3
3
3
30
26
22
22
22
3
30
30
3
32
8
8
10
H
21
4
20
20
31
4
9
9
1
24
21
24
20
5
28
28
28
9
9
31
4
31
6
1
9
20
24
1
5
1
I
6
25
13
11
31
17
31
2
11
19
6
19
13
7
7
23
23
2
2
17
17
31
6
34
2
19
20
11
7
34
J
27
4
20
20
31
4
9
9
1
24
4
24
20
1
28
33
33
9
9
31
4
31
33
1
9
20
20
1
1
31Table 11
Vehicle routes of Mode 2.
Plants
Routes for each plant
A
A-20-24-19-18-25-8-28-9-4-22-A-23-17-1-11-2-5-7-30-6-21-A-26-3-13-14-29-12-10-27-15-16-A
B
B-1-17-23-16-15-22-20-21-6-30-B-2-11-B-26-27-10-12-29-14-9-25-18-B-7-24-4-3-13-28-8-19-5-B
C
C-13-14-29-12-10-27-C-4-30-7-5-2-11-1-6-21-20-C-26-15-17-23-16-22-C-24-19-18-25-8-28-9-3-C
D
D-15-16-23-17-6-21-20-22-26-D-10-D-13-3-25-18-19-5-7-2-1-11-30-D-27-4-24-8-28-9-14-29-12-D
E
E-22-3-9-28-8-25-18-19-24-4-20-E-23-E-15-27-10-12-29-14-13-26-E-16-21-6-30-7-5-2-11-1-17-E
F
F-28-9-13-27-10-12-29-14-F-25-18-19-8-F-24-30-7-5-2-11-1-6-21-20-F-3-26-15-17-23-16-22-4-F
G
G-29-14-G-25-18-5-7-19-8-28-9-G-24-30-2-11-1-6-21-20-4-3-G-12-10-27-15-17-23-16-22-26-13-G
H
H-1-11-2-7-5-19-18-25-8-24-30-H-27-10-12-29-14-13-3-9-28-4-H-17-23-16-15-26-22-20-21-6-H
I
I-16-23-17-1-11-2-5-7-30-6-I-26-27-10-12-29-14-13-3-9-28-I-15-I-22-4-8-25-18-19-24-21-20-I
J
J-22-3-13-29-14-9-28-4-21-6-J-20-26-27-12-10-15-16-23-17-J-11-30-24-8-25-18-19-5-7-2-1-JTable 12
Vehicle routes of Mode 3.
Suppliers
Routes for each supplier
1
1-J-A-E-I-D-C-G-F-H-B-1
2
2-H-J-A-E-I-D-C-G-F-2-B-2
3
3-G-F-B-H-J-A-E-I-D-C-3
4
4-H-B-J-A-E-I-D-C-G-F-4
5
5-B-5-H-J-A-E-I-D-C-G-F-5
6
6-J-B-H-6-F-G-C-D-E-A-I-6
7
7-F-G-C-D-I-E-A-J-7-H-B-7
8
8-F-G-8-B-H-J-A-E-I-D-C-8
9
9-C-D-I-E-A-J-H-B-F-9-G-9
10
10-D-10-C-I-E-A-J-H-B-F-G-10
11
11-H-J-A-E-I-D-C-G-F-B-11
12
12-D-C-I-E-A-J-H-B-F-G-12
13
13-C-D-I-E-A-J-H-B-F-G-13
14
14-G-C-14-D-I-E-A-J-H-B-F-14
15
15-E-A-I-J-H-B-F-G-C-D-15
16
16-E-A-J-H-B-F-G-C-D-I-16
17
17-A-E-I-D-C-G-F-B-H-J-17
18
18-F-G-C-D-I-E-A-J-H-B-18
19
19-F-G-C-D-I-E-A-J-H-B-19
20
20-I-A-E-D-C-G-F-B-J-H-20
21
21-H-B-J-21-F-G-C-D-E-A-I-21
22
22-D-C-G-F-B-H-J-A-E-I-22
23
23-E-I-D-C-G-F-B-H-J-A-23
24
24-I-A-E-D-C-G-F-24-B-J-H-24
25
25-B-H-J-A-E-I-D-C-G-F-25
26
26-I-E-A-J-H-B-F-G-C-26-D-26
27
27-C-27-D-G-F-B-H-J-A-E-I-27
28
28-H-B-J-A-E-I-D-C-G-28-F-28
29
29-D-29-C-I-E-A-J-H-B-F-G-29
30
30-H-B-J-A-E-I-D-C-G-F-30Table13 presents the results for the three transportation plans. From the comparison between Mode 1 and Mode 2, although Mode 1 has a greater pipeline inventory cost (853.6) than Mode 2 does (803.9), the total cost and the total distance of Mode 1 (3036.2 and 2182.6, resp.) are both smaller than those of Mode 2 (3285.6 and 2481.7, resp.). Furthermore, Mode 1 needs three vehicles less than Mode 2 does. From the comparison between Mode 1 and Mode 3, all the four values of Mode 1 (objective function, total distance, pipeline inventory cost, and number of vehicle) are lower than those of Mode 3. From the analysis of Case Two, we can demonstrate the superiority of many-to-many transportation mode in the case of large-scale instances.Table 13
Results for the three transportation plans.
Transportation modes
Objective function
Total distance
Pipeline inventory cost
Number of vehicles
Mode 1
3036.2
2182.6
853.6
34
Mode 2
3285.6
2481.7
803.9
37
Mode 3
4082.7
5056.4
973.7
44


## 6. Conclusions

In this paper, we considered the problem of Milk-run route planning with multiple suppliers and plants and proposed and constructed a new mathematical model that minimizes transportation and pipeline inventory costs. To solve the proposed model, we developed a two-stage heuristic algorithm that comprehensively considers the characteristics of the problem. The first stage of the algorithm minimizes the total number of service nodes visited by all vehicles, while the second stage limits the search scope of SA to improve search efficiency. The results obtained from thirty small-, medium-, and large-scale numerical examples verified the efficacy of the algorithms proposed in this paper, and their performance proved to be more effective than the traditional SA. Finally, we studied the two cases and confirmed that then-n Milk-run transportation mode is superior to other modes.The work presented herein has the following limitations: (1) because of the complexity of the problem, the proposed algorithm was not compared with exact algorithms; (2) in this paper, we only considered the case in which the transportation frequency is one, whereas, in practice, the frequency of transportation within a day may be performed more than once to reduce the inventory of each plant. As a consequence, we plan to extend the model to take this into account. We also plan to consider long-distance transportation and extend the model to include the cross-docking strategy. Further, we will consider the effect of using different types of vehicles to perform the transportation tasks.



---

*Source: 428925-2015-08-31.xml*