Research Article  Open Access
Introducing the MCHF/OVRP/SDMP: Multicapacitated/Heterogeneous Fleet/Open Vehicle Routing Problems with Split Deliveries and Multiproducts
Abstract
In this paper, we analyze a realworld OVRP problem for a production company. Considering realworld constrains, we classify our problem as multicapacitated/heterogeneous fleet/open vehicle routing problem with split deliveries and multiproduct (MCHF/OVRP/SDMP) which is a novel classification of an OVRP. We have developed a mixed integer programming (MIP) model for the problem and generated test problems in different size (10–90 customers) considering realworld parameters. Although MIP is able to find optimal solutions of small size (10 customers) problems, when the number of customers increases, the problem gets harder to solve, and thus MIP could not find optimal solutions for problems that contain more than 10 customers. Moreover, MIP fails to find any feasible solution of largescale problems (50–90 customers) within time limits (7200 seconds). Therefore, we have developed a genetic algorithm (GA) based solution approach for largescale problems. The experimental results show that the GA based approach reaches successful solutions with 9.66% gap in 392.8 s on average instead of 7200 s for the problems that contain 10–50 customers. For largescale problems (50–90 customers), GA reaches feasible solutions of problems within time limits. In conclusion, for the realworld applications, GA is preferable rather than MIP to reach feasible solutions in short time periods.
1. Introduction
Open vehicle routing problems (OVRPs) have gained much attention recently since they represent a problem type that needs to be solved by many production companies. In most industries, companies choose to use a hired vehicle fleet for distributing their goods. In this way, they do not have to endure the extra cost for returning vehicles since they use the resources of a thirdparty logistics (3PL) provider such as trucks or TIRs [1]. As consequences of such benefits, however, the companies have to accept some restrictions defined by the 3PL providers. For example, the 3PL provider can determine some particular routes considering the experiences of the drivers or the highway conditions. In addition, the 3PL provider can restrict the number of customers visited by a particular vehicle. Such specifications require some changes in the classical OVRP structure to make them more applicable for realworld problems.
In this paper, we consider a realworld vehicle routing problem for a production company. The company produces two different types of products (multiple products—MP) with different volume and weight properties. The first product is a lightweight but large product, such as styrofoam, and the second product has opposite volumeweight characteristics, heavyweight but small product, such as tar. Since the company uses a fleet of vehicles of a 3PL provider, the vehicles do not need to return to the depot, and thus, the underlying problem becomes an open vehicle routing problem (OVRP). The 3PL provider has a heterogeneous fleet of vehicles, such as trucks and TIRs (heterogeneous fleet—HF), with different volume and weight capacities (multiple capacitated—MC). In addition, on any route, the vehicles can serve two customers at most, and the demand of a customer can be supplied by different vehicles (split delivery—SD). Based on these problem specifications, we classify our problem as a multicapacitated/heterogeneous fleet/open vehicle routing problem with split deliveries and multiproducts (MCHF/OVRP/SDMP), which, to the best of our knowledge, has not yet been considered in previous studies.
Since the OVRP consists of Hamiltonian paths, we should find the best Hamiltonian path for each set of customers assigned to a vehicle for the optimal solution [2]. Therefore, it can be concluded that the OVRP has an NPhard structure because of the NPhard Hamiltonian path subproblems. Many researchers have developed various heuristic approaches to solve the OVRP. Brandão [2] developed a tabu search algorithm for the OVRP to generate good solutions in short time periods where the researcher obtained better results using a random tabu tenure instead of a fixed one [2]. Zachariadis and Kiranoudis [3] developed a novel approach to produce initial solutions for OVRPs. The researchers tested their metaheuristic on some wellknown OVRP instances where improvements on several bestknown solutions from previous studies were presented [3]. Eksioglu et al. [4] presented an extensive literature review for the vehicle routing problem (VRP), classified the VRP studies according to their specifications, and built the taxonomy of the VRP literature. According to the study, the VRP literature was categorized into five main groups that contain 106 different subcategories. They discovered that the number of the VRP studies including load specific and heterogeneous vehicles has been three times less than the others. We note that there have been a few VRP studies that contain heterogeneous vehicles in the literature. In this study, we consider an OVRP with heterogeneous fleet where we evaluate the loading process using the volume and weight coefficients of the products which has not yet been studied in detail so far.
Capacitated open vehicle routing problems (COVRPs) have also been studied by various researchers. Among these, Letchford et al. [5] formulated the first exact algorithm, and Simbolon [6] developed a direct search algorithm for the COVRP. In our study, we additionally consider the volume and weight capacities (multicapacitated—MC) of a heterogeneous fleet (HF) of vehicles. Just like the multicapacitated problems, the heterogeneous fleet of vehicles is also common in realworld transportation problems since a smaller vehicle with sufficient capacity would always be preferable to a larger one which, beyond decreasing the transportation cost, is also a basic principle of green logistics. Some studies considered different variations of these problems, such as Gendreau et al. [7] and Taillard [8] who developed a tabu search algorithm and a column generation method to solve the heterogeneous fleet vehicle routing problem (HVRP), respectively. The study of TavakkoliMoghaddam et al. [9] additionally considered the concept of split service in the capacitated vehicle routing problem (CVRP) where they showed that splitting demand implies a higher capacity utilization. In split delivery, customers might be visited by more than one vehicle similar to some reallife scenarios especially in the existence of a 3PL provider. Many researchers have studied the vehicle routing problems with split delivery properties. The study of Archetti and Speranza [10] presented a survey of split delivery vehicle routing problems (SDVRP) that includes the description of the SDVRP, its properties, exact algorithms and heuristics, and its variants and applications. None of the above studies, however, considered the open vehicle routing problem (OVRP) with a multiple capacitated heterogeneous fleet and multiple products with split deliveries together as in our study.
This paper is organized as follows. In Section 2, a mixedinteger programming (MIP) model is formulated. In Section 3 the proposed hybrid geneticlocal search is described. Computational results and conclusions are presented in Sections 4 and 5, respectively.
2. MixedInteger Programming Model
In this section, we present a mixedinteger programming model for the MCHF/OVRP/SDMP. A realworld problem of production company with its reallife assumptions was considered for the model. As mentioned before, the company uses a heterogeneous fleet of vehicles of a 3PL provider for distribution of goods, and the 3PL provider allows the vehicles to serve two customers at most. The 3PL provider determines the transportation cost considering the target customer. In addition, if a vehicle supplies the demand of an intermediate customer along with the target customer, the 3PL provider charges a stopping cost which depends on the vehicle type. We therefore calculate the total cost of the distribution process considering the transportation and the stopping cost of the vehicles together, which is aimed to be minimized in the MIP model. In our experiments, we consider two different types of products, however we define the type of materials with “” index, and we allow MIP to model the problems with multiple types of materials (more than two) simultaneously. The notation of the MIP model can be found in Table 1. Consider

The objective function (1) aims for the minimization of total cost which includes the transportation and the stopping cost of vehicles. To reduce the size of the solution space, we determine the maximum number of vehicles for the problem, and use this knowledge as a parameter for the index of vehicles. Constraint (2) ensures that each vehicle has only one target customer at most. Constraint (3) shows that each customer is decided as a target by only one vehicle at most. Constraint (4) guarantees that a vehicle can serve two customers at most; one, if applies, is the intermediate customer and the other is the target customer. The relationship between decision variables has been set by constraint (5). Constraint (6) ensures that the demands of customers are met by the vehicles. Constraints (7) and (8) represent the capacities of the vehicles in terms of weight and volume, respectively. Finally, (9) and (10) present the definition space of the decision variables.
Since the definition of the problem is unique, no benchmark for this type of problem is available. We therefore generate random samples for evaluating the performance of the MIP model which are detailed in Section 4.1. Although we can find the solutions with smaller gaps of smallsize problems (10–50 customers) using MIP, it is not able to find any feasible solution of largesize realworldlike problems (60–90 customers) within time limits (7200 s). To handle largesize problems, we develop a genetic algorithm with local search which is detailed in Section 3.
3. Hybrid GeneticLocal Search Algorithm
A genetic algorithm (GA) is a search technique based on the biological process of evolution theory that mimics natural selection. One of the important issues in GA is the genetic representation (string of symbols) of each solution in a population. The string is referred to as chromosome and the symbols as genes [11]. After generation of the initial population and determination of the fitness function values for each chromosome, GA manipulates the selection process by operations such as reproduction, crossover, and mutation [11]. The introduction of GA dates back to 1970s [12]. As the searching technique of genetic algorithms (GAs) [13] became popular in the mid1980s, many researchers started to apply the approach to different types of problems. To date, GA has been applied to many different types of problems including vehicle routing problems. In this paper, we also have developed a hybrid geneticlocal search algorithm for the MCHF/OVRP/SDMP problem which is detailed in the next section. The proposed local search algorithms, implemented in C#, are performed during the fitness function calculation process. In order to eliminate the infeasible chromosome structure, we were mainly inspired by our previous algorithms related to scheduling problem in the literature [14, 15]. The novel characteristic of the studied problem is the capacity constraint. The algorithm is not only applied to logistic area but also included capacities of vehicles and orders. The developed new approach related to capacity issue is integrated into the algorithm successfully and is explained in the next section.
3.1. Genetic Algorithm for Order Splitting Property
3.1.1. Encoding Scheme
In the proposed GA, a chromosome is mainly designed by a string of random numbers that are uniformly generated between and . These random keys show the transportation route for each vehicle. In addition, each chromosome also carries the number of suborders for each customer. It would be necessary to note that this chromosome structure is used in our previous study [14]. Figure 1 illustrates a sample chromosome. In the first section of the chromosome, the string contains (number of customers) segment. Each segment is further divided into (number of vehicle types) part, and then, each part for vehicle types is divided into genes. The number of genes for vehicles is two in this problem, which shows the maximum number of suborders. Note that the chromosome structure is convenient to increase the maximum suborder quantities. To decide the exact number of suborders for each customer, it is necessary to check the second section of the chromosome. In the second section of the chromosome, the string contains (number of customers) segment which shows the number of suborders for each customer. These numbers might be or for this problem, and they are generated randomly to reflect the variable suborder environment. For the chromosome structure in Figure 1, we have 2 vehicle types and 3 customers. For , since the number of suborders equals 2 as in the second section’s first segment of the chromosome, the smallest two random numbers will be selected from the first section of the chromosome among all the numbers generated for . The selected values (0.2 and 0.3 for ) are in bold in Figure 1. It means that the first order of will be transported by and the second order of will be transported by . Similarly, the number of suborders is 2 for and they will be transported by . The number of suborders is 1 for and it will be transported by .
3.1.2. Fitness Function
In the proposed algorithm, the chromosome structure defines the routes for each vehicle type. But we also need to integrate capacity conditions into the algorithm. Figure 2 shows the outline of integration process. We consider the volume and weight constraints of vehicles as well as the demands of each customer. After obtaining the routes for the customers—which will be explained in detail shortly after as in Figure 3—we need to iterate the integration process in two phases as shown in Figure 2. In the first phase, if the demands of customers are greater than a full load vehicle, then the number of fully loaded vehicles and costs are calculated. In the second phase of the algorithm, the remaining demands of customers are trying to be consolidated to maximize vehicles’ utilizations. After assigning the first demand to the vehicle, if the capacity of the existing vehicle is enough, route is the same, and the number of customers on the route is less than two, then the next demand can be loaded into the vehicle. This loading demand might be the suborder of the same product, the demand of the second product of the same customer, or the demand of the next customer on the route. Otherwise, the new vehicle is required. To explain how the algorithm obtains the routes for the customers and calculate total costs, the example in Figure 1 is considered. According to the chromosome shown in Figure 1, the number of suborders for , , and is determined as two, two, and one, respectively. will supply the 1st suborder of and the suborder of . will supply the 1st and the 2nd suborders of and the 2nd suborder of . The sequence of customers on will be determined by the random key numbers (0.20 and 0.43). Increasing arrangement of these random numbers will designate ’s customer sequence. Similarly, the sequence of customers on will be determined by the random key numbers (0.11, 0.26, and 0.30). Figure 3 shows the vehicle types and the customers for the examining problem. According to Figure 2, we need to calculate the full load vehicles for each customer’s suborders, considering two types of products. Then for the remaining loads, we need to combine the vehicles on the same route. For the considered example, three customers and two vehicle types, we will ignore the vehicles which are fully loaded in order to simplify the explanation.
As mentioned earlier, two types of costs are considered in the problem: stopping costs and transportation costs for each type of vehicle. For the considered example, each vehicle type’s route and customers’ loads (information about splitting or not splitting) are shown in Figure 3. Assume that Table 2 shows the stopping costs for each type of vehicle, and Tables 3 and 4 present the transportation costs of and , respectively. Transportation costs from the depot to the corresponding customers are shown in the first rows of Tables 3 and 4. The total cost is calculated by utilizing the transportation and stopping costs. Considering the first vehicle of , the transportation cost from depot to is 60, and the stopping cost of is 40, so the partial total cost is 60 + 40 = 100, and then the same vehicle transports from to and the transportation cost is 20. Similarly the stopping cost of is 40. It means that the total cost of is 100 + 20 + 40 = 160. The total cost of is greater than the total cost of . Our aim is to minimize the maximum total cost. As seen in Figure 3, the fitness value of chromosome in Figure 1 is 395. The structure of chromosome will be modified with the aim of finding the routes with minimum cost utilizing local search. The details of the GAbased local search algorithm are explained in Section 3.2. After the algorithm runs, we obtain the total cost for the chromosome. The population size determines the number of different chromosomes and fitness function values.



3.1.3. Genetic Operators
After generating the initial population, the operations of selection, crossover, and mutation are iteratively used to search for the best solution.(i)Selection: chromosomes are selected into the mating pool based on the random selection method [16]. In this method, parents are randomly chosen from the population.(ii)Crossover: the crossover operator, with crossover rate , is a method for sharing information between chromosomes. We use single point crossover approach in our algorithm which randomly chooses the crossing point and exchanges the genes between two parents to create the offspring. We use separate crossover operations in the first and second phases. Figure 4 illustrates the crossover operation.(iii)Mutation: the mutation operator is used to prevent the algorithm converging to a local optimum. In our algorithm, a mutation operation is performed as follows. Mutation operation is applied to the randomly selected chromosomes of the population according to the mutation rate (). The value of randomly selected gene of the selected chromosome is replaced with a new random number. By applying the operation to all selected chromosomes, we can obtain new chromosomes with new routes and fitness function values. The fitness value might be better or worse or remain same after applying the operator. An example for our problem is shown in Figure 5.
To enhance the performance of the genetic algorithm, we use a local search algorithm in the literature [14], which is detailed in the next subsection.
3.2. Local Search
Integrating one of the local search techniques within GA will commonly generate more competitive results. For instance, the integrating dominance properties method that was originally developed by Chang and Chen [17] was afterward applied successfully to machine scheduling problem [15]. In the other paper of researchers, which includes job splitting property in scheduling problem [14], an advance approach is adopted. And in this paper, finally the algorithm adopted more complex aforementioned logistics problem which contains capacity constraints. By searching all possible alternatives in a chromosome, we aim to reduce the total cost. After this process, the routes will be determined.
We consider all alternative situations for the local search as summarized in the Appendix. The proposed intervehicle type (interchange) and intravehicle type (exchange) customer changes are explained in Section 3.2.1 and the calibration of random numbers are explained in Section 3.2.2. The notation used in the pseudocode and some explanations on the local search operation are as follows.(i): the selected genes from the chromosome are ordered by vehicle type and then by customer sequence. The value of shows the total number of these genes. Each individual gene contains vehicle type and customer features. Figure 6 explains a sample gene structure of an example with seven customers and two vehicle types. The structure also indicates the sequence of customers for each type of vehicle.(ii) and : they are the gene numbers which are compared to decide whether the customers encoded in these genes should be exchanged or not.(iii): it indicates the type of vehicle encoded in the th gene. According to the structure shown in Figure 6, the type of vehicle for the 1st gene () is 1.(iv): it indicates the customer encoded in the th gene. According to the structure shown in Figure 6, the customer for the 1st gene () is Customer D.(v): it is the transportation cost between and for .(vi): it is the stopping cost of for .(vii): it is the total transportation cost for .
3.2.1. Interchange and Exchange of Customers
(i)Interchange of customers: there are two cases to be considered within the intervehicle type interchange: adjacent customers interchange and nonadjacent customers interchange.(I)Adjacent interchange: in the intervehicle type, adjacent customers interchange section of the pseudocode, steps are as follows.(1)If precedes , then “a” represents the total difference of transportation costs between after interchange and before interchange situations.(2)If is the first customer on the vehicle, then “b” represents the total difference of transportation costs between after interchange and before interchange situations.(3)Either “a” or “b” must be equal to 0.(4)If precedes , then “c” represents the differences of transportation costs between after interchange and before interchange situations.(5)If the sum of “a,” “b,” and “c” is smaller than 0, it means that a lower total transportation cost is obtained. In this case, and are interchanged.(6)Calibrate the random key numbers. In Section 3.2.2, item (i) gives the description of the intervehicle type interchange.(II)Nonadjacent interchange: the steps for determining the exchange of and on the route are similar to adjacent interchange with one exception. During the calculation of “a” and “b,” the additional difference of transportation costs should be considered. The remaining steps are the same with the adjacent interchange situation.(ii)Intravehicle type exchanging: in the intravehicle type exchanging section of the pseudocode, steps are as follows.(1)Calculate the total transportation cost of and before exchange. “D1" represents the maximum transportation cost of them.(2)If precedes on , then “a1" represents the difference of transportation costs between after interchange and before interchange situations. If is the first customer on , the second “a1" in the pseudocode represents the total difference of transportation costs between after interchange and before interchange situations.(3)If precedes on , then “b1” represents the difference of transportation costs between after interchange and before interchange situations. If is the first customer on , the second “b1” in the pseudocode represents the total difference of transportation costs between after interchange and before interchange situations.(4)If precedes on , then “b2” represents the difference of transportation costs between after interchange and before interchange situations.(5)Calculate the total transportation cost of and after a possible exchange. “D2” represents the maximum of these two values.(6)If “D2” is smaller than “D1”, it means a lower total transportation cost is obtained. In this case, and are interchanged.(7)Calibrate the random key numbers. In Section 3.2.2, item (ii) gives the description of intravehicle type exchange.
3.2.2. Calibration of the Random Numbers on the Chromosome
If it is decided to change and , we need to calibrate the chromosome to represent the new situation by the following configurations. During the integration process, the positions of the random numbers are considered. In Figure 4, parent 2 is taken into account to explain the calibration phase. The orders of and are split into 2 suborders (no split for the orders of ) for the considered chromosome. Due to simplicity, the sections of the number of suborders are not shown on the chromosomes of Figures 7 and 8.(i)Intervehicle type interchange: Figure 7 shows the before and after interchange of the customers on the same vehicle type. The genes of a chromosome can be labeled as “selected genes from the chromosome” and “nonselected genes from the chromosome.” The realized interchanges are shown in bold on the after interchange section. The steps of interchanging and for the example are as follows.(I)Random numbers of the customers to be exchanged are 0.40 for and 0.62 for to be changed to 0.62 for , and 0.40 for .(II)After the interchange, 0.40 is smaller than 0.62, and hence there is no need to change other random numbers for . Because, the current situation guarantees to select .(III)After the interchange, 0.62 is bigger than 0.40, and hence other random numbers for , which are smaller than 0.62 (as 0.41) in the “nonselected genes from chromosome,” must be changed. A new random number (bigger than 0.62) will be randomly generated (as 0.96) and changed to 0.41. Another random number, smaller than 0.62, is 0.61 in the “nonselected genes from chromosome” for . Therefore, a new random number (bigger than 0.62) will be randomly generated (as 0.75) and changed to 0.61.(ii)Intravehicle type exchange: Figure 8 shows the before and after exchange of the customers on two different vehicle types. The realized exchanges are shown in bold on the after exchange section of Figure 8. The steps of exchanging and for the example of Figure 8 are as follows. Value is a maximum random number value for on . Similarly, Value is a maximum random number value for on .(I)The random numbers for and Value are exchanged in the chromosome. Also, the random numbers for and Value are exchanged in the chromosome.(II)Between “selected genes from chromosome” and “nonselected genes from chromosome,” customers for and Value are exchanged. And similarly customers for and Value are exchanged.(III)After the exchange, it must be checked if there is a random number in the “nonselected genes from chromosome” for that is smaller than the random number of (equal to 0.84) or other selected value of (equal to 0.34). For example, 0.41 is smaller than 0.84. Therefore, another random number (bigger than 0.84) will be randomly generated (as 0.95) and changed to 0.41.(IV)It must also be checked if there is a random number in the “nonselected genes from chromosome” for that is smaller than the random number of (equal to 0.40) or other selected value of (equal to 0.62). For example, 0.61 is smaller than 0.62. Thus, another random number (bigger than 0.62) will be randomly generated (as 0.87) and changed to 0.61.
As mentioned before, we consider a realworld MCHF/OVRP/SDMP problem. We first model the real system with mixedinteger programming, and then we develop a hybrid geneticlocal search algorithm. We determine the required parameters and their values for the algorithms. We compare these algorithms with randomly generated problems using realworld parameters and detail the computational results in the next section.
4. Computational Results
4.1. Dataset
We consider the distribution network of a production company in Turkey. The company produces two types of products and uses a fleet of vehicles of a 3PL provider for distribution of goods to customers. In order to evaluate the performance of MIP and GA, we use the realworld data which are supplied from the 3PL provider. We gather the information of transportation costs, vehicle types and specialties, and properties of products from the 3PL provider, which are presented in Tables 5, 6, and 7, respectively.



We investigate the realworld demands and determine the distribution of them for two products as seen in Table 8. We represent the connections between customers with a 01 “Route ()” matrix as in MIP model. The matrix is defined by the 3PL provider considering various constrains such as highway conditions or experiences of drivers. To generate different problems, we determine the matrix randomly for each problem. We develop matrix for customers and define the cells by distribution. If the value of customer and customer is 1, it means that there is a connection between the customer and customer . We also correct the diagonal of the matrix which should be 1. The generated examples can be found in [18].

4.2. Experimental Results
We generate different sizes of problems with different number of customers (varies from 10 to 90) which can be found in [18]. All experiments are run on a personal computer with an Intel Core i73612QM processor running at 2.10 GHz. We use C# for the GA implementation and CPLEX 12.4 for the MIP model. Before comparing MIP and GA, we first analyze the performance of MIP with different time limits. We limit MIP with 1800, 3600, and 7200 seconds and represent the results as MIP_18, MIP_36, and MIP_72, respectively. In the realworld application, the company has more than 50 customers. Although we prepare samples up to 90 customers, we could not use all of them since MIP_7200 could not find any feasible solution within time limits for the problem with 60 customers. Therefore, we restrict Table 9 with the samples up to 50 customers. We increase the number of samples for smallsize problem, but the gap% gets bigger after 30 customers and problems get harder to solve; thus we generate only one problem for more than 30 customers. The gap% for the MIP models is calculated by the branch and bound algorithm in CPLEX 12.4 as in (11), where the best solution is represented by and the best integer solution is represented by . The experimental results are detailed in Table 9, The first part of Table 9 gives the information about samples, the second part presents the objective function values found within different time limits, the third part shows the required time for solutions, and the fourth part indicates the gap% values which are calculated by CPLEX 12.4. The last part shows the improvements in percentage between time limits; the first column (18–36) presents the improvement in total cost if we run the MIP model within 3600 seconds instead of 1800 seconds, the second column (18–72) presents the improvement in total cost if we run the MIP model within 7200 seconds instead of 1800 seconds, and the last column (36–72) presents the improvement in total cost if we run the MIP model within 7200 seconds instead of 3600 seconds.

Table 9 indicates that the MIP model solves the smallsize problems (10 customers) optimally within minutes. When the number of customers increases to 15, MIP generates similar results with all time limits, but it could not prove the optimality of solutions within time limits. Considering the problems that contain 20–50 customers, the difference between time limits becomes clear, and MIP_7200 performs better than the other compared options. In addition, MIP_7200 reaches solutions with smaller gap which are close to optimality as seen in the gap% part. Although the solutions look similar with different time limits, the improvement part clarifies the contribution of 7200 time limits to model results. If we use 3600 seconds instead of 1800 seconds, we improve the results with an average of 1.15%. However, if we use 7200 seconds, we improve the results with an average of 2.38%. Moreover, the realworld applications usually contain more than 50 customers, and the MIP_7200 model would provide better feasible results for larger examples; thus we limit MIP within 7200 seconds for next comparisons.
We compare the MIP model with the GA by using same randomly generated samples in previous experiment [18]. We limit MIP within 7200 seconds and determine the generation limits for GA as 500 generations. We detail the comparison between MIP and GA in Table 10. The first column shows the number of customers in each problem, the second column indicates the number of problems, the third part demonstrates the objective/fitness function of the MIP and GA, and the fourth part indicates the solution time of algorithms. We calculate the gap% for GA similar with the branch and bound algorithm by using the MIP result as the best solution in (12). The last two parts demonstrate the number of vehicles used by the algorithms. We define an upper bound for the vehicles in order to limit the vehicle index in the MIP model; however, the GA does not need such an upper bound:
 
MIP could not find any feasible solution within time limits. **Gap% values for GA could not be calculated because of the deficiency of MIP cells. 
Experimental results show that MIP is able to solve the problems with 10 customers to optimality. However, for the problems that contain 20–50 customers, MIP could not reach the optimal solutions within time limits (2 hours). Moreover, MIP could not provide any feasible solution in 2 hours for largesize problems that contain 60–90 customers. On the other hand, GA is able to reach similar results with a 9.66% gap. In addition, the GA provides feasible solutions for larger problems which have similar size with realworld applications. The last three results of GA’s elapsed time are higher than 7200 seconds, which are shown in Italic. It is possible to decrease the elapsed time via decreasing the number of generations. For example, if the same algorithm would run with 250 generations instead of 500 generations for the problem of 70 customers, the elapsed time would be 3817.24 seconds and total cost would be 85.660. Therefore, we can change the parameters of GA in order to catch the flexibility of realworld applications, and we can arrange them to reach feasible solutions within desirable times. In conclusion, the number of customers increases, the problem becomes harder to solve for both MIP and GA. However for the realworld applications, GA is preferable rather than MIP to reach feasible solutions in short time periods, by arranging the algorithm’s parameters.
5. Conclusions and Directions for Future Research
In this paper, we introduce a vehicle routing problem, the MCHF/OVRP/SDMP, which has not yet been investigated in detail so far. With the aim of minimizing the total distribution cost, we have formulated a MIP model for the problem and developed a solution approach based on a geneticlocal search algorithm requiring much less computational time than the MIP. Since the problem type discussed in this paper is not exactly the same as those in the literature, we randomly generated different problems with different routes and demands. We used these problems to evaluate the MIP model and the GA. The experimental results showed that although MIP was able to find better solutions than GA for smallsize problems, it was not possible to obtain any feasible solutions of largesize problems within time limits. On the other hand, GA was able to reach feasible solutions of largesize problems in shorter time periods. However, the gap% between GA and MIP is greater than expected, and therefore, we continue our research to improve GA and to investigate new solution methods for this unique vehicle routing problem. Generating hybrid algorithms through the use of different metaheuristics might be a worthwhile avenue of research.
Appendix
The pseudocode for the local search: for brevity, only the key steps are detailed. C# notation is used see Algorithm 1.
