J. Cent. South Univ. (2012) 19: 433-442

DOI: 10.1007/s11771-012-1022-5

Solving material distribution routing problem in mixed manufacturing systems with a hybrid multi-objective evolutionary algorithm

GAO Gui-bing(高贵兵)^{1, 2}, ZHANG Guo-jun(张国军)¹,

HUANG Gang(黄刚)¹, ZHU Hai-ping(朱海平)¹, GU Pei-hua(顾佩华)^{1, 3}

1. State Key Laboratory of Digital Manufacturing Equipment and Technology,

Huazhong University of Science and Technology, Wuhan 430074, China;

2. School of Energy and Safety Engineering, Hunan University of Technology, Xiangtan 411201, China;

3. College of Engineering, Shantou University, Shantou 515063, China

? Central South University Press and Springer-Verlag Berlin Heidelberg 2012

Abstract: The material distribution routing problem in the manufacturing system is a complex combinatorial optimization problem and its main task is to deliver materials to the working stations with low cost and high efficiency. A multi-objective model was presented for the material distribution routing problem in mixed manufacturing systems, and it was solved by a hybrid multi-objective evolutionary algorithm (HMOEA). The characteristics of the HMOEA are as follows: 1) A route pool is employed to preserve the best routes for the population initiation; 2) A specialized best-worst route crossover (BWRC) mode is designed to perform the crossover operators for selecting the best route from Chromosomes 1 to exchange with the worst one in Chromosomes 2, so that the better genes are inherited to the offspring; 3) A route swap mode is used to perform the mutation for improving the convergence speed and preserving the better gene; 4) Local heuristics search methods are applied in this algorithm. Computational study of a practical case shows that the proposed algorithm can decrease the total travel distance by 51.66%, enhance the average vehicle load rate by 37.85%, cut down 15 routes and reduce a deliver vehicle. The convergence speed of HMOEA is faster than that of famous NSGA-II.

Key words: material distribution routing problem; multi-objective optimization; evolutionary algorithm; local search

1 Introduction

In the present era of market globalization, customization and agility, the mixed manufacturing system becomes very popular as it can meet both the diversity requirements of customers and the efficiency requirements of enterprises. In the mixed manufacturing system, the material distribution becomes an important approach to reduce the inventory cost, improve the production efficiency and enhance the benefits of enterprises [1]. In the material distribution, the vehicle routing problem (VRP) is the essence, and VRP is a hot issue in the field of operator research, applied mathematics, network analysis, graph theory, computer applications and so on during the last decade [2]. Solving VRP in mixed manufacturing systems can improve material delivery ability and decrease the inventory. It is crucial for developing an effective algorithm to optimize the material distribution routing problem (MDRP) in the mixed manufacturing system.

The MDRP is a combinatorial optimization problem which involves the optimization of the total cost (or the total travel distance) of vehicles, the number of routes and other objectives simultaneously while satisfying the given constraints, and it belongs to a type of VRPs. Due to the inherent variations in real world environment, the solutions to each VRP are often unique and satisfy an exclusive set of objectives and constraints according to the problem scenario. BODIN and BRUCE [3] proved that this type of problem was non-deterministic polynomial-time (NP)-hard. The exact algorithm can only solve small size instances [4], so many researchers have developed diversity approaches based on the heuristics and meta-heuristics methods [5-13]. These approaches sought for approximate solutions instead of exact solutions which would be at intolerably high time cost.

Various meta-heuristics methods such as generation algorithm (GA), simulated annealing (SA), neural network, Tabu search and ant colony algorithm can be found for VRPs [14-17]. For the large scale VRPs, the pure meta-heuristics have the insufficiency problems such as the intolerable time consumption or the poor quality solutions, while a hybrid approach typically outperforms pure approaches, so more and more hybrid approaches combining two or more different methods have been developed to solve the complex VRPs [5, 7, 9, 13-19]. Some successful applications strongly favor such hybrid approaches. For example, BAKER and AYECHEW [4] applied a SA to generate an initial evolution population for a GA, and their research result showed that the performance of hybrid methods was better than that of the SA and GA. TAN et al [11, 19] incorporated various heuristics for local exploration in the evolution for solving the vehicle routing problem with time windows (VRPTW). BENT and van HENTENRYCK [6] proposed a two stage hybrid algorithm that used a SA to decrease the number of routes in the first stage, while in the second stage a large neighborhood search technique was used to decrease the total travel distance. MARINAKIS and MARINAKI [7] presented a hybrid genetic-particle swarm optimization algorithm (HybGENPSO) which considered the individuals as a single swarm and a particle swarm optimizer was used to improve the physical movement of individuals following the basic principal of particle swarm optimization. As BAKER and AYECHEW [4] thought, although development of modern heuristics had resulted in considerable progress, the demands for improving performance continued.

Indeed, the single objective optimization algorithms may find the solutions with low travel cost or little number of vehicles, but it is difficult to find a solution with low travel cost and few vehicles concurrently, and not suitable for solving the complex multi-objective combinatorial optimization problems such as the MDRP. While the multi-objective evolutionary algorithms (MOEAs) can find a set of trade-off solutions for the conflicting, incompatible objectives, and it has been developed to solve the multi-objective optimization problems (MOPs). For the VRPs, the MOEAs can offer a family of Pareto optimal routing solutions covering the minimized travel cost, the number of vehicles or the other objectives. The mechanisms used in MOEA, such as the individual representation, crossover and mutation, maintaining the diversity of individuals and local search techniques, determine the performance of MOEAs. For most MOEAs, the common problem is that the used techniques are quite inefficient in terms of time consumption and solutions quality, which leads to premature convergence, inefficiency and low convergence speed [20]. Consequently, the improved MOEAs have shifted toward incorporating the heuristics and other search techniques. Additionally, the structure of the initial chromosomes population considerably impacts the performance of algorithm, and the diverse and well-structured initial population is significantly better than that generated randomly in solutions quality [21]. A MOEA and the classic heuristics are complementary, and a hybrid method is typically superior to a pure MOEA or a classic heuristics method [12]. In a hybrid method, a MOEA is used to explore the global space, whereas the heuristics are mainly employed to deal with the local search.

In this work, a hybrid multi-objective evolutionary algorithm (HMOEA) is proposed to deal with the MDRP. The overall structure of the algorithm is motivated by the recognition that the pure heuristics method may not be the most effective way to get the high quality solutions in VRPs. To solve the MDRP, the HMOEA incorporates the neighborhood search method for local exploitation and the Pareto non-dominated ranking [13] for searching the trade-off solutions. The HMOEA optimizes all objectives simultaneously without violating all the constants. It is really different from the conventional MOEAs that are characterized with simple and fixed-length code, and unchanged evolution process. The HOMEA is featured with a diverse and well-structured initial population, and the specialized genetic operators and the heuristics method are used to accommodate the sequence-oriented optimization problem in MDRP.

2 Problem formulation

2.1 Problem scenario

The MDRP model with detailed maneuver of the vehicles in a material distribution center is extended from the actual application of VRPs. It is the variance of vehicle routing problem with time window constraints (VRPTW). The additional constraints and conditions applied in MDRP show that the problems are quite complex. The MDRP model takes the rule of distribution, the material requirement plan, the movement of vehicle in the material distribution center (MDC) and working stations into account. Every material delivery route contains the material requirement plan, the location of working stations and the time and capacity constraints.

The problem modeled here is based on the daily basis requirements of MDC in the factory. MDC owns a number of vehicles to deliver the materials for the working stations, so the vehicles must leave the MDC with full loaded materials. After serving the required working stations, they return to the MDC to reload materials and deliver the required materials to the next working stations again. The task of MDC is to meet the material requirements in different working stations with lower cost and higher efficiency.

2.2 Problem modeling

1) Working stations and material delivery center

The set {W_i, i=1, 2, …, n} represents a set of working stations to be delivered, and the set {0, 1, 2, …, n} represents all the sites considered in the problem, in which 0 denotes the MDC. The travel distance between different sites (i, j) is denoted by d_{i, j}. Traveling distance satisfies the triangular inequality d_i,j+d_j,k>d_j,k. In addition, each working station has a service time s_i>0.

2) Vehicles

The set {v_k=1, 2, …, K} denotes the vehicles, and each vehicle has a capacity Q.

3) Delivery routes

The set R={0, w₁, w₂, …, w_j, 0} represents a delivery route which means a vehicle loaded with required material leaves from the MDC, visits j working stations at most once, and returns to the MDC. It is noted that W_j values are different. The travel distance of a route R, denoted by D_r, is the distance of visiting all its W_j, i.e.

       (1)

The load rate of a route R_m, denoted by L_r, is the ratio of the sum of requirements for all visited W_i to the capacity of vehicle, i.e.

                              (2)

For the nonempty set R_m: D_r>0, L_Dr>0.

4) Delivery plan

D_P represents the routing plan, a set of delivery routes {R₁, R₂, …, R_m} that every working station must be visited exactly once, i.e.

D_P={R₁, R₂, …, R_m}                          (3)

5) Time windows

The working stations and MDC have time windows. The interval [e_i, l_i] represents the time window of W_i, where e_i denotes the earliest arrival time and l_i denotes the latest leave time. As the space is limited in working stations, the materials need to arrival in the stipulated time windows. For the working station W_i, the vehicle must arrive at [e_i, l_i]:

e₀=0, l₀=0                                   (4)

                      (5)

where s_i denotes the service time of workstation i;  denotes the travel time from i^-1 to i.

After the MDRP is defined, some key objectives are considered:

Minimize the total travel distance or cost:

                           (6)

or

                          (7)

Minimize the number of routes:

f₂(R)=min(R)                                 (8)

Maximize the load rate of the vehicles:

f₃(R)=max(L_r)                                (9)

Minimize the average travel distance:

f₄(R)=min(D_r/j)                              (10)

Subject to the total material demands in a route to not exceed the vehicle capacity:

                                 (11)

The arrival time a_i satisfies the constraints of working stations:

e_iii                                    (12)

3 Hybrid multi-objective evolutionary algorithm

The proposed HMOEA is designed for solving the MDRP by minimizing both the total travel distance and number of routes simultaneously. The average travel distance and load rate of the route are used to measure the routes in the route pool.

3.1 Algorithm flowchart of HMOEA

The flowchart of the proposed algorithm is shown in Fig. 1. First, a database is built to store the information of working stations, including the distances among the working stations and material distribution central or the coordinates, the demanded materials and service time of every working station, as well as the capacity of the vehicles and service time of the material distribution center. An initial route pool is then built such that the starting working station is randomly selected, and the subsequent working stations satisfying the given constraints are selected by the heuristics method. The initial population is then built based on the route pool and heuristics method, and if the stopping criterion is not met, the individuals would be ranked by the Pareto non-dominated ranking mechanism, and the individuals are grouped into pairs randomly in order to be selected for the next specialized generation operator which includes the tournament selection, best-worst crossover and mutation. Simple heuristics are employed to improve the local exploitation at each generation, and the elite individuals measured by Pareto non-dominated ranking mechanism are preserved. The route pool is updated by keeping the best routes and discarding the worst routes measured by Eqs. (9) and (10). This evolution process iterates until no momentousness performance improvement is observed or the predefined number of generations is obtained.

Fig. 1 Flowchart of HMOEA

3.2 Chromosome representation

A variable-length chromosome representation [12, 20] is employed to build the chromosome, in which each chromosome is recorded via the path representation of the delivery routing. The chromosome is constructed by several routes which involve different working stations, and the routes are not constant but some specific sequences of working stations must be delivered. Such a representation allows the good routes or genes to be not destructed during the evolutionary operators.

3.3 Initial route pool

The initialization process begins by creating an empty route pool which is used to save or update the delivery routes found during the evolution process. At the first stage, the delivery routes are built by the heuristics method, and the starting and end points of R_r are the MDC. Then, a working station is selected as the second point of R_r. For the third and all subsequent nodes, the expanding neighborhood search method, proposed by MARINAKIS et al [22-24], is employed to select the remainder working stations. When working stations are violated by the constraints, they are deleted from the current route; if not, they are inserted to R_r. This circle of selection and insertion is continued until the route is finished. Then, the delivery route is stored in the route pool, and this process continues until the stop criterion is met. Figure 2 shows the procedure of building a feasible route. It must be made clear that there are n routes in the routes pool after the initialization, and their second nodes of these routes are the n working stations, respectively.

In order to decrease the computational burden and find out a more feasible delivery route, a small radius is used to search for a better R_i. For example, in the two-opt algorithm (see Fig. 3), if the distance of the current working station to MDC is equal to A, the initial circle has a radius of A/2. Then, the local search strategy is applied to search the next working station that satisfies the constraints for the next node of the route. If more than one candidate exist in the circle, the radius is decreased to A/4, else the radius is expanded to A, until a working station is found.

Fig. 2 Procedure of building a feasible route

Fig. 3 Expanding neighborhood search method

3.4 Initial population

Each individual in the population denotes a valid feasible D_P in which every route is feasible, denoted by D_P={R₁, R₂, …, R_m}. A heuristics method proposed by SOLOMON, called push forward insertion heuristics (PFIH), has been employed in this algorithm. In the method of PFIH, Eq. (13) denotes the first working station in each new delivery route. The detailed description of PFIH can be seen in Ref. [17]. This method is used to find the first working station for the delivery route. Once the first working station is found, the route in which the second node is the same with the found working station is selected from the route pool. If there are more than one route, the route with the highest load rate and the lowest average traveling distance is selected. Additionally, only the route in which all the working stations are un-routed is inserted to D_P:

               (13)

where α, β and γ are 0.7, 0.1 and 0.2, respectively; d_0i is the distance from the material delivery center to the working station i; b_i is the upper time; p_i is the polar coordinate angle of the working station i.

3.5 Best-worst route crossover

For the MDRP, the classical one-point crossover operator may generate infeasible solutions [3]. To overcome this insufficiency, many crossover modes have been proposed [4-9, 14-19]. The famous one is the two-point ordered crossover method proposed by ISHIBASHI et al [21]. It randomly selects two crossing points from parents, then divides the chromosomes into three parts, in which parts of the chromosomes will be inherited to the offspring. OMBUKI et al [15] introduced an improved uniform order crossover (UOX) called router crossover (RC). Further, they introduced the best cost rout crossover (BCRC) [12], which sough to find a route with the minimized number of vehicles and the lowest cost simultaneously. In BCRC, two routes R₁ and R₂ (with the lowest cost) are selected from the parents P₁ and P₂, respectively; Next, the route R₁ is inserted to P₂ to generate offspring Q₁, and R₂ is inserted to P₁ to generate Q₂. Then, the duplicated gene is removed from Q₁ and Q₂. GHOSEIRI and GHANNADPOUR [9] employed a method called the best cost-best route crossover (BCBRC), which was similar to BCRC. The best route is chosen by the average cost of each node, and the removed gene is relocated to the best possible location.

Evolution theory works on the idea that the fittest individuals survive and produce fit offspring. Generally, the crossover operation should meet the need of the exploration of unknown space and the refinement of the found region simultaneously. Most of evolutionary algorithms only consider the fittest individuals, without taking the good genes in the individuals into account, and the good genes should be inherited to the offspring. For the discussed problem, the chromosomes are constructed by a sequence of routes or genes. Motivated by the two-point crossover, BCRC and BCBRC, a specialized crossover mode called the best-worst route crossover is designed. It is constituted by the following steps: first, select the best and the worst routes from the parents P₁ and P₂ according to the criteria of load rate and average distance; then for the P₁, the worst route is replaced by the best one which is selected from the parent P₂ to generate offspring Q₁, and the same operator on P₂ to generate Q₂; third, relocate the nodes in the worst routes to the best possible location and remove the repetitive nodes in the routes. For the chromosomes which have only one route, select two segments of the route with the best and the worst average distance to exchange with the other chromosomes. This method allows the good genes in parents to be shared and inherited to the offspring, as shown in Fig. 4.

3.6 Mutation based on route pool

To complement the crossover operator, a specialized mutation mode allowing a larger search space to be explored is implemented.

TAN et al [11] introduced the partial swap, which swapped segments in chromosome randomly. The method employed in this work is similar to the partial swap. It swaps the routes between the chromosome and route pool, and the detailed processes are found in Fig. 5. First, it randomly selects a route R₁ from the chromosome to be mutated, and selects a best route R₂ in which the first working station is the same as R₁ from the route pool, then swaps these two routes; then a repair operator is used to remove the duplicated nodes and relocate the un-routed nodes to the best possible locations, and the details are as follows:

1) For the working stations appearing twice in the new swapped route and the other old routes, the new swapped route will be preserved and the working stations will be deleted from the old routes.

2) The un-routed working stations are inserted at the feasible place where the route has a lower average distance and a higher load rate, as well it must satisfy the constraints. If there is no feasible place to insert, the new chromosome will be discarded and the old parent chromosome is restored.

3) All the routes in the feasible solutions are stored in the route pool. If the route pool is full, the routes in the route pool will be ranked by the load rate and the  average distance, and the inferior routes which have lower load rate and larger average travel distance will be deleted.

Fig. 4 Illustration of best-worst routes exchange

Fig. 5 Illustration of mutation

4 Real applications involving manufacturing system material deliver

4.1 Practical case

The proposed algorithm is used to solve the practical material delivery problem in an automobile company. In this typical mixed manufacturing system, there are press shop, wedding shop, painting shop and assembly shop. The press shop and wedding shop have 42 and 50 working stations, respectively. As the required materials in the painting shop are supplied by the supplier directly, their material distribution is not considered. The assembly shop is divided into five sections: sunroof assembly line, upholstery assembly line, automotive chassis assembly line, final assembly line and inspection-test line. The inspection-test line is excluded as it is not involved in the final assembly and there are no materials for distributions. The other four assembly lines include 105 working stations which are scattered in different locations in the factory. The MDC owns a number of vehicles and its task is to meet the material requirements for all the working stations with the most cost-effective way. Currently, the scheduling of delivery routes is determined based on the experience of the dispatcher and the location of the assembly line, and the material distribution is determined based on experience of management and the clustering degree of material requirements in the working stations. For this real case, the material requirements are obtained from the MRP-II systems, for purposes of calculation. All the requirements of materials in working stations are converted into equivalent logistics, and the capacity of the vehicle is 200. The distances of working stations and MDC are measured by the real situations. The current distribution scheduling is given in Table 1 and the delivery routing is shown in Fig. 6.

It can be seen from the Table 2 and Fig. 4 that the current distribution mode is extremely ineffective. For the required services, a total of 4 762 m is traveled and 35 routes are needed. Furthermore, for the load rate, only the routes 3, 7 and 14 are fully loaded, total 17 routes are <50%, and the average load rate is only 50.4%. The low delivery efficiency is caused by the delivery method. When the working stations need the materials, the workers put down the ANDON and send the material requirements. The MDC immediately delivers the materials to the working stations once they receive the material demands.

Table 1 Route, corresponding distance and loading rate for each vehicle in case

Fig. 6 Diagram of current delivery routes

Table 2 Parameters set

This real case is solved by the proposed HMOEA. It is implemented in Matlab and tested on a computer with Intel Pentium VI 2.6G dual-processor and dual-core   64 bit under the Microsoft Windows XP operating system. As shown in Table 2, the evolutionary parameters are set to suitable values.

4.2 Computational result

For the multi-objective optimal problem discussed in this work, the solutions are a set of trade-off solutions with the minimized total travel distance and number of routes. The decision-maker can use these solutions to make a feasible material distribution routing schedule. In this practical case, before the optimization, the vehicles need to run back and forth for 35 times, and three delivery vehicles are very busy because they are expected to finish a round-trip in 40 min. After the optimization, the list of routes is given in Table 3 and the routes are simulated in Fig. 7. From Table 3, it can be seen that the total distance has a reduction of 51.6%, only two vehicles can meet the delivery requirements in which they finish a round-trip in 48 min, and the average load rate is increased to 88.25% from 50.4%.

Table 3 Comparisons between current schedule and optimal schedule

Fig. 7 Diagram of optimal delivery routes

The convergence trace of the total distance and the number of routes are plotted in Fig. 8 and Fig. 9, which show significantly the advantage of the proposed algorithm compared with the well-known multi-objective evolutionary algorithm NSGA-II [25] (the parameters set for NSGA-II are the same with HMOEA, see Table 2). It is observed that there are some distinctive disturbances in early stages of the convergence trace. These disturbances are caused mainly by the occurrence of local heuristics search. A new route is constructed based on the route pool and compared with the original one. The dominated one is preserved, and the solutions are reappraised, so the travel distance is changed with this reconstruction.

Furthermore, it can be seen from Figs. 8 and 9 that the convergence speed of HMOEA is significantly faster than NSGA-II. The main reason is that the HMOEA uses heuristics method to construct the population and BWRC to perform the crossover. On the other hand, at the generations 25-45, the convergence speed of the HMOEA is very fast and stable, because in this stage the good genes or routes are preserved after the crossover, and the routes in route pool are stable, so the mutation based on the route pool can find the better solutions.

Table 4 Optimal routes of practical case

Fig. 8 Convergence trace of HMOEA and NSGA for average travel distance

Fig. 9 Convergence trace of HMOEA and NSGA for average number of routes

5 Conclusions

1) A best-worst route crossover (BWRC) method is designed, and this method can preserve the better gene to be not destroyed in the evolution.

2) A route pool is constructed to preserve the good genes or routes, and the mutation is the routes swap based on the route pool.

3) The heuristics local search methods are used to build the population. The computational result shows that the proposed algorithm can solve the real application problem effectively. It can decrease the total travel distance by 51.66%, enhance the average vehicle load rate by 37.85%, cut down 15 routes and reduce a deliver vehicle.

References

[1] CHAKRAVORTY S S. Improving distribution operations: Implementation of material handling systems [J]. Int J Production Economics, 2009, 122: 89-106.

[2] LI Xiang-yang, TIAN Peng, LEUNG S C H. Vehicle routing problems with time windows and stochastic travel and services times: Models and algorithm [J]. Int J Production Economics, 2010, 125: 137-145.

[3] BODIN L, BRUCE G. Classification in vehicle routing and scheduling [J]. Network, 1981, 11: 97-108.

[4] BAKER B M, AYECHEW M A. A genetic algorithm for the vehicle routing problem [J]. Computers& Operations Research, 2003, 30: 787-800.

[5] TAN K C, CHEW Y H, LEE L H. A hybrid multi-objective evolutionary algorithm for solving vehicle routing problem with time windows [J]. Computational Optimization and Applications, 2006, 34: 115-151.

[6] BENT R, van HENTENRYCK P. A two stage hybrid algorithm for pickup and delivery vehicle routing problems with time windows [J]. Computers and Operations Research, 2006, 33: 875-893.

[7] MARINAKIS Y, MARINAKI M. A hybrid genetic-Particle swarm optimization algorithm for the vehicle routing problem [J]. Expert Systems with Application, 2010, 37: 1446-1455.

[8] WANG Chung-ho, LU Jiu-zhang. An effective evolutionary algorithm for the practical capacitated vehicle routing problems [J]. Journal of Intelligence Manufacturing, 2010, 21: 363-375.

[9] GHOSEIRI K, GHANNADPOUR S F. Multi-objective vehicle routing problem with time windows using goal programming and genetic algorithm [J]. Applied Soft Computing, 2010, 10: 1096-1107.

[10] CHEN J, PAN J C H, LIN C M. A hybrid genetic algorithm for the re-entrant flow-shop scheduling problem [J]. Expert System with Application, 2008, 34: 570-577.

[11] TAN K C, CHEONG C K, GOH C K. Solving multi-objective vehicle routing problem with stochastic demand via evolutionary computation [J]. European Journal of Operation Research, 2007, 177: 813-839.

[12] OMBUKI B, ROSS B, HANSHAR F. Multi-objective genetic algorithm for vehicle routing problem with time windows [J]. Applied Intelligence, 2006, 24: 17-30.

[13] MARINAKIS Y, MARINAKI M, DOUNIAS G. Honey bees mating optimization algorithm for the vehicle routing problem [C]// KRASNOGOR N, NICOSIA G, PAVONE M, PELTA D. Nature inspired cooperative strategies for optimization–NICSO 2007. Studies in Computational Intelligence, 2008, 129: 139-148.

[14] HOMBERGER J, GEHRING H. Two evolutionary meta-heuristics for the vehicle routing problem with time windows [J]. INFOR, 1999, 37, 297-318.

[15] OMBUKI B, NAKAMURA M, MAEDA O. A hybrid search based on genetic algorithms and tabu search for vehicle routing [C]// BANFF A B, LEUNG H. Proceedings of the 6th IASTED International Conference on Artificial Intelligence and Soft Computing (ASC2002). Marbella, Spain: 2002: 176-181.

[16] THANGIAH S R. A hybrid genetic algorithm simulated annealing and Tabu search heuristic for vehicle routing problems with time windows [M]// CHAMBERS L. Practical Handbook of Genetic Algorithms Complex Structures. vol.3. CRC Press: 1999, 347–381.

[17] Solomon M M. Algorithms for the vehicle routing and scheduling problems with time window constraints [J]. Operations Research, 1987, 35: 254-265.

[18] MARINAKIS Y, MIGDALAS A,PARDALOS P M. A hybrid genetic-GRASP Algorithm using Langrangean relaxation for the traveling salesman problem [J]. Journal of Combinatorial Optimization, 2005, 10: 311–326.

[19] TAN K C, LEE T H, CHEW Y H., LEE L H. A multi-objective evolutionary algorithm for solving vehicle routing problem with time windows [C]// Proceedings of the IEEE International Conference on Systems, Man and Cybernetics, Washington, DC, USA, 2003: 361-366.

[20] COELLO C C A. Evolutionary multi-objective optimization: A historical view of the field [J]. IEEE Computational Intelligence Magazine, 2006, 1: 28-36.

[21] ISHIBASHI H, AGUIRRE H, TANAKA K, SUGIMURA T. Multi-objective optimization with improved genetic algorithm [C]// Proceedings of the IEEE International Conference on Systems, Man, and Cybernetics (SMC). Nashville, 2000: 3852–3857.

[22] MARINAKIS Y, MIGDALAS A, PARDALOS P M. Expanding neighborhood GRASP for the traveling salesman problem [J]. Computational Optimization and Applications, 2005, 32: 231–257.

[23] MARINAKIS Y, MIGDALAS A, PARDALOS P M. Multiple phase neighborhood search GRASP based on Lagrangean relaxation and random back tracking Lin Kernighan for the traveling salesman problem [J]. Journal of Combinatorial Optimization, 2007, 18: 201–216.

[24] MARINAKIS Y, MIGDALAS A, PARDALOS P M. A new bi-level formulation for the vehicle routing problem and a solution method using a genetic algorithm [J]. Journal of Global Optimization, 2007, 38: 555-580.

[25] DEB K, PRATAP A, AGARWAL S, MEYARIVAN T. A fast and elitist multiobjective Genetic algorithm: NSGA-II [J]. IEEE Transactions on Evolutionary Computation, 2002, 6(2): 182-197.

(Edited by YANG Bing)

Foundation item: Project(50775089) supported by the National Natural Science Foundation of China; Project(2007AA04Z190, 2009AA043301) supported by the National High Technology Research and Development Program of China; Project(2005CB724100) supported by the National Basic Research Program of China

Received date: 2010-12-23; Accepted date: 2011-03-27

Corresponding author: GAO Gui-bing; Tel: +86-731-58290040; E-mail: gaoguibing@163.com