An Enhanced Approach for the Multiple Vehicle Routing Problem with Heterogeneous Vehicles and a Soft Time Window

: The vehicle routing problem (VRP) is a challenging combinatorial optimization problem. This research focuses on the problem under which a manufacturer needs to outsource materials from other suppliers and to ship the materials back to the company. Heterogeneous vehicles are available to ship the materials, and each vehicle has a limited loading capacity and a limited travelling distance. The purpose of this research is to study a multiple vehicle routing problem (MVRP) with soft time window and heterogeneous vehicles. Two models, using mixed integer programming (MIP) and genetic algorithm (GA), are developed to solve the problem. The MIP model is ﬁrst constructed to minimize the total transportation cost, which includes the assignment cost, travelling cost, and the tardiness cost, for the manufacturer. The optimal solution can present multiple vehicle routing and the loading size of each vehicle in each period. The GA is next applied to solve the problem so that a near-optimal solution can be obtained when the problem is too difﬁcult to be solved using the MIP. A case of a food manufacturing company is used to examine the practicality of the proposed MIP model and the GA model. The results show that the MIP model can obtain the optimal solution under a short computational time when the scale of the problem is small. When the problem becomes non-deterministic polynomial hard (NP-hard), the MIP model cannot ﬁnd the optimal solution. On the other hand, the GA model can obtain a near-optimal solution within a reasonable amount of computational time. This paper is related to several important topics of the Symmetry journal in the areas of mathematics and computer science theory and methods. In the area of mathematics, the theories of linear and non-linear algebraic structures and information technology are adopted. In the area of computer science, theory and methods, and metaheuristics are applied.


Introduction
In today's competitive business environment, firms need to optimize their supply chains, and thus supply chain management is very important for firms' sustainability. Logistics is one of the main areas in supply chain management. It is often related to production and inventory decisions, and the delivery cost usually accounts for a major part of the total logistic costs [1,2]. The vehicle routing problem (VRP) is a core issue in logistics, and it refers to a class of combinatorial optimization problems in which customers are to be served by a number of vehicles [3]. VRP was first introduced by Dantzig and Ramser [4], and the main objectives of the problem are to minimize the total travelling cost, time, or distance with a fleet of vehicles, starting and ending their routes at the depot while A multi-objective hybrid approach, which integrated multi-objective particle swarm optimization and adapted a multi-objective variable neighborhood search, was proposed. The results were found to outperform some genetic algorithm (GA)-based methods. Madankumar and Rajendran [29] studied the green VRP with pickup and delivery, and used the semiconductor supply chain as an example. Three MILP models were constructed. The first model was to minimize the total routing cost and to obtain the routes and schedules for the alternative fuel vehicles while satisfying a set of requests. The second model extended the first model and considered different fuel prices at different refueling stations. The third model considered the basic pickup and delivery problem without the consideration of the green supply chain, and it was compared with the first two models. Dalmeijer and Spliet [30] studied the time window assignment VRP, in which time windows for delivery are assigned to clients before the demand volume of clients is known. To deal with demand uncertainty, the problem assumes that the demand scenarios and their probabilities of occurrence are known. A flow formulation was presented, and a branch-and-cut algorithm was proposed, to obtain solutions quickly, and to solve large instances. Xu et al. [31] studied a dynamic VRP and proposed an enhanced ant colony optimization model, which incorporated improved K-means and crossover operation into the traditional ant colony optimization. The improved K-means aimed to divide the region with the most reasonable distance, and the crossover operation tried to extend the search space and tried to avoid falling into a local search prematurely.
GA has been applied in solving VRP. Some works are reviewed here. Baker and Ayechew [32] stated that unlike other heuristics, which has had considerable progress in solving VRP, GA had not made a great impact. The authors developed a conceptually straightforward GA to solve the basic VRP, and a hybrid heuristic, which incorporated a neighborhood search into the GA. A comparison with other heuristics showed that the GA was competitive in terms of computing time and solution quality. Ombuki et al. [11] studied the VRP with hard time windows, and treated it as a multi-objective problem with two objectives: the number of vehicles and the total distance traveled. A multi-objective GA approach using the Pareto ranking technique was developed so that the weighted sum, which might cause solution bias, did not need to be applied. Alvarenga et al. [33] studied the VRP with time windows and proposed a two-phase approach that incorporated GA and a set partitioning formulation. The objective was to minimize the travel distance, and the authors stated that the approach outperformed past heuristic methods in terms of the minimal travel distance. Vidal et al. [34] studied three variants of VRP, the multi-depot VRP, the periodic VRP, and the multi-depot periodic VRP, under an environment of capacitated vehicles and constrained route duration. Based on the GA, the authors proposed a hybrid genetic search with an adaptive diversity control metaheuristic. Vidal et al. [35] further proposed a hybrid genetic search with advanced diversity control to tackle four types of VRP: VRP with time windows, periodic VRP with time windows, multi-depot VRP with time windows, and site-dependent VRP with time windows. Some evaluation techniques were introduced to the GA to evaluate and prune neighborhoods, and to decompose large instances efficiently. Xiao and Konak [36] studied the green vehicle routing and the scheduling problem, in which a fleet of heterogeneous vehicles travel within a time-varying traffic environment to minimize the total CO 2 emissions. An MIP was first proposed for small-sized problems, and an exact dynamic programming algorithm was presented next for large-sized problems. A hybrid GA, combined with the exact dynamic programming algorithm, was then developed as an efficient solution approach.

Proposed Model
In this section, a model for the MVRP that considers a soft time window and heterogeneous vehicles is proposed by the MIP and the GA.

Assumptions
Some assumptions are made in the problem as follows: 1.
The decision model considers multiple periods, multiple suppliers, and multiple vehicles. The distance between each two suppliers is fixed and known. 4.
The travelling time between each two suppliers is fixed and known.

5.
Different types of vehicles have different unit travelling costs. 6.
The unloading time in each supplier is fixed and known. 7.
The assigning cost for each vehicle is fixed and known. Different assigning costs occur for different types of vehicles. 8.
Each vehicle has a limited loading capacity. Different types of vehicles have different loading capacities. 9.
Each vehicle has a limited travelling distance in a period. Different types of vehicles have different travelling distance limits. 10. Each supplier has a soft time window. When a vehicle arrives to a supplier after the latest soft time, a tardiness cost will be charged by the supplier based on the tardiness time. 11. Multiple vehicles can be used in a period. A supplier can only be visited by at most a vehicle in a period. 12. No shortage of outsourced materials is allowed.

Construction of the MIP Model
In this study, the total transportation cost includes the fixed cost for assigning vehicles, the travelling costs among the suppliers, and the tardiness cost when a vehicle arrives a supplier late. It is as follows: where τ v is the fixed cost for assigning vehicle v, X t,0,j,v is a binary variable, which indicates whether vehicle v departs from the depot (i = 0) to supplier j in period t, π i,j is the travelling distance from supplier i to supplier j, ρ v is the travelling cost per unit of distance, X t,i,j,v is a binary variable, which indicates whether vehicle v travels from supplier i to supplier j in period t, ρ i is the tardiness cost per unit of time, and L t,i,v is the tardiness time of vehicle v when arriving at supplier i in period t. An MIP model is presented to solve the transportation problem, which considers the vehicle assigning cost, travelling cost, and tardiness cost. The objective is to minimize the total transportation cost in the system. Minimize Subject to: Objective function (4) is to minimize the total transportation cost in the system. Constraint (5) ensures that at most, V vehicle(s) depart from the depot (node 0) in period t. Constraint (6) ensures that at most, V vehicle(s) return to the depot (node 0) in period t. Constraint (7) assures that exactly one vehicle leaves from supplier i in period t, and that i = 0 and i = j. Constraint (8) assures that exactly one vehicle enters to supplier j in period t, and that j = 0 and i = j. Constraint (9) states that if vehicle v visits supplier j in period t, it must also depart from supplier j in period t. Constraint (10) makes sure that the total loading size of vehicle v in period t must be less than or equal to the maximum loading size for that vehicle. Constraint (11) ensures that the total travelling distance of vehicle v in period t must be less than or equal to the maximum travelling distance for that vehicle. Constraint (12) assures that if vehicle v visits supplier j after supplier i in period t, the service start time for supplier j cannot begin earlier than the service start time for supplier i, plus the unload time at supplier i and the travelling time from supplier i to supplier j. Constraint (13) computes the tardiness time of vehicle v when visiting supplier i in period t. Constraint (14) is standard sub-tour elimination constraints. Constraint (15) states that X t,i,j,v is a binary variable, 1 indicating that vehicle v travels from supplier i to supplier j in period t, and 0 indicating that no travel is incurred.

Genetic Algorithm Model
In this research, the GA is adopted to solve the MVRP with a soft time window and heterogeneous vehicles, so that near-optimal solutions can be obtained in a short computational time for large-scale problems. The procedures are as follows [16,37,38]: Step 1. Code Scheme In the travelling salesman problem, there are basically five different vector schemes to represent a tour among cities: a path representation, an adjacency representation, an ordinal representation, a position listing representation, and an adjacency listing representation [39]. Due to its simplicity, the path representation has been frequently adopted in the GA to solve the travelling salesman problem.
The constraints used in the GA in this study are as follows: X t,i,j,v is a binary number, Q t,v , R t,v and S t,j,v are integer numbers. Constraint (16) ensures sure that the total loading size of vehicle v in period t must be less than or equal to the maximum loading size for that vehicle. Constraint (17) assures that the total travelling distance of vehicle v in period t must be less than or equal to the maximum travelling distance for that vehicle. Constraint (18) makes sure that if vehicle v visits supplier j after supplier i in period t, the service start time for supplier j cannot begin earlier than the service start time for supplier i plus the unload time at supplier i and the travelling time from supplier i to supplier j.

Step 2. Initial Population of Chromosomes
The initial population is generated randomly.

Step 3. Fitness Function
The fitness function is to minimize TC, where TC is the sum of the vehicle assignment cost, travelling cost, and tardiness cost. Min TC is the minimum cost among all the chromosomes across the population.
Step 4. Crossover Operation The standard two-cut-point crossover operator is adopted.

Step 5. Mutation Operator
A mutation operator is applied to change a randomly selected gene in the genetic code (0-1, 1-0).

Step 6. Selection of Subsequent Population
Individuals are sorted by their fitness values, and those with higher fitness values are more likely to be selected for the mating pool.

Step 7. Elitism Selection
In the GA applications, the concept of elitism, or the "best-so-far" technique, is widely adopted so that elitist carries the best chromosome from the parent generation to the subsequent generation [40].
Step 8. Termination Repeat the crossover, selection and replacement process until the objective function is optimized or the stop criterion is met.

Case Study
The case study is based on a manufacturer in the food industry. Assume that there are one manufacturer, several suppliers, and several vehicles. The objective is to minimize the total transportation cost, which includes the assignment cost, the travelling cost, and the tardiness cost. The software packages LINGO 10 (LINGO System Inc., Chicago, IL, USA) [41] and MATLAB 2015 (The MathWorks, Inc., Natick, MA, USA) [42] are used to solve the problem by the MIP and the GA, respectively.

Case Information
In this research, three cases are considered. There are multiple vehicles that need to be assigned to multiple suppliers in multiple periods. Three types of vehicles are available, large (L), medium (M), and small (S). The travelling distance limit, loading capacity, and assigning cost of each type of vehicles are different. The planning horizon covers multiple periods. The relevant information of the three cases is as follows. Table 1 shows the fixed cost for assigning a vehicle in a period, and the travelling cost per distance unit. Table 2 shows the loading and travelling distance limits of each vehicle. Table 3 shows the data of various suppliers: u i , p i , and l i . u i is the unload time (in minutes) required at supplier i, p i is the tardiness cost per minute charged by supplier i when a vehicle arrives after the latest soft time, and l i is the latest soft time (in minutes) to start the service at supplier i. Table 4 shows the travelling distance from one supplier (or depot) to the other. Table 5 shows the travelling time required from one supplier (or depot) to the other.    0  150  100  120  210  100  135  150  127  85  175  120  115  1  150  0  160  310  120  195  140  158  215  227  40  175  255  2  100  160  0  150  145  65  55  250  245  190  147  40  173  3  120  310  150  0  285  95  200  256  178  82  265  124  35  4  210  120  145  285  0  205  70  270  310  290  32  169  309  5  100  195  65  95  205  0  115  270  235  160  207  61  125  6  135  140  55  200  70  115  0  260  275  225  100  80  225  7  150  158  250  256  270  270  260  0  111  195  215  283  254  8  127  215  245  178  310  235  275  111  0  80  305  265  143  9  85  227  190  82  290  160  225  195  80  0  290  198  55  10  175  40  147  265  32  207  100  215  305  290  0  180  289  11  120  175  40  124  169  61  80  283  265  198  180  0  168  12  115  255  173  35  309  125  225  254  143  55  289  168  0 Unit of measure: km.  Table 6 lists the problem configurations for the three cases, including the number of periods, number of suppliers, number of vehicles, size of the transportation network, number of variables for the MIP, and number of constraints for the MIP. In the first case, there are five periods, the number of suppliers is five, the number of vehicles is three, the transportation network is 6 × 6, and the numbers of variables and constraints for the MIP model are 931 and 1171, respectively. In the second case, there are seven periods, the number of suppliers is nine, the number of vehicles is four, the transportation network is 10 × 10, and the numbers of variables and constraints for the MIP model are 3903 and 3977, respectively. In the third case, there are nine periods, the number of suppliers is 12, the number of vehicles is five, the transportation network is 13 × 13, and the numbers of variables and constraints for the MIP model are 9758 and 9919, respectively.

Case Results and Analysis
Case I In case I, there are five suppliers. There are three vehicles, one small vehicle and two medium vehicles. Vehicle 1 is a small truck, while Vehicles 2 and 3 are medium trucks. The planning horizon contains five periods. Table 7 shows the demand from each supplier in each period. The MIP model is solved by LINGO 10, and the results are as follows. The transportation route and loading sizes of each vehicle in each period are shown in Table 8. Because of the travelling distance and the loading size limits of each vehicle, three vehicles are required in each period. For example, Vehicle 1 travels from the depot to Supplier 3 and ships 40 units of materials back to the depot in period 1. Vehicle 2 travels from the depot to Supplier 1 to pick up 21 units of materials, to Supplier 4 to pick up 13 units of materials, and then back to the depot. Vehicle 3 travels from the depot to Supplier 2 to pick up 22 units, to Supplier 5 to pick up 26 units, and then back to the depot. The vehicle routing in Period 1 is depicted in Figure 1.
The GA model is next solved by MATLAB (2015). A two-cut-point crossover is adopted, and an inversion mutation operator is applied, to avoid a solution being trapped in a local optimum and to approach the global optimum. The initial population size is set as 50. The crossover rate is set as 0.85, that is, around 85% pairs of individuals are involved in the reproduction. The mutation rate is 0.15, indicating that each gene of a newly created solution is mutated, with a probability of 0.15. The algorithm is set to terminate at the 1000th generation. In case I, the best generation occurs at the 94th generation, as shown in Figure 2. The chromosome is coded as a string of i integer digits, and the ith gene is a binary number. When vehicle v travels from supplier node i to supplier node j in period t, the value for X t,i,j,v equals one; otherwise, the value equals zero. In Period 1, the result from the GA shows X 1,0,3,1 = 1 and X 1,3,0,1 = 1 for Chromosome 1, X 1,0,1,2 = 1, X 1,1,4,2 = 1 and X 1,4,0,2 = 1 for Chromosome 2, and X 1,0,2,3 = 1, X 1,2,5,3 = 1 and X 1,5,0,3 = 1 for Chromosome 3. The routing can also be observed in Figure 3. For example, in Period 1, Vehicle 1 travels from the depot to supplier 3 (X 1,0,3,1 = 1), and then back to the depot (X 1,3,0,1 = 1). The vehicle routing in Period 1 is depicted in Figure 1.  The GA model is next solved by MATLAB (2015). A two-cut-point crossover is adopted, and an inversion mutation operator is applied, to avoid a solution being trapped in a local optimum and to approach the global optimum. The initial population size is set as 50. The crossover rate is set as 0.85, that is, around 85% pairs of individuals are involved in the reproduction. The mutation rate is 0.15, indicating that each gene of a newly created solution is mutated, with a probability of 0.15. The algorithm is set to terminate at the 1000th generation. In case I, the best generation occurs at the 94th generation, as shown in Figure 2. The chromosome is coded as a string of i integer digits, and the th i gene is a binary number. When vehicle v travels from supplier node i to supplier node j in period t, the value for Xt,i,j,v equals one; otherwise, the value equals zero. In Period 1, the result from the GA shows X1,0,3,1 = 1 and X1,3,0,1 = 1 for Chromosome 1, X1,0,1,2 = 1, X1,1,4,2 = 1 and X1,4,0,2 = 1 for Chromosome 2, and X1,0,2,3 = 1, X1,2,5,3 = 1 and X1,5,0,3 = 1 for Chromosome 3. The routing can also be observed in Figure 3. For example, in Period 1, Vehicle 1 travels from the depot to supplier 3 (X1,0,3,1 = 1), and then back to the depot (X1,3,0,1 = 1). The results of Case I, generated by the MIP model and by the GA model are the same, as shown in Table 9. The assigning cost is $4000, the travelling cost is $4888, the tardiness cost is $143, and the total transportation cost is $9031. The computational times using the MIP and the GA are 0.5 s and 105 s, respectively. Thus, the MIP is effective in this case study. The results of Case I, generated by the MIP model and by the GA model are the same, as shown in Table 9. The assigning cost is $4000, the travelling cost is $4888, the tardiness cost is $143, and the total transportation cost is $9031. The computational times using the MIP and the GA are 0.5 s and 105 s, respectively. Thus, the MIP is effective in this case study.    The results of Case I, generated by the MIP model and by the GA model are the same, as shown in Table 9. The assigning cost is $4000, the travelling cost is $4888, the tardiness cost is $143, and the total transportation cost is $9031. The computational times using the MIP and the GA are 0.5 s and 105 s, respectively. Thus, the MIP is effective in this case study.

Case II
In Case II, there are nine suppliers and four vehicles, including one small vehicle, two medium vehicles and one large vehicle. Vehicle 1 is a small truck, Vehicles 2 and 3 are medium trucks, and Vehicle 4 is a large truck. The planning horizon contains seven periods. The demand information from nine suppliers is as shown in Table 10. The transportation route and loading sizes of each vehicle in each period under the MIP are shown in Table 11. Three vehicles are required in each period.      Figure 5 depicts the vehicle routing in Period 1 under the GA. Table 13 shows that the results from the MIP and the GA. Under the MIP, the assigning cost is $8400, the travelling cost is $11,521.5, the tardiness cost is $193.6, and the total transportation cost is $20,115.1. Under the GA, the assigning cost, travelling cost, tardiness cost, and the total transportation cost are $8400, $12,013.1, $505.4, and $20,918.5, respectively. The total transportation cost under the GA Symmetry 2018, 10, 650 13 of 20 is $803.4 higher than that under the MIP, with a difference of 3.99%. The computational time is 19 s under the MIP. Figure 6 shows the execution results from the GA. The size of the initial population is 100. The crossover rate and the mutation rate are 0.85 and 0.15, respectively. The algorithm is set to terminate at the 1000th generation. In the end, the computational time is 95 s, and the optimal solution is generated at the 133rd generation.  Table 13 shows that the results from the MIP and the GA. Under the MIP, the assigning cost is $8400, the travelling cost is $11,521.5, the tardiness cost is $193.6, and the total transportation cost is $20,115.1. Under the GA, the assigning cost, travelling cost, tardiness cost, and the total transportation cost are $8400, $12,013.1, $505.4, and $20,918.5, respectively. The total transportation cost under the GA is $803.4 higher than that under the MIP, with a difference of 3.99%. The computational time is 19 s under the MIP. Figure 6 shows the execution results from the GA. The size of the initial population is 100. The crossover rate and the mutation rate are 0.85 and 0.15, respectively. The algorithm is set to terminate at the 1000th generation. In the end, the computational time is 95 s, and the optimal solution is generated at the 133rd generation.

Case III
There are 15 factories in Case III. There are three kinds of vehicles, with one small vehicle, two medium vehicles, and two large vehicles. Vehicle 1 is a small truck, vehicle 2 and 3 are medium trucks, and vehicle 4 and 5 are large trucks. The planning horizon contains nine periods. The demand information from 12 suppliers is as shown in Table 14.

Case III
There are 15 factories in Case III. There are three kinds of vehicles, with one small vehicle, two medium vehicles, and two large vehicles. Vehicle 1 is a small truck, vehicle 2 and 3 are medium trucks, and vehicle 4 and 5 are large trucks. The planning horizon contains nine periods. The demand information from 12 suppliers is as shown in Table 14.  1  22  19  20  13  21  23  22  24  13  26  24  24  2  26  14  20  17  25  26  24  21  20  23  12  19  3  24  16  9  26  17  28  20  9  31  22  25  24  4  30  16  12  11  24  30  16  24  18  15  21  23  5  22  12  15  28  23  22  13  25  25  23  22  17  6  21  17  19  21  23  21  16  9  28  11  18  16  7  22  13  30  28  13  22  12  21  20  22  27  11  8  21  17  19  21  23  21  16  9  31  21  18  16  9  22  13  20  28  13  22  12  11  29  22  27  11 Because of the increase in the parameters and the decision variables, the problem becomes more complicated to be solved by the MIP. In fact, the problem in Case III is NP-hard. Therefore, only the GA is adopted to solve the problem. The results are shown in Tables 15 and 16. Table 15 shows the transportation route and loading sizes of each vehicle in each period under the GA, and all five vehicles are required in each period. Table 16 shows that the relevant costs results from the GA. The assigning cost is $14,400, the travelling cost is $17,200, the tardiness cost is $332.6, and the total transportation cost is $31,952.6. Figure 7 depicts the vehicle routing in Period 1 under the GA. Figure 8 shows the execution results from the GA. The size of the initial population is 150, the crossover rate is 0.85, and the mutation rate is 0.15. The algorithm is set to terminate at the 1000th generation. The computational time is 203 s, and the optimal solution is generated at the 312nd generation.

Conclusions
Vehicle routing problem (VRP) has been studied abundantly to determine the optimal set of routes for a fleet of vehicles to travel so that certain objectives can be achieved, such as reducing the total cost for a firm. To the authors' knowledge, this is the first instance of research that constructs a multiple vehicle routing problem (MVRP) model with a soft time window and heterogeneous vehicles. The contributions of this research are as follows. First, the MIP is used to construct a general formulation of the MVRP model. The model comprises the costs, including assigning cost, travelling cost, and tardiness cost for each planning period, and the goal aims to minimize the total transportation cost. Second, to tackle complex instances efficiently, a genetic algorithm (GA) model is proposed. Third, to examine the proposed MIP and GA models, two small-scale cases are studied, Figure 8. Results from the GA using MATLAB for Case III.

Conclusions
Vehicle routing problem (VRP) has been studied abundantly to determine the optimal set of routes for a fleet of vehicles to travel so that certain objectives can be achieved, such as reducing the total cost for a firm. To the authors' knowledge, this is the first instance of research that constructs a multiple vehicle routing problem (MVRP) model with a soft time window and heterogeneous vehicles. The contributions of this research are as follows. First, the MIP is used to construct a general formulation of the MVRP model. The model comprises the costs, including assigning cost, travelling cost, and tardiness cost for each planning period, and the goal aims to minimize the total transportation cost. Second, to tackle complex instances efficiently, a genetic algorithm (GA) model is proposed. Third, to examine the proposed MIP and GA models, two small-scale cases are studied, and both models can generate the optimal solutions. Fourth, a complicated case is used to examine the GA model. While the MIP model cannot solve the problem or it may require a very long time to perform the calculations if a case becomes too complicated, the proposed GA model can obtain near-optimal solutions in a short computational time. Therefore, the GA model can be effective in searching for near-optimal solutions. The proposed models, especially the GA model, can be valuable techniques for managers in real business settings. For example, based on the outcomes of the models, managers can determine the optimal or near optimal method for assigning the routings of multiple vehicles in each period, and for allocating loading sizes for each vehicle in each period, while aiming to minimize the total transportation cost.
In the future, a more comprehensive case for MVRP in manufacturing can be researched. Some assumptions are made in this research, such as fixed demand, fixed lead time, and one kind of product. A model that considers the probability demand, fuzzy lead time, fuzzy-demand joint replenishment, safety stock, and different priorities of orders can be constructed. In addition, the ordering cost for each purchase may be variable. Furthermore, a fuzzy multi-objective problem can be studied, and both a fuzzy multi-objective programming model and a GA model can be proposed to tackle the problem.
Other metaheuristics techniques, such as artificial immune system, particle swarm optimization, ant colony system and revised aforementioned metaheuristics, may also be adopted. The methods can be compared and contrasted to examine which method is more appropriate for solving the MVRP.