Abstract

We present the vehicle routing problem with potential demands and time windows (VRP-PDTW), which is a variation of the classical VRP. A homogenous fleet of vehicles originated in a central depot serves customers with soft time windows and deliveries from/to their locations, and split delivery is considered. Also, besides the initial demand in the order contract, the potential demand caused by conformity consuming behavior is also integrated and modeled in our problem. The objective of minimizing the cost traveled by the vehicles and penalized cost due to violating time windows is then constructed. We propose a heuristics-based parthenogenetic algorithm (HPGA) for successfully solving optimal solutions to the problem, in which heuristics is introduced to generate the initial solution. Computational experiments are reported for instances and the proposed algorithm is compared with genetic algorithm (GA) and heuristics-based genetic algorithm (HGA) from the literature. The comparison results show that our algorithm is quite competitive by considering the quality of solutions and computation time.

1. Introduction

The vehicle routing problem (VRP) is a well-known combinatorial optimization problem arising in transportation logistic, which was first introduced by Dantzig and Ramser in their study on the truck dispatching problem [1]. The classical VRP (CVRP) is generally defined by a graph , where is the set of vertices comprised of a central depot and customers and is the arc set. The CVRP aims to determine a set of routes for identical and capacitated vehicles by minimizing the total routing cost (distance, time, etc.), such that each customer is visited exactly once. It mainly focuses on deterministic demand, which means customer demands are fixed and known in advance, and there is no requirement for the service time.

Numerous variants of CVRP have been put forward and extensively studied which capture the complications of real-world problems [24], for example, VRP with time windows (VRPTW) [521], in which customers have to be visited within a predefined time interval which can be described by the earliest arrival time and the latest arrival; VRP with multiple depots (VRPMD) [20, 2225]; VRP with simultaneous delivery and pickup (VRPSDP) [10, 26], where customers require simultaneous pickup and delivery service; VRP with stochastic demands (VRPSD) [19, 2735], which involves demands that are random variables with known distributions; VRP with split delivery, where the demand of a customer can be split and delivered by multiple vehicles [36]; and the heterogeneous fleet VRP with several vehicle types [37].

Besides, lots of surveys were dedicated to studying solution methods of VRPs which can be classified into exact algorithms, heuristics, and metaheuristics [3, 4, 22]. We refer the interested reader to the most common exact algorithms, namely, a branch-and-cut by Hà et al. [38]; a branch-and-price by Christiansen and Lysgaard [28] and Hadjar and Soumis [5]; a branch-cut-and-price by Gauvin et al. [29] and Lysgaard and Wøhlk [39]; Lagrangian relaxation by Kallehauge et al. [6]; Dantzig-Wolfe decomposition by Qureshi et al. [7] and Chabrier [40]; and a branch-and-bound by Carpaneto et al. [23]. However, exact algorithms are usually designed for the static and deterministic VRPs, where all input is known and keeps unchanged once vehicles are in service. Since VRP is a typical NP-hard problem, many heuristics-based studies are proposed for various VRPs. The work of Pang [8] and Figliozzi [9] applied the improved route construction heuristics to solve VRPTM, and Liu et al. [24] and Z. Wang and Z. Wang [41] raised the novel heuristics based on the traditional two-phase method and the results demonstrated the effectiveness of the proposed heuristics. Most of the recent researches for VRPs paid extensive attention to the development of metaheuristics. This can be summarized as (1) genetic algorithms, (2) Ant Colony Optimization (ACO), (3) simulated annealing (SA), (4) Variable Neighborhood Search (VNS), (5) tabu search (TS), (6) local search algorithm, (7) Artificial Bee Colony (ABC), and (8) Particle Swarm Optimization (PSO). For further details on VRP metaheuristics, we refer the reader to the work of [3, 22, 42].

Conformity behavior was firstly introduced by M. Sherif, a social psychologist, and proposed to study the how peers’ behavior affects an individual’s behavior in a group [43]. Many studies have focused on herd behavior [44, 45]. For example, Jiang and Du analyzed the homebuyers’ conformity behavior [46] and Cipriani and Guarino studied herd behavior in financial market [47]. However, there is currently very little methodology in the literature to model and optimize the vehicle routing problem with conformity behavior. Our VRP is motivated by a real-life problem, where the depot receives orders from a set of customers and then assigns vehicles to deliver products after signing the contract. The customer demand is generally confirmed in the contract, which is confirmed as the initial demand. However, once a customer knows the initial demand of other customers, they could make some adjustments for their initial demand due to conformity behavior. Thus, customer’s ordering behavior is prone to be affected by other customers [47], which is named as conformity consuming behavior in this paper. This may lead to generation of new orders in the delivery, which is consistent with the rush order; the reduction or canceling of the initial demand is also possible. In the following sections, we call it the potential demand to differ it from the initial demand in the contract. Once the potential demand occurs, the manager has to reschedule the routes of vehicles. Henceforth, how to assign vehicles and give the routing plan is an issue that needs to be solved, which is what our study concerned. To the best of our knowledge, there is no literature devoted to the VRP with potential demands and time windows.

Our paper has three main contributions, which are as follows: the first is the establishment of the vehicle routing problem with potential demands and time windows, in which the concept of conformity consuming behavior is taken into consideration and the route of vehicles is constructed by multiple cycles; the second is the use of parthenogenetic algorithm for saving the computation time caused by searching/repairing the invalid solutions after crossover operator in GA, in which a novel coding technique using 2D vector is proposed to address our problem; and then heuristics is designed for the generation of the initial population for the parthenogenetic algorithm (PGA) such that a hybrid algorithm is proposed, which is our third contribution. The remaining outline of this work is organized as follows. Section 2 reviews the literature; Section 3 introduces the modeling assumptions and the mathematic model of the VRP-PDTW; in Sections 4 and 5, HPGA by combining parthenogenetic algorithm with heuristics is designed and developed to solve this problem; this is followed by computational results and comparison in Section 6; Section 7 concludes the paper and the further research is also given.

2. Literature Review

In the deterministic capacitated vehicle routing problem (DCVPR), as the variation of the CVRP, the demand of each customer is deterministic and there is no limit of time windows for visiting customers [39]. The objective in the work of Lysgaard and Wøhlk is to minimize the sum of arrival times at the customer which is referred to as cumulative VRP, and then a set partitioning formulation and a vehicle flow formulation are proposed and combined to apply into a branch-and-cut-and-price algorithm [39]. De Oliveira et al. considered a multidepot DCVPR (MD-DCVRP) where the maximum route duration time is set, and a coevolutionary approach is introduced through decomposing the MD-DCVRP into a number of single-depot DCVRPs [25]. Tasan and Gen developed a genetic algorithm based approach for the DCVRP with simultaneous delivery and pickup (DCVRPSDP) [26]. Wang and Chen modeled for the DCVRPSDP with time windows and in order to solve this problem, the cheapest insertion method is introduced into the genetic algorithm to speed up the solution procedure [10]. Kallehauge [11], Baños et al. [12], Qi et al. [13], Yu et al. [14], Vidal et al. [15], and Koç et al. [16] considered the case of the DCVPR with hard time windows, in which the vehicles have to arrive at the customers before the latest arrival time; otherwise, the vehicles may wait at a customer location if they arrive earlier than the earliest arrival time. In the literature [1416], neighborhood search is also introduced to enhance the local search ability. A Pareto-based hybrid algorithm that combines evolutionary computation and simulated annealing is reported by Baños et al. for solving the multiobjective formulations of the DCVRPTW [12]. The DCVPR with soft time windows is investigated by Ghoseiri and Ghannadpour [17], in which the immediate service from late arrival vehicles is accepted by customers but a penalized cost is imposed on the latter arrival vehicles. An adaptive genetic algorithm incorporating push-forward insertion heuristic and -interchange mechanism is then constructed to search for better routing solutions. Related works in the soft time windows area are also proposed by Lau et al. [18] and Qureshi et al. [7], in which an exact algorithm arising from the Dantzig-Wolfe decomposition method and a tabu search approach are presented, respectively.

However, in most real applications, one or more parameters (time, demand, etc.) of the CVRP tend to be random or stochastic, giving rise to the stochastic vehicle routing problem (SVRP). We found that most researches mainly focused on VRPSD, where customer demands are known only when the vehicle arrives to location [2835]. In the following, some studies will be mentioned. Christiansen and Lysgaard developed a branch-and-price algorithm [28] for the minimal distribution cost. Marinakis et al. [30] proposed a new hybrid algorithmic approach, which combined the PSO algorithm with the 2-opt and 3-opt local search algorithms and with the path relinking strategy. A paired cooperative optimization strategy is developed for the VRPSD with split delivery [31]. But, the constraint for the visiting time of customers is not considered in these studies [2834]. However, to be closer to the practice, Tan et al. [35] and Lei et al. [19] applied a soft time windows constraint for the VRPSD, which means that vehicles may be violated if the arrival time at the location of customers is beyond time windows. Tan et al. [35] defined the time constrained VRPSD as a multiobjective optimization problem such that an evolutionary algorithm with enhanced genetic operators, Pareto fitness ranking, and local search heuristics is presented. The simulation results illustrated that the solutions obtained from the proposed algorithm are robust to the stochastic nature of the problem. The VRPSD is modeled as a stochastic program with recourse by Lei et al. and for its solution an adaptive large neighborhood search heuristic is designed, in which several remove-insert operations are combined to efficiently explore the solution space [19].

There are some researches on the fuzzy VRP (FVRP) which arises whenever some elements of the problem are uncertain, subjective, and vague. Common examples are fuzzy travel times and fuzzy demands [20, 21, 48]. Generally, the fuzzy variable is firstly handled using the possibility measure theory. For example, the fuzzy demand on each customer and the travel times between nodes and time windows are often described as triangular fuzzy numbers. Tang et al. considered the vehicle routing problem with fuzzy time windows, in which the deviation of the servicing time from the customer’s time window is assumed to be associated with the level of customer’s satisfaction, and proposed a two-stage algorithm to obtain a Pareto solution [20]. Asl et al. extended work [20] to multidepot and multilevel VRP with the fuzzy time windows and sequentially a three-section algorithm is proposed [21]. Cao and Lai designed a fuzzy chance constrained program model based on fuzzy credibility theory and a hybrid intelligent algorithm with stochastic simulation and differential evolution algorithm is presented to solve this problem [48].

3. Model Formulation

3.1. Problem Description and Assumptions

We focused on the vehicle routing problem with a central depot, denoted by “,” and a set of customers to be serviced, which is defined by a directed graph . is the node set including the central depot “” and customers and is the set of arcs , , . Arc represents the possibility of traveling from to with an associated distance, duration, or cost. The customers are serviced by homogeneous vehicles with a limited capacity and the same speed, and all vehicles are located at the central depot such that they have to depart from the depot and return to the depot after visiting some customers. In the following, some assumptions are given for the modeling.(1)Consider a single type of products stored at the central depot, which is delivered by vehicles and ordered by customers.(2)The duration for a vehicle from departing from depot to ending at depot after going through several customers is defined as a cycle.(3)The split delivery for the customer demands is allowed; thus each customer may be visited more than once and serviced by different vehicles. This implies that the route of a vehicle may be subject to multiple cycles. Taking the route of a vehicle “” as an example, there are two cycles “” and “.” Particularly, the number of the cycle in a route of vehicle is represented by the maximum delivery frequency in our following model.(4)The initial demand of customer , is , but his/her potential demand may occur before a vehicle visits him/her for the first time. relies on the difference between the initial demand of customer and that of other customers’ . Thus, the potential demand is defined as , in which it has a wide choice of function forms for including linear function, logarithmic function, and semilogarithmic function.(5)Time windows , in which we let denote the earliest arrival time and let be the latest arrival time that customer can be serviced by a vehicle, are imposed on customer . These are soft bounds that can be violated, so the penalty cost incurred by waiting for servicing and late servicing at customers arises. However, there is no constraint of time windows for the central depot.(6)The service time at node is , and the travel time from customer to customer with the travel cost per unit time is associated with each arc in the directed graph.(7)The product is always available for the central depot, that is, no shortage.

3.2. Mathematical Model

The objective function of the proposed VRP-PDTW is to minimize the total cost including the travel cost caused by the delivery of vehicles and the penalty cost due to violating the time windows, which is given by the following:where and are set of all vehicles and the number of cycles for the route of vehicle , respectively. Let denote the arrival time of a vehicle in the th cycle at customer , and returns the maximum value from and . There are two decision variables, that is, and . if vehicle travels from customer to customer in the th cycle; otherwise, ; and is the demand of customer which is satisfied by vehicle in the th cycle.

In the objective function (1), and are the penalty cost per unit time as related to violating time windows, respectively. Constraints (2) and (3) represent that a vehicle starts its route at the depot, visits a number of customers, and returns to the depot in the th cycle. Constraint (4) imposes the notion that, in the th cycle, a vehicle visits customer only once. The sequence of a route “” for vehicle in the th cycle is imposed by constraint (5). The demand of all customers met by the vehicle in the th cycle does not exceed the maximum capacity, as mentioned in constraint (6). Constraint (7) indicates that both the initial and potential demands of customer are satisfied. The time constraint is controlled by (8) to ensure the time of a vehicle in the th cycle reached to customer must be not earlier than , in which the equality holds if fails to the interval ; otherwise, .

The VRP with potential demands and time windows is a typical problem in combinatorial optimization, which is also an NP-hard problem [4953]. It is theoretically and practically impossible to use exact methods for the large-scale VRP; thus a hybrid parthenogenetic algorithm mixed with heuristics is designed in our view to find approximate solution, which is described in the following section.

4. Framework of the Heuristics-Based Parthenogenetic Algorithm

The genetic algorithm based on the drawing concept of evolution, survival of the fittest, which was proposed by Holland in 1975, has been extensively adapted to solve the NP-hard problems. It and the improved algorithms based on GA have been widely implemented on VRP and its variants over the past decades (see [10, 1517, 26]). In GA, the candidate solutions of solving problem are represented as a population of chromosomes (individuals) and each chromosome is comprised of a string of genes. The search starts from the initial population, which is either randomly generated or generated using heuristics, and runs in generations where the individuals evaluated by the fitness function keep improving by applying genetic operators. Commonly, genetic algorithm iterates over generations till a predefined stop criterion is met. The steps of the genetic algorithm are summarized as follows.

Step 1 (problem representation). In GA the problem is in the form of a genotype representation by using a chromosome regarding a solution. The genes in chromosomes can be encoded in binary for bit string representation or represented as a string of real numbers.

Step 2 (initialization). Randomly generate the initial population of solutions with the effective population size to cover the search space.

Step 3 (fitness evaluation). The fitness value of each chromosome in the population is calculated and ranked by the fitness function.

Step 4 (genetic operators). Evolution mechanisms subject to selection, crossover, and mutation are used to create fitter individuals for the next generation.

Step 5. Repeat Steps 3 and 4 until terminating condition is met, in which the stop criterion can be either the maximum number of generations or a satisfactory fitness level.

The crossover operator, as the major genetic operator, is used to generate new offspring by mixing two parents such that some characteristics from the parents can be inherited by newer individuals. The traditional crossover operators, that is, one-point crossover, two-point crossover, scattered crossover, and a partially mapped crossover, can be found in the VRPs, and some improved crossover operators are also proposed depending on the various problems [3, 4]. The crossover operator is intended to generate solutions in the whole search space. However, how to design a feasible crossover operator for a special VRP is still a key challenge since invalid offspring could be obtained through crossover and it is time-consuming to repair them.

Parthenogenetic algorithm is fortunately proposed [54] to deal with the above issue by the removal of the crossover operator, which greatly simplifies the procedure of GA and improves the effectiveness and performance of searching optimal solutions. This is because in PGA there is the shift operator only performed on a single chromosome. It avoids the offspring obtained by crossover operator jump to the invalid solution area. As we mentioned in the previous section, that is, (1), it is impossible or difficult to create the initial solution randomly that constitutes the initial population of our problem in PGA. A heuristic algorithm is therefore designed to construct the initial population, which will be shown in the following section. Figure 1 shows the framework of the heuristics-based parthenogenetic algorithm.

5. Design of the Heuristics-Based Parthenogenetic Algorithm

5.1. Heuristic for Initial Solutions

From the problem description mentioned in Section 3, we design a heuristic algorithm for the generation of initial solutions, as described in Figure 2, in which two rules are firstly developed to perform search. For more details, the procedure can be found in Appendix.

Rule 1. For a set of customers with the initial demand , a vehicle chooses the nearest customer to provide service. However, the customer with the largest demand is visited by vehicle prior to others once there are two or more customers with the same traveling distance.

Rule 2. The initial demand for each customer is always considered to be preferentially satisfied against his potential demand.

The quality of the initial population size is also a key issue since it directly affects the performance and efficiency of HPGA. It is commonly preset as an integer depending on the problem to be solved. Without exception, we give the initial population size as an integer which is associated with the number of vehicles and customers, that is, .

5.2. Solution Coding

The important step for HPGA is to determine the solution coding technique to represent chromosome. Binary representation, permutation representation, integer representation, and sector representation are commonly applied in the literature on VRPs. Assume that if we have obtained an initial solution for our problem, which is described in Table 1, it is impossible to code it using the coding approaches mentioned above. It is because for a vehicle there may be multiple cycles, rather than delivery only once. So, we propose here a novel permutation representation based on a 2D vector , where and are the indices of vehicle and node, respectively, such that is defined as the vehicle visits node . The generalizable chromosome representation is then described in Figure 3.

It can be seen from Figure 3 that each vehicle leaves the central depot 0 and comes back to the depot through some customers after completion of its delivery. Using the solution coding technique in Figure 3, the initial solution in Table 1 can be expressed asObviously, there are multiple cycles for vehicles 1 and 3.

5.3. Shift Operator

Different from genetic operators in GA, shift operator is usually applied to generate offspring in PGA, as defined in the following.

Definition 1. Shift operator exchanges some genes with a certain probability and works on a chromosome, in which the location of genes exchanged is randomly chosen.

Considering the length difference of each cycle in a vehicle route, that is, the number of customers visited by vehicle in a cycle, different scenarios could be applied based on the int  -point exchange policy, in which int  rounds a number down to a nearest integer. It implies that there are two different shift operators subject to single-point exchange and multipoint exchange, as shown in Figures 4(a) and 4(b), respectively. Single-point exchange corresponds to the location of a pair of genes that mutually interchange, referred to as genes B and D in Figure 4(a). However, multipoint exchange is defined as two or more pairs of genes are interchanged. From Figure 4(b), it can be obviously seen that there are two pairs of genes exchanged, that is, and .

Shift operator in Definition 1 is utilized to produce newer individual, but in our context, it is restricted to a cycle of each vehicle. This is because exchange across cycles will result in invalid individuals and cost more repair time. For the chromosome in Table 1, Figure 5 gives the process of shift operator that is used to the route of vehicle 2. Similarly, it is also implemented on the routes of vehicles 1 and 3.

5.4. Evaluation and Selection

Using the objective function shown in (1), the fitness is defined as the reciprocal of TC; that is, . Thus, the fitness value of each individual from the initial population and its offspring can be calculated to scale the performance of each individual. The roulette wheel selection is used to select fitter individuals, as defined in (10). It implies that the probability of individuals with lower cost to be selected for the next generation is higher. Considerwhere is the number of individuals to be selected.

6. Experimental Results and Comparisons

We coded the above procedure using MATLAB R2012b and ran it on a PC with 3.60 GHz CPU and 1 T memory, and all computing times are reported in CPU seconds. This section presents comparison results and analyses for results carried out to investigate the effectiveness of the proposed HPGA through a number of computational experiments. We first apply the proposed HPGA to a real-world instance of a distribution company in Handan to verify its applicability in real life, and then simulation experiments are conducted to assess the performance of HPGA by comparing with the GA and HGA.

6.1. Parameter Settings

There are 12 customer nodes and they are serviced by 4 homogeneous vehicles which are located at the central depot and the vehicle capacity is 80 units. The traveling time between customers is set in terms of the corresponding Euclidean distance between them, and some parameters values used in the mathematical model are also given in Table 2. The potential demand in Table 2 is obtained by simplifying assumption (4), in which only one related customer of is considered in this simulation and a positive/negative number indicates the increase/decrease of the initial demand, respectively. The following cost parameter values are used in our experiments: , , and .

Also, the proposed algorithm HPGA is applied with population size , where we let denote the number of customers with nonzero initial demand, that is, . In our experimental studies, the run is terminated if the fitness is not improved; that is, the solution value is within a preset percentage of the optimum . The shift probability on shift operator is set to 0.7 by a pilot study.

6.2. Computation Results

Depending on the parameters in Section 6.1 and the proposed algorithm, the optimal routes and work time for vehicles are obtained, as shown in Table 3. More details on the arrival/departure time (AT, DT) of vehicles at each node are given in Table 4, and the corresponding quantity unloaded at nodes is also output. From Tables 3 and 4, it can be seen that vehicles 1 and 3 firstly end their first cycle at 54 and 58, respectively, and then start their second cycle to service customers after uploading products at the central depot “.” Furthermore, Figure 6 shows the best objective function value through the evolution process in this instance, which reveals the optimal solution yields when the proposed algorithm runs for 250 CPU seconds approximately. However, it should be noted that different results could be obtained for other parameters or the regeneration of potential demands.

6.3. Comparative Analyses

We here perform comparative analyses of the results of the HPGA with the GA and HGA to evaluate the efficiency between them, using the previous modeling parameter settings. Heuristics proposed in Section 5.1 is further integrated into the GA to create the initial solution. We expect HGA can save computational time since it is not easy to obtain feasible solution by random generation for our VRP with potential demands and time windows. Note that in GA and heuristics-based GA some initial parameters including population size, the crossover and mutation probabilities, and the termination condition should be first set. To remove any factor that could impact the performance, the population size and the termination condition are the same in GA, heuristics-based GA, and HPGA. Since what is concerned is the effect of the crossover probability in GA and HGA on the performance of solution, the crossover probability varies from 0 to 1 with a step of 0.2. However, the mutation probability is set to 0.3 in our instances by using the proposed method by Greenwell et al. [55].

The comparison results are summarized in Table 5, where the proposed algorithm is tested in six instances and in the first column each instance is numbered. These instances are available by randomly generating potential demands, keeping other parameters unchanged. The efficiency of the algorithm is measured by the best objective function value and the computation time, which are denoted by Obj and Time in Table 5, respectively. Also, the gap is calculated in terms of the relative deviation from the best known solution to measure the performance of solutions with the following equation:where Obj# is the value of one of the other methods (# is either GA or HGA) and is the solution value of HPGA. In Table 5, the optimal solutions for each instance are denoted with bold face letters, and the corresponding gap values are also marked in bold.

As Table 5 demonstrates, our proposed algorithm HPGA is more effective and very competitive compared to GA and HGA. It is noted that, for each instance, there is only one best solution in HPGA. This can be explained as the crossover operator is not used in HPGA, so the change of does not work on HPGA. We can see that the computational time of GA and HGA decreases with the increase of , which results from the fact that a larger speeds up the convergence. For the scenario in instance 2, the computation time of GA is represented as “—,” which implies that it runs beyond the preset time threshold, 3 hours. Obviously, HPGA searches the optimal solution with the smallest time while GA costs the longest time, HGA between them. This is probably because heuristics improves the search speed of HGA and HPGA by generating the reasonable initial solutions as soon as possible. It may be also because the crossover operator in GA and HGA increases the time of searching the optimal solution and repairing the invalid solution.

In Table 5, it can be concluded from the results of Obj and gap that the quality of our proposed algorithm outperforms those of GA and HGA, which is what we expected. The gap values obtained from these instances are positive, varying between 0.31 and 0.02. This provides full and strong support to illustrate the solution quality of HPGA. Additionally, it can be found that the optimal solution for different instances is obtained with different crossover probabilities, which is out of the scope of this study. However, the proposed algorithm is not proved to be the most effective compared to all other algorithms ever published in the literature.

7. Conclusions

In this paper, we investigated the vehicle routing problem with potential demands and time windows inspired by a real-world problem, which considered the initial demand from customers and the possible potential demand induced by conformity consuming behavior. A fleet of vehicles originated in a central depot serves customers and the split delivery is allowed in our problem. The VRP-PDTW is formulated as a single objective model to minimize the total cost and find the optimal routes of vehicles. The VRP-PDTW is also an NP-hard problem; thus it is difficult and time-consuming to obtain the solutions using exact algorithms. Then, a heuristics-based parthenogenetic algorithm for solving VRP-PDTW was proposed in this study, in which the initial solution can be created using heuristics. The shift operator is only performed on a chromosome in our proposed algorithm, which differs from that in GA. Also, the coding technique based on a 2D vector is presented and applied in our algorithm. Finally, some experiments are carried out to evaluate the performance of the proposed algorithm. The results show that our proposed algorithm produces competitive solutions and is the fastest when compared with the best results of GA and HGA. The benefit of HPGA comes from its ability to produce the initial solution by heuristics and avoid the occurrence of invalid offspring which is generated by the crossover operator.

Our further studies will be focused on other real-world routing problems, for example, considering multiple types of products, the vehicle availability and the shock from environment; and relaxing the fixed number of routes will be further optimized.

Appendix

See Algorithm 1.

Step  1: Set the model and algorithm parameters , , , and initialize some transient variables
for heuristics , , , , , .
 We let if customer i has not been visited by vehicles yet, otherwise, ; is defined as the remaining quantity
of products loaded by a vehicle k when it departs from customer in its lth cycle; is the cost due to the time windows
is violated, that is, the second term of the objective function, (1).
Step  2: Initially, each vehicle starts from the central depot “” for delivery, and update and for simplicity
where and ;
Step  3: Vehicle k is assigned to service customer j in terms of Rule 1, and update and ;
If  
   Generate the potential demand of customer j, update , and go to Step  4;
Else go to Step  4;
End;
Step  4: Calculate the decision variable in terms of Rule 2, that is, the demand of customer j met by a vehicle in its lth cycle,
and the penalty cost induced by the time windows.
 (a) If   
   If    then and update
         , , ;
      Else    and update , ;
      End;
   Else
    If    then and update , ;
    Else    and update , ;
    End;
   End;
 (b) It can be concluded that from assumption (5):
   If    then ,;
   Elseif then ;
   Else   , ;
   End;
Step  5: Decision whether to continue to delivery or return to the central depot.
If  
   , , go to Step  6;
Else
    If    then , and update , ;
     Go to Step  3;
   Else go to Step  3;
   End;
End;
Step  6: Output the route of all vehicles according to , and the decision variable such that a solution for VRP is created.
Step  7: Repeat Steps  1–6 and terminate until M individuals are generated.

Competing Interests

The authors declare that they have no competing interests.

Acknowledgments

This research is partially supported by the Natural Science Foundation of Hebei Province (no. G2014402027), Science and Technology Project of Hebei Province (nos. 15454707D and 16457402D), Humanities and Social Sciences Research Project for Universities in Hebei Province (no. 2014068), Three Modernization and Collaborative Development Base Project of Hebei Province (no. SH2015QR01), and High Level Talent Support Project of Hebei Province (no. A201500112).