- Email: [email protected]
Hybrid multiobjective evolutionary algorithm with fast sampling strategy-based global search and route sequence difference-based local search for VRPTW
Hybrid Multiobjective Evolutionary Algorithm with Fast Sampling Strategy-based Global Search and Route Sequence Difference-based Local Search for VRPTW
Journal Pre-proof
Hybrid Multiobjective Evolutionary Algorithm with Fast Sampling Strategy-based Global Search and Route Sequence Difference-based Local Search for VRPTW Wenqiang Zhang, Diji Yang, Guohui Zhang, Mitsuo Gen PII: DOI: Reference:
S0957-4174(19)30868-1 https://doi.org/10.1016/j.eswa.2019.113151 ESWA 113151
To appear in:
Expert Systems With Applications
Received date: Revised date: Accepted date:
18 August 2019 3 December 2019 19 December 2019
Please cite this article as: Wenqiang Zhang, Diji Yang, Guohui Zhang, Mitsuo Gen, Hybrid Multiobjective Evolutionary Algorithm with Fast Sampling Strategy-based Global Search and Route Sequence Difference-based Local Search for VRPTW, Expert Systems With Applications (2019), doi: https://doi.org/10.1016/j.eswa.2019.113151
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.
Highlights • HMOEA-GL combines FSS-GS and RSD-LS. • FSS-GS quickly explores the entire solution space. • RSD-LS further enhances the search ability of HMOEA-GL. • Designing suitable coding method and proper genetic operators. • Simple insertion search is used to reduce the number of vehicles.
1
Hybrid Multiobjective Evolutionary Algorithm with Fast Sampling Strategy-based Global Search and Route Sequence Difference-based Local Search for VRPTW Wenqiang Zhanga,∗, Diji Yanga , Guohui Zhangb , Mitsuo Genc a College
of Information Science and Engineering, Henan University of Technology, China of Management Engineering, Zhengzhou University of Aeronautics, China c Fuzzy Logic Systems Institute/Tokyo University of Science, Japan
b School
Abstract The vehicle routing problem with time windows (VRPTW) is an important and widely studied combinatorial optimization problems. This paper aims at VRPTW with the objectives of reducing the number of vehicles and minimizing the time-wasting during the delivery process caused by early arrival. To solve this NP-hard problem, a hybrid multiobjective evolutionary algorithm with fast sampling strategy-based global search and route sequence difference-based local search (HMOEA-GL) is proposed. Firstly, fast sampling strategy-based global search (FSS-GS) of HMOEA-GL extensively explores the entire solution space to quickly guide the search direction towards the center and edge areas of Pareto frontier. Secondly, route sequence difference-based local search (RSD-LS) is executed on the individuals with poor performance in the population obtained by FSS-GS to enhance the search ability of HMOEA-GL. In addition, the suitable coding method and proper genetic operators are designed, especially, a simple insertion search is used to reduce the number of vehicles in VRPTW. Comparing with NSGA-II, SPEA2, and MOEA/D, experimental results on 12 Solomon benchmark test problems indicate that the proposed HMOEA-GL is effective, and more excellent in convergence, while maintaining a satisfying distribution ∗ Corresponding
author Email addresses: [email protected] (Wenqiang Zhang), [email protected] (Diji Yang), [email protected] (Guohui Zhang), [email protected] (Mitsuo Gen)
Preprint submitted to Expert Systems with Applications
December 20, 2019
performance. HMOEA-GL could be an effective intelligent algorithm for expert and intelligent decision support system to help logistics companies to make decisions. Keywords: multiobjective evolutionary algorithm, vehicle routing problem, global search, local search, transportation
1. Introduction Under the environment of economic globalization and information technology, people’s online shopping and long-distance ordering is more and more convenient and popular, and logistics plays an increasingly obvious role in mod5
ern economic activities, at the same time, the logistics problem becomes more complex. It is very necessary to design and use efficient intelligent decision support system(IDSS) to help decision-makers deal with logistics problem. IDSS is a decision support system that makes extensive use of artificial intelligence (AI) techniques, which works like a human consultant to support the decision-
10
makers. In-depth analysis of the logistics problem, it can be found that the core is the vehicle routing problem (VRP). A good solution of VRP can take into account many aspects, not only can complete the customer service on time and on demand, but also can reduce the cost of distribution center consumption, to achieve a win-win. With the progress of the Internet, the type of customers
15
of VRP is more likely to be transformed from large-volume customers like production base and supermarkets to small-volume customers such as family and individuals. Unlike large-volume customers who can receive service 24 hours a day, small-volume customers are generally only allowed to receive service within a few hours. Therefore, the consideration of time window is very necessary. At
20
the same time, small-volume customers generally do not have a parking space, so the waiting cost of vehicles in the city may increase significantly, should not be ignored. With the efforts of researchers, great progress has been made in solving this problem, but due to the complexity of the problem and the differ-
3
ent research background, the research on this problem is still very important 25
and ongoing. VRP is one of the most important and widely studied combinatorial optimization problems, which is relevant to transportation logistics such as post, parcel and distribution services. Its objective is to obtain the lowest-cost set of routes. In the beginning, truck dispatching problem was proposed by Dantzig
30
and Ramser Dantzig & Ramser (1959) as a generalization of the travelling salesman problem (TSP), cost has mostly been associated with the travel distance. In 1964, Clarke and Wright summarized this problem as vehicle routing problem that is well-known: i.e., how to use a fleet of trucks to serve a set of customers which have different locations and needs. The reason VRP can be one of the
35
most widely studied topics in the field of operations research is that VRP have important practical significance and with high complexity when this problem integrated into real-life. However, the current VRP model is very different from the original, because when VRP is combined with real-life, there are many factors that cannot be ignored. Like the number of vehicles, delivery time, parking
40
cost, makespan, and workload balance, etc. are important to consider. Considering these realistic factors, VRP have several common variants that involve different constraints, like time-dependent travel times, with time windows, support delivery while pickup, can update customer information in real time, etc. Among them, the variant with time windows (VRPTW) has vehicles with lim-
45
ited capacity and the specific delivery time windows, and is particularly relevant to practical applications. The main problem studied in this paper is the vehicle routing problem with time windows which objective to minimize the number of vehicles and the waste time caused by the premature arrival. Because the limitations of many constraints of VRP, so it is often to con-
50
sider the VRP as a multiobjective problem. VRP is also considered as an NP-hard problem Lenstra & Kan (1981), when solving large-scale real-world VRPs, heuristics and meta-heuristics are often more suitable, exact algorithms are only efficient for small problem. The same is true for the variant VRPTW, optimal solutions for small instances of the VRPTW can be obtained by exact 4
55
methods, but the computation time required increases considerably for larger instances Desrochers et al. (1992). For this reason, in recent years, various heuristic and meta-heuristics algorithm for solving VRPTW have become the hotspot of research. Algorithms such as simulated annealing (SA), tabu search (TS), ant colony optimization (ACO), genetic algorithm (GA), particle swarm
60
optimization (PSO), differential evolution (DE), have proven to be effective in solving complex multiobjective problems. These algorithms are also used to solve complex VRPTW. Heuristic algorithms can greatly reduce the computational complexity and have many successful applications. However, with the increasing of the complexity and scale of the VRPTW, there may suffer some
65
problems. For example, the speed of the algorithm or the final solution set is still not satisfactory. Therefore, research on VRPTW is still ongoing. According to the No Free Lunch Theorem Wolpert et al. (1997), no algorithm can solve all the problems. Each of the classic algorithms has its own unique advantages, and inevitably there are some drawbacks. Proper combination of
70
different algorithms with different characteristics may further enhance the search effectiveness by adopting the advantages of each algorithm, and consequently may overcome the inherent limitations of single algorithm Tang & Wang (2013). Therefore, it is naturally to mix different techniques and utilize their advantages for solving complicated problems Zhou et al. (2011).
75
As we know, GA is one of the most popular type of evolutionary algorithms (EA), which can easily be adapted to various types of problems Holland (1992)Gen & Cheng (2000)Gen et al. (2008)Yu & Gen (2010)Tasan & Gen (2012). DE is a direction-based stochastic search technique which can improve the convergence of the algorithm by lead the poor solution to close to the good
80
solution Storn & Price (1997). But both of these algorithms have some drawbacks, which is easy to converge to local optima or even arbitrary points rather than the global optimum, resulting in premature convergence. This challenge is particularly hard when the dimension of problems is high and there are a lot of local optima.
85
Therefore, it is an interesting but challenging to hybridize GA and DE re5
garding the evolutionary behavior and potential advantages for VRPTW. Inspired by the successful hybrid strategies and application in previous research, in this paper, hybrid multiobjective evolutionary algorithm with fast sampling strategy-based global search and route sequence difference-based local search 90
(HMOEA-GL) is proposed to solve complicated multiobjective VRPTW. The objectives of the problems one is to reduce the number of vehicles and another is to minimize the time-wasting during the delivery process caused by early arrival. HMOEA-GL combines fast sampling strategy-based global search (FSS-GS) and route sequence difference-based local search (RSD-LS). FSS-GS is
95
designed as the global search strategy which includes center area-based elite sampling strategy and edge area-based sampling strategy. Center area-based elite sampling strategy according to Pareto dominating and dominated relationshipbased fitness function (PDDR-FF) can fast picking out the individuals located in the central area of the Pareto frontier Zhang et al. (2014). Edge area-based
100
sampling strategy based on vector evaluated genetic algorithm (VEGA) which prefers the individuals located in the edge area of the Pareto frontier Schaffer (1985). By mixing these two strategies and executing suitable genetic operators as well as simple insertion search (SIS), FSS-GS can rapidly improve the convergence and distribution performances toward the center and edge regions
105
of Pareto frontier. As a follow-up step, drawing on the evolutionary strategy of DE, RSD-LS is used to further improve the quality of individuals. RSD-LS guides individuals with poor performance to move toward better performing individuals. By hybridizing FSS-GS and RSD-LS, HMOEA-GL can achieve fast convergence and sufficient distribution. As the core of IDSS, HMOEA-GL
110
could be used to solve various scheduling problems such as VRPTW, which can provide multiple appropriate solutions for decision-makers by different weights of the objectives according to decision-makers’ needs. With HMOEA-GL, an IDSS can be established to effectively support decision-makers in a series of scheduling decision-making processes such as vehicle scheduling.
115
The rest of this paper is organized as follows. Section 2 introduces the related work to solve the vehicle routing problem with time windows. Section 6
3 provides definitions of VRPTW. A detailed description of the HMOEA-GL and the main strategies are shown in Section 4. Experimental design and the analysis of results are detailed in Section 5. Conclusions and future work are 120
provided in Section 6.
2. Related Works Over the past decades, the VRP and its variants have grown ever more popular in the academic literature. Braekers et al. present a taxonomic review of the VRP literature published between 2009 and June 2015 Braekers et al. 125
(2016). For VRPTW, Dixit et al. reviewed some of the recent advancements in the VRPTW using meta-heuristic techniques Dixit et al. (2019). SA inspired from annealing in metallurgy is often combined with other algorithms to solve the VRPTW. Pratiwi et al. proposed a natural inspired algorithms for dealing with constraints of VRPTW which involves bat algorithm and
130
cat swarm optimization. Among the proposed algorithms, bat algorithm is hybridized with SA, the worst solution of bat algorithm is replaced by the solution from SA. Algorithm which is based on behavior of cats, cat swarm optimization, is improved using crow search algorithm to make simpler and faster convergence Pratiwi et al. (2018). Gonz´ alez et al. proposed a memetic algorithm (MA) with
135
SA as trajectory-based method for solving the capacitated vehicle routing problem with time windows (CVRPTW). In MA, a novel crossover operator called the single breaking-point sequence based crossover (SBSBX) is introduced. One of the principles behind the design of SBSBX is to reduce the disruptive behavior of SBX, with the aim of providing additional intensification Gonz´ alez et al.
140
(2017). Baˇ nos, Ra´ ul et al. proposed a multiobjective procedure based on SA and the multiple temperature Pareto SA to solve VRPTW of considering the travelled distance and the imbalance of the routes Baˇ nos et al. (2013). TS is a meta-heuristic search method employing local search methods used for mathematical optimization, it was used to solve VRPTW long ago. Schnei-
145
der proposed a tabu search algorithm for the VRP with time windows and
7
driver-specific times, a variant of the classical VRPTW that uses driver-specific travel and service times to model the familiarity of the different drivers with the customers to visit Schneider (2016). Wang et al. focused on the heterogeneous multi-type fleet vehicle routing problem with time windows and an incompatible 150
loading constraint, developed a mathematical model, proposed a ruin-recreate heuristic algorithm, and a threshold tabu search method Wang et al. (2015). Yu et al. proposed a hybrid approach, which consists of TS and ACO to solve VRPTW. In proposed algorithm, a neighborhood search is introduced to improve the performance of ACO, TS is used to maintain the distribution and
155
explore new solutions when close to the convergence Yu et al. (2011). Liu and Lee proposed a mathematical model of inventory routing problem with time windows to solve the VRPTW and the inventory control decision problem simultaneously. The proposed method is divided into two phases, first is to find the initial solution and the second phase is to improve the solution adopting the
160
variable neighborhood TS Liu & Chen (2011). ACO is a probabilistic technique for solving computational problems which can be used to find good paths through graphs, has a good performance in computer science and operations research, it is also widely used in VRPTW. To solve the vehicle routing problem with service time customization (VRPTW-STC),
165
which the aim is to find an optimum solution to the smallest fleet size, the lowest travelling distance as well as the largest total service time of all customers, Wang et al. designed a multi ant system based hybrid heuristic algorithm which inspired by decomposing a multiobjective problem into several single objective ones. ACO algorithms are applied to every single-objective problem. Several
170
local search algorithms are also incorporated into MAS to improve the solution quality Wang & Xing (2018). Gupta and Saini proposed an improved ant colony optimization to solve VRPTW, which uses a new pheromone reset and update function to enhance the route searching and a 2-opt method to improve the optimized path. The experimental results verify the effectiveness of the algorithm
175
and confirms it as the better path optimizing method for VRPTW, in comparison with the traditional ACO and other heuristics Gupta & Saini (2017). Zhang 8
et al. proposed a solution strategy based on ant colony optimization and three mutation operators for multiobjective VRP with flexible time windows. In this special problem, a fleet of vehicles can service a set of customers earlier and later 180
than the required time with a given tolerance, in the hope of enabling a logistics company to save distribution costs at the expense of customer satisfaction Zhang et al. (2019). GA based on recombination, mutation and natural selection also have many related researches, which applied GA on VRPTW. G¨ oc¸ken et al. proposed a
185
hybridized version of genetic algorithm which is a meta-heuristic solution technique with constructive heuristic methods for solving VRPTW. By using sweep algorithm in initial population generation phase of GA, the search starts with high quality solution sets. Comparing with the results of GA based on nearest neighbor, the proposed GA reaches more effective solutions by beginning with
190
sweep based initial population G¨ oc¸ken et al. (2017). Abidi et al. focused on a class of problems known as rich vehicle routing problem (RVRP) with dynamically changing orders and time windows constraints. The proposed GA with a simple heuristic to solve the dynamic vehicle routing problem with time windows Abidi et al. (2018). Shi et al. proposed a hybrid GA integrated with stochas-
195
tic simulation methods to solve VRPTW with uncertain demand. Among, the uncertain demand is considered as a fuzzy variable, a fuzzy chance constraint model is designed Yong et al. (2016). Moon et al. studied VRPTW and extend two cases that allows overtime for drivers and the possibility of using outsourced vehicles. A mixed integer programming model, a hybrid algorithm based on SA
200
and GA is developed to solve this problem Moon et al. (2012). Long et al. proposed a Pareto-based evolutionary algorithm for solving prize-collecting vehicle routing problem (PCVRP). In PCVRP, besides the vehicle assignment and visiting sequencing, customer selection also must be determined. The proposed algorithm combined GA and a local search strategy, a decomposition strategy
205
is designed to decompose PCVRP into multiple sub problem Long et al. (2019). PSO is a population based stochastic optimization technique which inspired by social behavior of bird flocking or fish schooling, can also be used to solve 9
VRPTW. Marinakis et al. proposed a new variant of PSO for the solution of the VRPTW. Three different adaptive strategies are used in the proposed 210
multi-adaptive particle swarm optimization (MAPSO) algorithm. First, adaptive strategy concerns the use of a greedy randomized adaptive search procedure (GRASP). Second, new adaptive strategy concerns the adaptiveness in the movement of the particles from one solution to another. Third, the adaptiveness is in all the parameters of the particle swarm optimization Marinakis et al.
215
(2019). Zhang et al. proposed an evolutionary scatter search particle swarm optimization algorithm (ESS-PSO) to solve the VRPTW. In ESS, a genetic algorithm and a new ’route+/-’ evolutionary operator are introduced in scatter search template. Also Zhang et al. proposed a discrete PSO that sets the route-segment as the velocity of particles and in which the velocity and position
220
updating rules are designed based on the concept of ’ruin and recreate’. A new solution representation called ’auxiliary code’ based on customer allocation is designed to maintain the distribution of the reference set Zhang et al. (2018). Jiang et al. combined fields of logistics and unmanned aerial vehicle (UAV) created a model of task assignment for UAV in logistics. The model takes into
225
account multi-constraints (such as weight coefficients, time windows constraints, the constraints of the UAV and soon), and proposed an improved PSO to solve the task assignment problem with multiple constraints Jiang et al. (2017). Based on the differences between individuals, the DE can guide individuals with poor fitness to close to excellent one, which also has many research
230
works in solving the VRPTW. Pu et al. proposed a hybrid differential evolution (HDE) optimization for the vehicle routing problem with time windows and driver specific times (VRPTW-DST). In this problem, adopts service times and driver-specific travel times to model the degree of familiarity to different drivers with the service for customers. By comparing the results of experiment
235
and simulation, the knowledge for drivers in the routing problem may effectively enhance the efficiency for vehicle routesPu et al. (2017). Salamatbakhsh-Varjovi et al. developed a bi-objective mathematical model to evaluate a periodic vehicle routing problem with time windows under uncertainty for companies, and 10
an improved differential evolution algorithm was used to identify efficient solu240
tionsSalamatbakhsh-Varjovi et al. (2018). After many years of efforts by many researchers, there are already many methods to solve VRP, but still face some problems, such as the speed of the algorithm or the quality of the final solution is not satisfactory. Moreover, different researchers concern the VRP under different environmental backgrounds.
245
The influence of different environments leads to different conditional constraints, so the proposed method may only be applied to certain specific situations. In addition, the combination of different technologies to form a hybrid algorithm has become a new trend. Hybrid algorithm may have unexpected effects for solving VRP, and there are many combination methods of hybrid algorithm
250
still have not been attempted yet. Therefore, research on VRP issues is still very important and is still ongoing.
3. Problem Definitions The VRPTW is one of a variants of VRP which is NP-hard. VRPTW can be seen as a combination of the TSP and the bin packing problem (BPP). In 255
this paper, the objectives of the problem are to find set of routes with minimumwaste time and less number of vehicle to service customers. The customers have specific delivery time windows, and the vehicles are identical and have limited capacity Garcia-Najera & Bullinaria (2011). An instance of VRPTW can be defined as follows. There is a set V =
260
{0, . . . , N } of vertices, which the vertex 0 is the depot, and the other vertices
V 0 = V \{0} are defined as customers. Each customer i ∈ V 0 located at coordinates (xi , yi ), has a demand of goods gi > 0, a delivery time window [bi , ei ], bi is the beginning time of customer allow receive goods, and ei is the ending time of time windows for allow receive goods. Each customer also has a service 265
time si for unloading goods. The depot is positioned at (x0 , y0 ), has demand g0 = 0, and time window [0, e0 ≥ max{ei : i ∈ V 0 }]. The vehicles are the same and with the capacity Q ≥ max{gi : i ∈ V 0 }. 11
The distance dij and the passing time tij between the vertices i and j are generally considered to be the Euclidean distance. The problem is to find a set 270
of k routes R = {r1 , . . . , rk }, each route begins and ends at the depot, and each customer is serviced once. So each vehicle is assigned to a set of customers, and the vehicle capacity Q must be able to satisfy the sum of these customers’ demands. If rk = (u(1, k), . . . , u(Nk , k)) is the sequence of Nk customers supplied in the k-th route, where u(i, k) is the i-th customers to be visited in route k,
275
then Vk0 = u(1, k), . . . , u(Nk , k) is the set of customers serviced. There is no need to display the depot, but when calculating the fitness function, it is necessary to add the depot at the beginning and the end of the route. Dk is total distance of routing rk can be expressed as:
Dk =
Nk X
d(u(i,k),u(i+1,k))
(1)
i=0
Suppose a(u(i, k)) is the arrival time of the vehicle k at vertex i and l(u(i, k)) 280
means the time to leave the vertex. The arrival time can be calculated as follows.
a (u (i, k)) = l (u (i − 1, k)) + t(u(i−1,k),u(i,k))
(2)
Late service is not allowed, that means the route is invalid if the arrival time is later than the end of the customer’s time window. However, arriving early is allowed, but the vehicle has to wait until the beginning of the customer time window, so the vehicle has a waiting time. 0, ifa (u (i, k)) ≥ b u(i,k) w (u (i, k)) = b − a (u (i, k)) , otherwise
(3)
u(i,k)
285
The leaving time from the i-th customer in route k is described as:
l(u(i, k)) = a(u(i, k)) + w(u(i, k)) + su(i,k)
(4)
The total time of route rk can be expressed as:
Tk =
Nk X
(t(u(i,k),u(i+1,k)) + w u (i + 1, k) + su(i+1,k)
i=0
12
(5)
This paper will concentrate on two objectives, one is to reduce the number of vehicles and another is to minimize the time-consuming waste during the delivery process caused by the premature arrival. The first objective number of 290
vehicles can be expressed as:
f1 (R) = |R| = K
(6)
The second objective is total waiting time in each route k. It is necessary to mention that the waiting time of the first vertex which be serviced of each vehicle is set to zero.
f2 (R) =
Nk K X X
w (u (i + 1, k))
(7)
k=1 i=0
The relationship between customer’s demand and vehicle capacity in route 295
rk can be expressed as:
Gk =
X
0 i∈vk
gi ≤ Q, ∀k = 1, . . . , K
(8)
The arrival times need to meet customer’s time window constraints.
a (u (i, k)) ≤ eu(i,k) , ∀k = 1, . . . , K, 1 ≤ i ≤ Nk
(9)
4. Algorithm Description As the VRPTW belongs to NP-hard problem with the high computational complexity and lots of conditional constraints, in this paper HMOEA-GL is pro300
posed for solving VRPTW. In this section, section 4.1 introduces the overview of HMOEA-GL. Section 4.2 describes the details of FSS-GS. Section 4.3 shows the genetic representation of HMOEA-GL, including encoding and decoding. The genetic operators of HMOEA-GL will be described in section 4.4. Section 4.5 introduces the process of SIS and RSD-LS will also be introduced in section
305
4.6.
13
4.1. Overview of the HMOEA-GL The main advantages of HMOEA-GL are hybridizing FSS-GS and RSD-LS, FSS-GS has fast convergence and sufficient distribution, and subsequently using RSD-LS to further improve the performance of the algorithm. Among FSS-GS, 310
center area-based elite sampling strategy according to a fitness function called PDDR-FF Zhang et al. (2014). PDDR-FF is used to evaluate the individuals that can fast determine which individuals are nondominated or dominated but have large domination area. Combining with edge area-based sampling strategy based on VEGA Schaffer (1985), FSS-GS increases the attention to individuals
315
which located on the edge and center regions of the Pareto front. Therefore, FSS-GS can effectively improve the convergence performance of HMOEA-GL while ensuring an acceptable distribution performance. FSS-GS is designed as the global search strategy of HMOEA-GL. In addition, RSD-LS is a local search which can guide the individuals with poor performance move closer to
320
those who perform better. By this way, the search ability of HMOEA-GL can be further enhanced. In addition, in order to reduce the number of vehicles, a single objective optimization method SIS is used in HMOEA-GL, and other suitable genetic operators are also used to recombine chromosomes. HMOEAGL can be divided into four steps roughly, and the framework of HMOEA-GL
325
is shown as follows, and the framework of HMOEA-GL is shown in figure 1, and the pseudo code is shown in algorithm 1. Step 1: Using the edge area-based sampling strategy, individuals are selected from the population Pt by binary tournament selection. Individuals who perform well on objective 1 will be put into SP1 . Similarly, the individuals
330
which perform well on objective 2 will be placed in SP2 . Step 2: Mixing SP1 , SP2 and elite population At as a mating pool, and pick out the parent individuals from it. Performing crossover operator, and mutation operator on parent individuals to generate the offspring. Executing 0 . SIS on offsprings to form temporary population Pt+1
335
Step 3: Performing RSD-LS operations on individuals in the temporary 0 population Pt+1 . Based on the difference between the two individuals, the
14
Start
Initializing Population Pt and Elite Population At
Yes
Stop
Output Pareto solution set
No FSS-GS: edge area-based sampling strategy
Step1
Creating Sub-Population SP1 and SP2 according to VEGA 1. Selecting individuals in Pt who perform well on objective 1 into SP1 2. Selecting individuals in Pt who perform well on objective 2 into SP2
Mixing SP1, SP2 and At as Mating pool Genetic operators: Selection Crossover Mutation
Step2
Simple Insertion Search
Creating temp. Population P’ t+1 1. Creating offSpringSet by perform genetic operators on Mating pool 2. Performing simple insertion search on offSpringSet
Route Sequence Difference-based Local Search
FSS-GS: center area-based elite sampling strategy
Step3
Creating Population Pt+1 1. Creating RSD-LS Population RSD-LSP’ t+1 by performing RSD-LS on P’ t+1 2. Selecting individual from RSD-LSP’ t+1 and P’ t+1 according to PDDR-FF
Step4
Updating Elite Population At to Elite Population At+1 1. Uniting Pt+1 and At as tempEliteSolutionset 2. Selecting individual from tempEliteSolutionset according to PDDR-FF
FSS-GS: center area-based elite sampling strategy
Figure 1: Basic framework of the HMOEA-GL.
15
Algorithm 1 HMOEA-GL Input: data set of VRPTW (for each instance); EA parameters: PopSize, MaxGenerationsSize, ElitePopSize, subPopSize; Output: Pareto SolutionSet; 1: t = 0; 2: Initializing population Pt by encoding routine; 3: Calculating objectives f1 , f2 by decoding routine; 4: Calculating fitness eval by PDDR-FF routine; 5: Initializing elite population At by center area-based elite sampling strategy; 6: while (t ≤ MaxGenerationsSize) do 7:
Creating subPopulation1 SP1 by edge area-based sampling strategy;
8:
Creating subPopulation2 SP2 by edge area-based sampling strategy;
9:
Mixing SP1 , SP2 , and At as Mating pool;
10:
Creating offspringSet by genetic operators on Mating pool;
11:
Performing simple insertion search on offspringSet;
12:
0 Temp population Pt+1 = offspringSet;
13:
0 0 ; Performing RSD-LS on Pt+1 to create RSD-LS Pt+1
14:
0 0 and RSD-LS Pt+1 Creating temp. solutionSet tempSet union of Pt+1 ;
15:
Creating Pt+1 from tempSet by center area-based elite sampling strategy;
16:
Calculating objectives f1 , f2 by decoding routine;
17:
Calculating fitness eval by PDDR-FF routine;
18:
Creating a temp. elite solutionSet tempEliteSet union of Pt+1 and At ;
19:
Creating At+1 from tempEliteSet by center area-based elite sampling strategy;
20:
t = t + 1;
21: end while 22: Creating Pareto SolutionSet by selecting nondominated in At ;
16
individual with poor fitness function can obtain a guiding evolution. PDDR-FF will be used to evaluate the new individual generated by RSD-LS, which judges the new individual whether can be accepted or not. The population optimized 340
by RSD-LS is determined as the new population Pt+1 . Step 4: Using the center area-based elite sampling strategy to update the elite population. The individuals both in the new population Pt+1 and elite population At are evaluated. Best | A | individuals will be selected and formed the new elite population At+1 .
345
HMOEA-GL can be divided into two parts. One is to use the FSS-GS to fast improve the convergence and distribution performance of the HMOEA-GL, and another is to use RSD-LS to guide individuals with poor fitness to further evolve, thus the performance of HMOEA-GL can be improved. The specific details of FSS-GS, RSD-LS, SIS, coding method and various genetic operators
350
in HMOEA-GL will be described later. 4.2. Fast Sampling Strategy-based Global Search FSS-GS can fast converge toward the center and edge regions of Pareto frontier while ensuring sufficient distribution by combining the center area-based elite sampling strategy and the edge area-based sampling strategy. Among,
355
a fitness function called PDDR-FF plays a key role in center area-based elite sampling strategy Zhang et al. (2014). PDDR-FF is used to evaluate individuals according to the dominant relationship among individuals, and each individual can obtain an evaluation value. With the evaluation value of PDDR-FF, FSSGS can fast identify which individuals are nondominated. Moreover, individuals
360
who are dominated but have a larger domination area (have great potential to become nondominated) are also evaluated as good. The specific formula of PDDR-FF for an example of calculating the fitness function of individual Si will be shown in following.
eval (Si ) = q (Si ) +
1 , i = 1, 2, . . . , popSize p (Si ) + 1
17
(10)
where q(Si ) is the number of the individuals that dominate Si , p(Si ) is the 365
number of the individuals who are dominated by Si , and popSize is the size of population. Obviously, the smaller the eval, means the individual is better. The ability of converging to the center area of Pareto frontier can be effectively improved by the center area-based elite sampling strategy. At the same time, in order to maintain the distribution performance of the algorithm,
370
FSS-GS combines center area-based elite sampling strategy with the edge areabased sampling strategy, which is according to the sampling strategy of VEGA Schaffer (1985). The sampling strategy of VEGA can be simply described as decomposing multiobjective problem into multiple single-objective problems, then set a sub-population for each objective. Each sub-population contains individ-
375
uals who prefer one certain objective. The edge area-based sampling strategy prefers to the individuals located in the edge area of Pareto frontier, that can make up for the deficiency of distribution performance of center area-based elite sampling strategy effectively. Taking a bi-objective minimization problem as an example, the main mechanism of FSS-GS is shown in figure 2. min f2
Edge area-based sampling strategy
Individuals prefer objective 1
Center area-based elite sampling strategy
Individuals with a good eval Individuals prefer objective 2 min f1
Figure 2: The main mechanism of Fast Sampling Strategy-based Global Search.
380
According to the center area-based elite sampling strategy, FSS-GS can fast pick out the individuals located around the central area on the Pareto frontier,
18
and the edge area-based sampling strategy prefers the edge area on the Pareto frontier. Mixing two sampling strategies, FSS-GS can effectively improve the convergence and distribution performance of the HMOEA-GL. 385
4.3. Genetic Representation 4.3.1. Encoding and Initialization The encoding of HMOEA-GL is simple and intuitive, which uses permutation code. The length of the chromosome is equal to the number of the customer vertices. Each gene of chromosome represents the one of the customer number.
390
In one chromosome, one customer number will only appear once and recurring is not allowed. In order to easily apply genetic operators on chromosome, there is no depot in the chromosome. But the depot needs to be taken into account when decoding and calculating fitness. When initializing, the population is formed of chromosomes which are randomly generated.
395
4.3.2. Decoding When decoding, the genes of chromosome will be traversed by the order from the first one to the last. Cutting the chromosome when exceeded the conditional limit, then the chromosome is separated into multiple sub-routes. One vehicle is designated to serve customers according to one of the sub-route. It should
400
be noted that each sub-route needs to add the depot at the beginning and the end. An example of decoding is shown in figure 3, and pseudo code of decoding in HMOEA-GL is also be shown in algorithm 2. Taking figure 3 as an example, it can suppose that there is a depot 0, and a group of customers numbers is 1-6. Each customer has a time window indicated
405
by square brackets and the distance between the customers determines the travel time of the vehicle. The service time for each customer is set to 1. For customer 1, the begin time of receiving goods is 1, and the end time is 6. The distance between the depot 0 and customer 1 is set to 1, then distance between the depot 0 and customer 6 is approximately considered to be 1.3. A chromosome {1, 5,
410
4, 2, 3, 6} will be decoded into two sub-routes {0, 1, 5, 4, 2, 0} and {0, 3, 6, 19
[5, 10] [3, 7]
4
1
0.7 [9, 15]
5 Chromosome:
1
5
4
2
3
6
2
1 [2, 6]
Vertex sequence:
1
5
4
2
3
1
6
1
0
1
exceeded conditional limit
Sub-route:
0
1
5
4
Vehicle_1
2
0
0
3
6
0
1
3
1.3
Vehicle_2
[2, 7] 1.5
6 [2, 5]
Figure 3: Decoding method in HMOEA-GL.
0}. When the vehicle finishes serving the customer 2, the end time is 10, it has no enough time to complete the service requirements of customer 3, so the vehicle is returned to the depot and use another vehicle to serve the customer 3. Conditional limit is not only limited by time, but also the load of the vehicle. 415
4.4. Genetic Operators 4.4.1. Selection Operator In HMOEA-GL, the binary tournament selection is used as the selection operator. Since the purpose of using the edge area-based sampling strategy based on VEGA is to improve the distribution of the algorithm, so in the process
420
of forming sub-populations, it is not necessary to divide all the individuals which have excellent performance on one particular objective into one sub-population. Instead, select two individuals randomly and choose the one who performs well on particular objective, and put it into the corresponding sub-population. The center area-based elite sampling strategy according to PDDR-FF is used to
425
improve the convergence of the algorithm. So is often to order the individuals of the population by their fitness values, then picks out the best part of them as the elite population.
20
Algorithm 2 Decoding procedure of Chromosome Input: VerticesData x, y, b, e, s, g; Distance d(i, j)=Euclidean(vi , vj ); VehicleData velocity V =1, capacity Q; Chromosome Sequence C; Output: Routes R = r1 , . . . , rk ; 1: Vn =0; // Vn is the vehicle number 2: cg =0; // cg is the current load of vehicle 3: at = bC(0) ; // at is the arrival time of vertex t 4: lt = bC(0) + sC(0) ; //lt is the leave time of vertex t 5: bdt ; // bdt is the time from vertex t back to depot 0 6: add v0 into r0 ; // add the depot as the beginning of route 7: add C(0) into r0 ; // add the first vertex in C into first route r0 8: cg = cg + gC(0) ; 9: for i = 1 to length of C do 10:
at = lt + d(C(i − 1), C(i))/V ;
11:
cg = cg + gC(i) ;
12:
bd1t = bC(i) + sC(i) + d(C(i), 0)/V ; //for early arrival
13:
bd2t = at + sC(i) + d(C(i), 0)/V ; //for no early arrival
14:
if at < bC(i) and cg ≤ Q and bdt ≤ e0 then
15: 16: 17:
add C(i) into rV n ; //C(i) belongs to the current vehicle lt = bC(i) + sC(i) ; else if bC(i) ≤ at ≤ eC(i) and cg ≤ Q and bd2t ≤ e0 then
18:
add C(i) into rV n ; //C(i) belongs to the current vehicle
19:
lt = at + sC(i) ;
20: 21:
else if at > eC(i) or cg > Q or bd1t > e0 or bd2t > e0 then //C(i) isn’t served by the current vehicle
22:
add v0 into rV n ; //add depot as the ending of route
23:
Vn = Vn + 1; //C(i) is served by the new vehicle
24:
add v0 into rV n ; //add depot as the beginning of route
25:
add C(i) into rV n ;
26:
lt = bC(i) + sC(i) ;
27: 28:
cg = gC(i) ; end if
29: end for
21
4.4.2. Crossover Operator Since the permutation code is adopted, HMOEA-GL uses order crossover 430
(OX) as the crossover operator to perform genetic recombination process Davis (1991). An example of order crossover for permutation code is shown in figure 4. Parent 1 Offspring 1
1 2 3 4 5 6 7 8 9
2 1 3 4 5 6 7 8 9 5 4 6 9 2 1 7 8 3 Parent 2
Figure 4: An example of order crossover.
4.4.3. Mutation Operator Exchange mutation (EM) is the mutation operator of HMOEA-GL Banzhaf 435
(1990), which aims to increase the search ability of the algorithm by causing mutations of the individual. In order to increase the intensity of variation, this paper turns the exchange of gene position from the single into fragments. An example of EM for permutation code is shown in figure 5. 1 2 3 4 5 6 7 8 9
1 7 8 4 5 6 2 3 9
Figure 5: An example of exchange mutation.
4.5. Simple Insertion Search 440
As we all know, the number of vehicles is one of the most important objectives in VRP. Whether in the initial VRP or in variants of various VRP problems, it is necessary to reduce the number of vehicles as much as possible. In reality, adding one vehicle is very likely to make other objectives getting worse such as a significant increase in cost, and it is meaningless that a solution
445
which requires a large number of vehicles. Therefore, this paper uses SIS to 22
reduce the number of vehicles used in one solution. The use of SIS for individuals generated by genetic operation can be seen as a secondary mutation of the individual. Figure 6 shown an example of the main mechanism of SIS. [5, 10] [3, 7]
[3, 7]
5 1.3
[9, 15]
1.5 2
0.7 [9, 15]
2
1 [2, 6]
1
1
1
0
1
4
1
5
2
[2, 6]
1
[5, 10]
4
1
1
3
1.3
0
1
1
3
1.3
[2, 7] 1.5
[2, 7] 1.5
6
6
[2, 5]
[2, 5]
Figure 6: An example of the main mechanism of simple insertion search.
Among them, 0 is the depot, and 1 to 6 are six customers. Each customer 450
has a time window indicated by square brackets, and needs 1 service time. The distance between the customers determines the travel time of the vehicle. Taking figure 6 as an example, figure 7 shows the detailed process of SIS, and the pseudo code is shown in algorithm 3. SIS is used to insert the next customer into a certain route of the current
455
vehicle. Chromosome {1, 2, 4, 5, 3, 6}, which direct decoding will result in a solution that requires 3 vehicles. SIS judges whether the customer 2 can be served before the customer 1, and does not affect the service request of customer 1. Assume that the vehicle first serves customer 2, the end time is 10, and takes 2 units of time to go to customer 1. When the vehicle arriving at customer
460
1, the time is already 12, which has exceeded the time window of customer 1. Obviously, before the customer 2 is placed in the customer 1, the service of the customer 1 will be abnormal, so give up inserting the customer 2 between the depot and the customer 1. Customer 1 is already the last service customer in the current sub-route, then judges if customer 2 can be served after customer
465
1. Similarly, judges in turn whether the customer 4 can be inserted before the customer 1 or before the customer 2. It can be seen when the vehicle completes 23
Chromosome:
1
2
4
5
3
6
1
2
4
5
3
6
0
0
4
5
0
0
Direct decoding
0
1
2
Vehicle_2
Vehicle_1
3
0
6
Vehicle_3
Simple insertion search
x 0
2
1
0
√
1
x 0
1
2
0
1
2
0
1
4
√
2
0
4
1
2
4
2
x
0
1
0
1
4
2 √
5
4
2
x x 3
5
0
0
1
5
4
1
5
4
2
3
6
2
0
0
3
.. . 0
1
5
4
2
0
0
Vehicle_1
0
Vehicle_2
Figure 7: An example of the detailed process of Simple Insertion Search.
servicing customer 1, the end time is 2, after that, takes 2 units of time to go to customer 4, and spends 1 unit time waiting for customers 4 to start receiving the goods, the end time after the service for the customer 4 is 6. At this time, 470
take 0.7 units of time to go to customer 2, and it can be found that the service of the customer 2 is not affected. Therefore, it is feasible to insert the operation of the customer 4 between the customer 1 and the customer 2. In the same way, insert customer 5 before customer 4. In the end, the solution obtained through a simple insertion search requires only 2 vehicles and the number of vehicles has
475
been reduced. 4.6. Route Sequence Difference-based Local Search Differential evolution is an efficient evolutionary algorithm that guides the evolution of individuals based on differences among individuals. HMOEA-GL 24
Algorithm 3 Simple Insertion Search Input: VerticesData x, y, b, e, s, g; Distance d(i, j)=Euclidean(vi , vj ); VehicleData velocity V =1, capacity Q; Chromosome Sequence C; Output: Chromosome after Simple Insertion Search simC;
1: cg =0; // cg is current route load 2: aC(0) = bC(0) ; // at is arrival time of vertex t 3: lC(0) = bC(0) + sC(0) ; //lt is leave time of vertex t 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20: 21: 22: 23: 24: 25: 26: 27: 28: 29: 30: 31: 32: 33: 34: 35: 36: 37: 38:
Vn = 0; // Vn is the vehicle number add v0 into r; // add the depot as the beginning of route Vn = Vn + 1;
rVn = v0 ;
add C(0) into r; // add the first vertex in C into first route r Vn = Vn + 1;
rVn = C(0);
cg = cg + gC(0) ;
for i = 1 to length of C do f lag = 0; //flag to judge insert or not,1 is yes,0 is no for j = 1 to Vn − 1 do ctemp = cg + gC(i) ; g
atemp = lrj + d(rj , C(i))/V t
if atemp ≤ bC(i) and ctemp ≤ Q then g t lttemp = bC(i) + sC(i) ;
else if bC(i) < atemp ≤ eC(i) and ctemp ≤ Q then g t lttemp = atemp + sC(i) ; t
else if atemp > eC(i) and ctemp ≤ Q then g t continue; else if ctemp > Q then //This vertex cannot be serviced in r g break; end if atemp2 = lttemp + d(C(i), rj+1 )/V ; t if atemp2 ≤ brj+1 then t Vn = Vn + 1 add C(i) into r between rj and rj+1 ; cg = cg + gC(i) ;
lC(i) = lttemp ;
flag = 1;
break; else if atemp2 > brj+1 then t continue; end if end for if flag = 1 then
continue;
else if flag = 0 then //Determine if C(i) can be serviced at the end of the r if C(i) can be serviced in r then add C(0) into the end of r; else if C(i) can not be serviced in r then record the sequence of r into simC;
39: use a new vehicle to serve customer C(i), 40: end if 41: end if 42: end for
25
borrowed the idea of DE to further improve the performance of the algorithm. 480
With appropriate evolutionary operations, FSS-GS can find effective solutions fast. To further improve the quality of solutions, HMOEA-GL uses a local search called route sequence difference-based local search. RSD-LS is an additional search operation on the population generated by FSS-GS during the iteration of HMOEA-GL. The main idea of RSD-LS can be roughly divided into four
485
steps. • Step 1: Randomly pick an individual from the population. • Step 2: Select another individual. If it is better or worse than the individual selected in Step 1, then execute Step 3, otherwise repeat Step 2.
490
• Step 3: Use exchange order to determine the difference of route sequence between individuals selected in Step 1 and Step 2. • Step 4: Taking a certain proportion of the difference sequence obtained by Step 3, then apply it to the individual who performs poorly. Edge area-based sampling strategy
min f2
Center area-based elite sampling strategy
DLS
(S1 – S2)
S2
F(S1 – S2) S2 DLS
S 2` S 2`
S1
DLS S1 DLS DLS min f1
Figure 8: An example of Route Sequence Difference-based Local Search.
26
Taking figure 8 as an example, according to the above expression, RSD-LS 495
process can be approximated by a formula. S20 = S2 + F (S1 − S2 )
(11)
where S1 and S2 are two individual selected by Step 1 and Step 2. Suppose S2 is the poor one, that is, S2 is dominated by S1 . (S1 -S2 ) is the difference of route sequence between individuals S1 and S2 . F is the proportion, and S20 is the new individual that evolved from S2 by RSD-LS. Each new individual 500
will be evaluated by PDDR-FF, and new individuals which superior to previous evolution will be accepted. Otherwise the original individual has to put into the next generation of the population. By using RSD-LS, the population can be closer to the Pareto front. Because HMOEA-GL uses permutation code, so HMOEA-GL uses exchange
505
sequence to determine the difference of route sequence between individuals. Taking figure 9 as an example to describe the process of how to use exchange order to define the differences of route sequence between two individuals for permutation code. S1 1 2 3 4 5 6 S2 3 2 4 5 1 6
S1
1 2 3 4 5 6
S 2'
3 2 4 5 1 6
S1
1 2 3 4 5 6
S 2'
1 2 4 5 3 6
S1
1 2 3 4 5 6
S 2'
1 2 3 5 4 6
S 2'
1 2 3 4 5 6
{(0,4)}
{(0,4),(2,4)}
{(0,4),(2,4),(3,4)} Over
Figure 9: An example of Exchange Order.
27
Suppose S1 = {1, 2, 3, 4, 5, 6}, S2 = {3, 2, 4, 5, 1, 6}, and calculating process 510
of the difference of route sequence between S1 and S2 is as follows. Starting from the first gene position of the chromosome, compare each gene of the two individuals. Mark the first position subscript (0, x) when there has a difference in S2 . Then find the target gene 1 in S2 , and mark the second position subscript (0, 4). Combining two position subscripts becomes an exchange sequence
515
{(0, 4)}. Subsequently, exchange genes in S2 to obtain a temporary chromosome
S20 . Continue to compare S20 and S1 until the last set of genes of the chromosome
is detected. Sequence {(0, 4), (2, 4), (3, 4)} is the exchange order belong to S2 , represents the difference of route sequence between S1 and S2 . The max length of the exchange sequence for S2 is 3, which can be processed by proportion F . 520
Proportion F is a random number which F ∈ [0, 1]. Suppose F = 0.5, then the length of the exchange sequence executed is 1 (multiply max length of exchange by proportion F and round down). That is, the actual execution sequence is {(0, 4)}, the new individual S20 is generated as S20 = {1, 2, 4, 5, 3, 6}. Based on the route sequence differences between individuals, RSD-LS guides
525
the individuals with poor performance close to the better one, and the search ability of the algorithm can be further improved. 5. Experimental Evaluation In this part, section 5.1 describes the overall setup of the experiment. Section 5.2 introduces the evaluation indicators used in this experiment. Section
530
5.3 discusses the parameter setting of the algorithms. Section 5.4 verifies the effectiveness of SIS. Section 5.5 shows the experimental results and some conclusions. 5.1. Experimental Setup To verify the performance of proposed HMOEA-GL, a publicly available
535
benchmark generated by Solomon is used in this experimental studies Solomon (1987). The problem number used in this experiment are R101-R104, C101C104, RC101-RC104, which is random and with short scheduling horizon. Each 28
benchmark gives information about the depot and 100 customers, including coordinate locations, demand, time window, and service time. In this experiment, 540
indices C, GD, IGD, HV , Spacing and Spread are used to evaluate the performance of the algorithms. Wilcoxon rank-sum test is used for statistical test between proposed algorithm and competing algorithms. In order to make the comparison of experiments more fairer, the parameters of each algorithm are adjusted and optimized before the final experiment. In addition, since the SIS
545
is used for each algorithm to reduce the number of vehicles in this experiment, the effect of the SIS is also verified. Experiments are conducted to compare the effects of the proposed algorithm, NSGA-II, SPEA2 and MOEA/D on 12 Solomon benchmark problems, each algorithm is executed 30 times and record the solutions in the Pareto approximation each time. The population size is 50
550
and the max number of generations is set to 500. Selection operator uses binary tournament selection, crossover operator uses order crossover, mutation operator is exchange mutation, the SIS applied to all algorithms. The experimental results are marked, the darker the cell background color, the better the result. 5.2. Performance Indices
555
Coverage (C). C(S1 , S2 ) is the percent of the individuals in S2 who are weakly dominated by S1 . The larger C(S1 , S2 ) is, the better S1 outperforms S2 in C. Generational Distance (GD). GD evaluates a solution set according to the shortest distance of each solution from the P F ∗ . The smaller GD for one
560
solution set is, the better it is in approaching P F ∗ . In this experiment, the P F ∗ of the problem consists of the best results that all algorithms run multiple times. Inverted Generational Distance (IGD). IGD calculates the distance from P F ∗ to the solution set. IGD mainly evaluates the convergence ability, but
565
when the solution set is approaching P F ∗ , smaller IGD requires the individuals being distributed similarly to P F ∗ . The smaller IGD might mean better convergence and distribution. In IGD, the P F ∗ of the problem is also consists 29
of the best results that all algorithms run multiple times. Hypervolume (HV ). HV is the size of the dominated space of the obtained 570
solution set. Larger HV means better convergence performance. The reference point of HV in this experiment is set as combining the worst fitness values in each objective, which found by all algorithms run multiple times. Spacing (SP ). SP is a distribution performance evaluation indicator which according to the closest distance between each solution and other solutions. A
575
smaller SP might mean better distribution. Spread (Spread). Spread calculates the average of the distance between consecutive solutions. At the same time, consider the smallest distance between edge points and boundary solutions. A smaller Spread means a better distribution and spread.
580
5.3. Parameter Settings This part mainly compares the influence of different crossover rate Cr , mutation rate Mr , and proportion F of RSD-LS on the HMOEA-GL. Table 1 shows the performance indices of the final solution set obtained by running different parameter HMOEA-GL under the Solomon benchmark R101 problem. Based
585
on previous experience, the crossover rate Cr is generally set large, and the mutation rate Mr is generally small. If the proportion F of RSD-LS is too large, it will easily produce the exact same solution, and resulting to fall into local optimum. So the value of Cr , Mr and F is set to: Cr ∈ {0.5, 0.6, 0.7, 0.8, 0.9}, Mr ∈ {0.1, 0.2, 0.3}, and RSD-LS proportion F ∈ {0.5, 0.6, 0.7, 0.8}.
590
It can be seen that the better performing combinations are mainly concentrated between experiments 25-28 and experiments 37-40. In experiments 25-28 the crossover rate Cr is 0.7, and the mutation rate Mr is 0.1. The crossover rate Cr of experiments 37-40 is 0.8, the mutation rate Mr is 0.1. Among them, from the number of times of get top three of evaluation indicator, experiments 37-40
595
is more than experiments 25-28. In experiments 25-28 and experiments 37-40, it can see that when the proportion F of RSD-LS is set to 0.8, HMOEA-GL can perform better. Therefore, the parameters of HMOEA-GL are set as follows: 30
Table 1: The result of HMOEA-GL in different parameter in R101 problem. Problem: R101 Experiment 1 Experiment 2 Experiment 3 Experiment 4 Experiment 5 Experiment 6 Experiment 7 Experiment 8 Experiment 9 Experiment 10 Experiment 11 Experiment 12 Experiment 13 Experiment 14 Experiment 15 Experiment 16 Experiment 17 Experiment 18 Experiment 19 Experiment 20 Experiment 21 Experiment 22 Experiment 23 Experiment 24 Experiment 25 Experiment 26 Experiment 27 Experiment 28 Experiment 29 Experiment 30 Experiment 31 Experiment 32 Experiment 33 Experiment 34 Experiment 35 Experiment 36 Experiment 37 Experiment 38 Experiment 39 Experiment 40 Experiment 41 Experiment 42 Experiment 43 Experiment 44 Experiment 45 Experiment 46 Experiment 47 Experiment 48 Experiment 49 Experiment 50 Experiment 51 Experiment 52 Experiment 53 Experiment 54 Experiment 55 Experiment 56 Experiment 57 Experiment 58 Experiment 59 Experiment 60
Cr
Mr 0.1
0.5
0.2
0.3
0.1
0.6
0.2
0.3
0.1
0.7
0.2
0.3
0.1
0.8
0.2
0.3
0.1
0.9
0.2
0.3
F 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8
GD 21.3388 23.8032 16.6762 20.0018 25.4418 23.7405 29.8896 23.0882 30.9756 46.1087 36.6538 33.4476 16.3872 16.8135 18.4615 20.8670 29.9787 23.9819 19.7299 23.5509 29.0993 28.0303 30.8049 30.9219 12.9247 14.1581 13.9804 14.2389 16.4505 20.6195 19.6137 20.1798 29.3523 20.9034 25.5161 26.9648 12.1302 19.7369 12.9911 14.5651 16.7877 16.8810 23.2771 18.0375 22.3219 23.4938 26.1277 25.1435 15.6213 16.1051 16.3387 15.6077 21.6154 20.6990 18.8295 15.9019 22.0410 30.1262 19.4764 24.5973
31
IGD 13.3350 13.7668 12.8990 14.6975 16.1033 15.1378 15.7111 16.7063 19.0228 19.9408 19.2888 17.5324 11.0455 10.9044 11.7491 12.2983 13.8058 14.1500 12.0492 13.5544 16.4089 15.7600 16.9685 17.3720 9.7634 9.9506 10.5774 9.0088 12.5619 12.3818 12.0495 12.9600 17.5648 15.3338 16.1236 14.0746 9.9663 10.2646 8.9186 9.3273 11.3039 12.0827 12.3263 11.3698 14.6250 15.6468 15.0392 15.1042 10.8904 9.8199 10.1662 10.2055 13.2986 12.9988 12.1309 12.4055 15.5344 16.6866 14.8417 16.0166
HV 19887.0763 20119.8256 20162.2998 19898.8391 19658.3691 19838.2241 19880.7286 19988.1011 19254.2900 19086.5119 19414.6001 19657.7081 20471.0208 20438.4427 20623.4639 20461.0787 20169.9439 20300.1074 20568.4782 20205.0562 19662.4597 19937.2317 19589.2903 19822.9547 20626.8089 20773.6158 20637.2858 21121.5829 20367.7308 20265.9197 20582.3587 20539.8800 19733.4465 20125.2448 19838.5651 20196.2071 20758.1035 20726.7729 20928.8921 21046.1259 20503.8451 20537.7070 20433.4750 20580.5499 20228.5421 19907.7994 20195.4369 20370.9136 20834.4784 20818.3227 20704.2913 21045.9968 20397.1851 20390.4750 20397.9926 20550.8149 19961.9319 19946.6018 20063.2120 19966.5271
Spacing 17.9338 19.3704 17.1619 18.7967 23.1657 19.4627 25.3967 16.8518 17.9079 34.9161 31.6266 31.1964 17.3784 18.4281 22.7160 17.6778 24.8959 17.3863 18.9225 18.1104 20.8584 22.6456 25.3131 24.0197 17.8313 14.6415 13.1678 21.6098 12.5120 18.8373 19.8368 17.8454 22.2331 17.4680 21.9535 26.0599 13.7202 24.9500 15.8211 16.2298 20.1938 16.5157 28.3139 17.3354 18.3644 13.0944 19.4031 21.9715 16.5970 20.0627 17.6052 20.0819 19.2909 18.6517 24.8851 14.6071 21.0580 23.0481 19.8246 22.2405
Spread 0.9233 0.9530 0.9434 0.9659 0.9802 0.9228 1.0051 0.9214 0.9589 1.0237 1.0553 1.0098 0.9672 0.9902 1.0284 0.9682 0.9934 0.9368 0.9630 0.9191 0.9543 0.9839 0.9893 1.0090 1.0008 0.9678 0.9599 0.9994 0.9471 0.9599 0.9531 0.9441 0.9623 0.9326 0.9525 1.0178 0.9659 0.9946 0.9666 0.9884 1.0280 0.9791 1.0006 0.9751 0.9876 0.9165 0.9588 0.9461 0.9838 0.9753 0.9754 1.0121 0.9761 0.9818 0.9969 0.9025 0.9306 0.9196 0.9389 0.9628
crossover rate Cr =0.8, mutation rate Mr =0.1, RSD-LS proportion F =0.8. Table 2: The result of NSGA-II in different parameter in R101 problem. Problem: R101
Cr
Mr
GD
IGD
HV
Spacing
Spread
0.1
45.8989
24.9881
16632.5620
28.0686
0.9587
0.2
49.4310
28.7195
16619.2598
25.9404
0.9556
Experiment 3
0.3
54.5880
28.5689
16229.2180
24.7489
0.9457
Experiment 4
0.1
46.9625
28.5377
16873.7573
20.3073
0.9208
0.2
55.2580
30.1363
16406.4515
20.9027
0.9321
0.3
74.7008
34.9182
15672.5058
27.7712
0.9915
0.1
69.6310
37.5529
16079.3407
28.8825
0.9664
0.2
71.8652
39.6551
15995.1530
23.8254
0.9201
0.3
91.9884
47.0863
15185.3113
17.9976
0.8885
Experiment 1 Experiment 2
Experiment 5
0.7
0.8
Experiment 6 Experiment 7 Experiment 8
0.9
Experiment 9
Table 3: The result of SPEA2 in different parameter in R101 problem. Problem: R101
Cr
Mr
GD
IGD
HV
Spacing
Spread
0.1
22.3897
13.5210
18759.3233
14.4591
0.9351
0.2
29.4389
16.8471
18268.8967
17.8439
0.9539
Experiment 3
0.3
41.8934
18.5128
18010.6933
22.3111
0.9996
Experiment 4
0.1
27.0167
14.6118
18541.1259
11.2576
0.8937 0.9195
Experiment 1 Experiment 2
Experiment 5
0.7
0.2
27.6127
16.4009
18586.4625
15.7125
Experiment 6
0.3
46.3519
19.9279
17565.6125
14.3577
0.9243
Experiment 7
0.1
32.4969
15.4946
18548.8857
14.3954
0.9058
0.2
36.3992
16.8760
18342.2598
17.4023
0.9537
0.3
42.5555
17.7649
17957.3059
12.8516
0.9428
Experiment 8 Experiment 9
0.8
0.9
Similarly, for the competing algorithms, table 2 to 4 shows the performance 600
of NSGA-II, SPEA2 and MOEA/D with different parameters under Solomon benchmark R101 problem. It can be seen that for NSGA-II and SPEA2, the better parameters are also set to crossover rate Cr =0.8 and mutation rate Mr =0.1. For MOEA/D, can be found that MOEA/D performs better in the case of crossover rate Cr =0.9 and mutation rate Mr =0.1. In order to make a more fair
605
comparison, MOEA/D under two parameters MOEA/D-a(Cr =0.8, Mr =0.1) and MOEA/D-b(Cr =0.9, Mr =0.1) is set in the following experiments. 5.4. Effect of the SIS Component Section 4.5 refers to a local search method called SIS that can reduce vehicle number. In this part, experiments are conducted between HMOEA-GL with SIS 32
Table 4: The result of MOEA/D in different parameter in R101 problem. Problem: R101
Cr
Mr
GD
IGD
HV
Spacing
Spread
0.1
15.4365
15.2818
14697.7947
14.0840
0.9604
0.2
13.4091
14.2141
15218.6083
15.0232
0.9909
0.3
11.6031
12.6066
16047.6508
12.0886
0.9466
0.1
12.6256
14.8540
15002.5816
12.3173
0.9770
0.2
10.4925
12.1173
16308.3994
10.6865
0.9208
Experiment 6
0.3
9.5854
11.9201
16612.0227
11.3378
0.9655
Experiment 7
0.1
7.8974
12.9815
16285.5243
5.8513
0.9284
0.2
9.1806
12.9264
16228.4371
10.3583
0.9744
0.3
9.3730
12.1153
16516.0961
11.8013
0.9653
Experiment 1 Experiment 2
0.7
Experiment 3 Experiment 4 Experiment 5
Experiment 8
0.8
0.9
Experiment 9
610
and without SIS to verify the effectiveness of SIS. The parameters of this part are set as follows: Cr =0.8 , Mr =0.1, F =0.8. Table 5 shows the average value and best value of the vehicle number of the solution obtained by HMOEA-GL with SIS and without SIS. Table 5: The vehicle number of HMOEA-GL with SIS and without SIS HMOEA-GL with SIS
HMOEA-GL without SIS
Vehicle Nmuber Mean
Best
Mean
Best
R101
38.6498
31
39.1064
32
R102
34.3476
29
33.6790
29
R103
27.2037
23
27.0000
24
R104
21.3469
19
22.0000
20
C101
21.1860
16
23.7532
18
C102
20.8333
14
21.4593
15
C103
17.9293
13
18.1917
15
C104
13.3077
12
13.7805
12
RC101
33.3000
30
32.8857
30
RC102
30.3056
27
29.5366
27
RC103
26.1190
24
26.0000
24
RC104
22.9444
21
22.9744
21
It can be seen from table 5 that, for the average of vehicle number, in the 615
half of experimental problems HMOEA-GL with SIS is better than HMOEA-GL without SIS. In six problems, HMOEA-GL with SIS is able to find a solution with fewer vehicle number, and on the remaining problems it is the same as with HMOEA-GL without SIS. It can be shown that the addition of SIS can 33
effectively reduce the number of vehicles in the solution. In the following exper620
iments, SIS is added to all the algorithms participating in the comparison. 5.5. Comparisons with other Multiobjective Algorithms Experiment will compare HMOEA-GL, NSGA-II, SPEA2, MOEA/D-a, MOEA/Db and multiobjective evolutionary algorithm with FSS-GS but no RSD-LS (MOEAG). According to the results discussed in section 5.3, the parameters of algo-
625
rithms are set as table 6. Table 6: The parameters of algorithms Parameters
HMOEA-GL
MOEA-G
NSGA-II
SPEA2
MOEA/D-a
MOEA/D-b 50
Population Size
50
50
50
50
50
Elite Population Size
25
25
50
50
/
/
Sub-population Size
25
25
/
/
/
/
Crossover Rate
0.8
0.8
0.8
0.8
0.8
0.9
Mutation Rate
0.1
0.1
0.1
0.1
0.1
0.1
DLS proportion F
0.8
/
/
/
/
/
Number of Neighborhood T
/
/
/
/
5
5
Number of Sub-problems N
/
/
/
/
50
50
The average and standard deviation of convergence indices, distribution indices, and spread indices is shown in table 7 to 11, the result of the hypothesis test between HMOEA-GL and competing algorithms are expressed in the H.T. column. 630
From the GD index as shown in table 7, compared with the other algorithms, among the 12 questions tested, HMOEA-GL performed best on 8 problems and ranked second on 2 problems. MOEA-G performed best on 2 problems and ranked second on 4 problems. This shows that the final solution set obtained by HMOEA-GL and MOEA-G are closer to P F ∗ than the solution set obtained
635
by NSGA-II, SPEA2 and MOEA/D, which means that the convergence performance is better. The result of Wilcoxon rank-sum test can also support these conclusions. From the IGD index and the result of Wilcoxon rank-sum test as shown in table 8, HMOEA-GL performed best on 8 problems again and ranked second 34
Table 7: The experimental results of GD index HMOEA-GL
MOEA-G
NSGA-II
SPEA2
MOEA/D-a
MOEA/D-b
GD Mean
S.D.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
R101 1.46E+01 1.10E+01 1.30E+01 9.64E+00
0
4.70E+01 3.03E+01
1
2.70E+01 2.12E+01
1
1.26E+01 7.14E+00
0
1.09E+01 6.27E+00
0
R102 1.59E+01 1.44E+01 1.55E+01 1.15E+01
0
4.74E+01 2.66E+01
1
2.47E+01 1.55E+01
1
1.01E+01 4.79E+00
0
1.06E+01 6.85E+00
0
R103 1.18E+01 1.85E+01 1.56E+01 1.96E+01
0
4.01E+01 1.98E+01
1
1.66E+01 1.54E+01
1
1.87E+01 1.33E+01
1
1.51E+01 1.37E+01
1
R104 4.21E+00 5.31E+00 7.03E+00 7.28E+00
0
1.64E+01 1.34E+01
1
8.16E+00 5.77E+00
1
1.12E+01 9.40E+00
1
9.79E+00 8.43E+00
1
C101 1.75E+02 1.35E+02 2.21E+02 1.77E+02
0
1.17E+03 3.26E+02
1
4.52E+02 2.41E+02
1
1.96E+03 4.89E+02
1
2.13E+03 5.02E+02
1
C102 8.27E+01 6.79E+01 8.40E+01 5.79E+01
0
5.40E+02 2.29E+02
1
1.73E+02 1.43E+02
1
1.01E+03 3.15E+02
1
9.15E+02 4.99E+02
1
C103 4.51E+01 4.93E+01 8.31E+01 7.23E+01
1
1.68E+02 9.60E+01
1
7.71E+01 6.50E+01
1
3.52E+02 3.41E+02
1
3.74E+02 3.48E+02
1
C104 3.08E+01 4.21E+01 4.02E+01 2.78E+01
1
9.74E+01 7.20E+01
1
3.63E+01 3.94E+01
0
1.12E+02 1.72E+02
0
1.18E+02 1.66E+02
1
RC101 7.66E+00 5.55E+00 6.74E+00 7.10E+00
0
1.49E+01 7.37E+00
1
1.23E+01 5.42E+00
1
1.03E+01 3.14E+00
1
1.02E+01 2.66E+00
1
RC102 6.70E+00 7.13E+00 9.91E+00 1.16E+01
0
1.26E+01 1.52E+01
1
7.84E+00 4.41E+00
1
7.86E+00 3.81E+00
1
8.16E+00 3.08E+00
1
RC103 4.39E+00 5.32E+00 3.82E+00 3.13E+00
0
1.26E+01 1.24E+01
1
8.63E+00 4.17E+00
1
8.17E+00 5.54E+00
1
1.04E+01 7.24E+00
1
RC104 2.60E+00 2.23E+00 5.96E+00 9.09E+00
0
6.74E+00 4.50E+00
1
7.43E+00 3.89E+00
1
6.63E+00 1.53E+00
1
7.17E+00 2.84E+00
1
Table 8: The experimental results of IGD index HMOEA-GL
MOEA-G
NSGA-II
SPEA2
MOEA/D-a
MOEA/D-b
IGD Mean
S.D.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
R101 9.33E+00 4.37E+00 9.24E+00 3.29E+00
0
2.85E+01 1.01E+01
1
1.46E+01 5.10E+00
1
1.49E+01 2.16E+00
1
1.25E+01 1.79E+00
1
R102 5.03E+00 1.95E+00 6.56E+00 3.00E+00
1
1.51E+01 3.16E+00
1
9.63E+00 3.18E+00
1
9.96E+00 2.21E+00
1
9.06E+00 1.39E+00
1
R103 2.75E+00 1.50E+00 3.84E+00 2.09E+00
1
8.93E+00 2.37E+00
1
6.57E+00 1.36E+00
1
7.98E+00 1.57E+00
1
7.62E+00 1.87E+00
1
R104 1.51E+00 1.02E+00 2.22E+00 1.33E+00
1
4.97E+00 1.23E+00
1
5.48E+00 1.17E+00
1
5.89E+00 1.26E+00
1
5.76E+00 1.27E+00
1
C101 2.04E+02 8.98E+01 2.77E+02 1.41E+02
1
9.50E+02 2.40E+02
1
3.80E+02 1.92E+02
1
6.49E+01 2.43E+01
1
6.78E+01 2.89E+01
1
C102 3.99E+01 2.03E+01 4.55E+01 1.86E+01
0
3.23E+02 1.37E+02
1
9.80E+01 1.01E+02
1
5.29E+01 1.83E+01
1
5.54E+01 2.30E+01
1
C103 2.93E+01 1.38E+01 3.34E+01 1.53E+01
0
4.04E+01 1.49E+01
1
2.51E+01 1.18E+01
0
4.71E+01 1.00E+01
1
4.39E+01 1.22E+01
1
C104 1.65E+00 1.05E+00 3.07E+00 1.25E+00
1
4.03E+00 1.20E+00
1
3.27E+00 1.11E+00
1
4.64E+00 1.04E+00
1
4.62E+00 9.88E-01
1
RC101 4.68E+00 1.06E+00 4.24E+00 1.81E+00
0
1.06E+01 1.78E+00
1
8.79E+00 1.42E+00
1
9.69E+00 1.77E+00
1
9.85E+00 1.89E+00
1
RC102 2.60E+00 1.35E+00 2.99E+00 1.43E+00
0
6.49E+00 1.31E+00
1
6.37E+00 1.44E+00
1
7.24E+00 1.66E+00
1
7.05E+00 1.46E+00
1
RC103 2.17E+00 1.04E+00 2.19E+00 9.68E-01
0
5.87E+00 1.21E+00
1
7.38E+00 1.49E+00
1
6.90E+00 1.45E+00
1
7.40E+00 1.34E+00
1
RC104 1.99E+00 9.73E-01 2.35E+00 9.94E-01
0
5.57E+00 1.16E+00
1
6.60E+00 1.19E+00
1
6.45E+00 1.39E+00
1
6.47E+00 1.33E+00
1
35
640
on 4 problems. MOEA-G also performed best on 2 problems and ranked second on 7 problems. IGD index can not only judge the convergence performance of an algorithm, but also consider the performance of distribution performance. The results of IGD experiments show that the final solution set obtained by HMOEA-GL and MOEA-G is not only closer to P F ∗ , but more consistent with
645
the entire P F ∗ . Table 9: The experimental results of HV index
HMOEA-GL
MOEA-G
NSGA-II
SPEA2
MOEA/D-a
MOEA/D-b
HV Mean
S.D.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
R101 2.10E+04 9.38E+02 2.09E+04 6.80E+02
0
1.69E+04 9.81E+02
1
1.85E+04 9.77E+02
1
1.50E+04 9.71E+02
1
1.62E+04 1.04E+03
1
R102 1.52E+04 6.71E+02 1.55E+04 9.02E+02
0
1.22E+04 1.03E+03
1
1.30E+04 1.04E+03
1
1.08E+04 1.29E+03
1
1.15E+04 6.78E+02
1
R103 6.06E+03 4.22E+02 5.99E+03 4.30E+02
0
4.61E+03 3.95E+02
1
4.82E+03 4.16E+02
1
4.02E+03 5.10E+02
1
4.12E+03 5.90E+02
1
R104 1.00E+03 1.12E+02 9.80E+02 9.91E+01
0
6.48E+02 1.19E+02
1
5.66E+02 1.01E+02
1
5.28E+02 1.20E+02
1
5.25E+02 1.30E+02
1
C101 2.80E+05 7.63E+03 2.74E+05 8.59E+03
1
2.18E+05 1.20E+04
1
2.61E+05 1.02E+04
1
1.37E+05 2.23E+04
1
1.45E+05 2.00E+04
1
C102 1.65E+05 5.57E+03 1.64E+05 6.38E+03
0
1.34E+05 7.94E+03
1
1.56E+05 6.26E+03
1
9.25E+04 1.13E+04
1
9.66E+04 1.11E+04
1
C103 6.99E+04 2.57E+03 6.92E+04 1.77E+03
0
5.98E+04 2.81E+03
1
6.62E+04 2.43E+03
1
4.45E+04 3.73E+03
1
4.61E+04 4.36E+03
1
C104 6.91E+03 5.50E+02 7.03E+03 5.00E+02
0
5.43E+03 6.69E+02
1
6.20E+03 5.71E+02
1
3.15E+03 8.74E+02
1
3.29E+03 1.06E+03
1
RC101 1.63E+03 1.11E+02 1.67E+03 1.59E+02
0
1.09E+03 1.34E+02
1
1.22E+03 1.37E+02
1
1.10E+03 1.79E+02
1
1.08E+03 1.88E+02
1
RC102 2.56E+03 2.55E+02 2.52E+03 2.49E+02
0
1.66E+03 2.32E+02
1
1.64E+03 2.84E+02
1
1.41E+03 3.62E+02
1
1.49E+03 3.45E+02
1
RC103 1.45E+03 1.29E+02 1.48E+03 1.39E+02
0
9.86E+02 1.48E+02
1
7.81E+02 1.80E+02
1
8.09E+02 1.81E+02
1
7.81E+02 1.73E+02
1
RC104 6.38E+02 8.91E+01 6.16E+02 7.76E+01
0
3.30E+02 8.79E+01
1
2.28E+02 1.15E+02
1
2.35E+02 1.21E+02
1
2.43E+02 1.07E+02
1
From the HV index and the result of statistical test as shown in table 9, HMOEA-GL still ranked first on 8 problems and ranked second on 4 other problems. MOEA-G achieved the first place on 4 problems and ranked second of other 8 problems. The larger the HV value, the larger the dominance area, 650
which means that the final solution set obtained by HMOEA-GL and MOEA-G have more larger dominance area than the solution set obtained by NSGAII, SPEA2 and MOEA/D. That is, HMOEA-GL and MOEA-G have better performance in convergence and distribution. It can be seen that HMOEA-GL and MOEA-G are not outstanding in the
655
Spacing index and Spread index as shown in table 10 and 11. Spacing index calculates the sum of the variance between the distance of successive solutions and the shortest distance in the entire solution set. It is mainly used to judge the distribution performance of the solution set, whether the distribution of the solution is uniform. Spread not only calculates the distance between successive
36
Table 10: The experimental results of Spacing index HMOEA-GL
MOEA-G
NSGA-II
SPEA2
MOEA/D-a
MOEA/D-b
Spacing Mean
S.D.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
R101 1.62E+01 1.59E+01 1.64E+01 2.08E+01
0
2.03E+01 1.94E+01
0
1.13E+01 7.63E+00
0
1.23E+01 1.38E+01
0
1.26E+01 1.19E+01
0
R102 1.72E+01 2.20E+01 1.29E+01 1.42E+01
0
1.55E+01 1.34E+01
0
1.11E+01 1.06E+01
0
5.58E+00 6.30E+00
1
6.88E+00 1.18E+01
1
R103 1.47E+01 3.05E+01 1.77E+01 3.60E+01
0
2.40E+01 2.45E+01
1
1.03E+01 1.34E+01
0
1.66E+01 2.37E+01
0
1.49E+01 1.70E+01
0
R104 1.68E+00 5.76E+00 3.63E+00 6.31E+00
1
7.83E+00 1.43E+01
1
3.29E+00 9.64E+00
0
6.11E+00 1.15E+01
1
5.21E+00 1.30E+01
0
C101 1.13E+02 8.51E+01 1.18E+02 8.62E+01
0
1.80E+02 1.79E+02
0
1.01E+02 5.28E+01
0
3.39E+02 3.93E+02
1
2.52E+02 1.33E+02
1
C102 1.12E+02 1.05E+02 1.19E+02 9.68E+01
0
1.46E+02 9.59E+01
0
8.30E+01 6.15E+01
0
4.92E+02 5.79E+02
1
4.25E+02 5.45E+02
1
C103 1.01E+02 9.23E+01 1.12E+02 1.18E+02
0
1.45E+02 1.32E+02
0
6.60E+01 5.71E+01
0
3.95E+02 4.58E+02
0
4.39E+02 4.16E+02
1
C104 1.78E+01 2.96E+01 4.38E+01 4.47E+01
1
9.52E+01 1.09E+02
1
2.31E+01 3.67E+01
0
9.41E+01 1.61E+02
0
6.46E+01 1.18E+02
0
RC101 3.89E+00 5.87E+00 3.78E+00 6.08E+00
0
5.63E+00 1.10E+01
0
1.42E+00 4.46E+00
1
1.31E+00 3.19E+00
0
3.11E-01 7.65E-01
1
RC102 5.48E+00 1.03E+01 8.36E+00 1.22E+01
0
9.05E+00 2.91E+01
0
2.44E+00 4.94E+00
0
3.92E+00 6.09E+00
0
1.08E+00 2.64E+00
0
RC103 2.55E-02 8.47E-02 2.18E+00 4.53E+00
0
4.62E+00 1.75E+01
0
1.03E+00 2.06E+00
0
1.32E+00 3.97E+00
0
4.41E-01 1.76E+00
0
RC104 3.45E+00 7.71E+00 1.03E+00 3.67E+00
0
4.06E+00 1.10E+01
0
4.22E-02 8.44E-02
0
0.00E+00 0.00E+00
0
0.00E+00 0.00E+00
0
Table 11: The experimental results of Spread index HMOEA-GL
MOEA-G
NSGA-II
SPEA2
MOEA/D-a
MOEA/D-b
Spread Mean
S.D.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
R101
9.88E-01 1.66E-01 9.55E-01 2.02E-01
0
9.21E-01 1.17E-01
0
8.94E-01 1.14E-01
1
9.77E-01 1.41E-01
0
9.34E-01 1.30E-01
R102
9.31E-01 2.01E-01 9.20E-01 2.03E-01
0
9.53E-01 1.14E-01
0
8.59E-01 1.35E-01
0
8.81E-01 1.03E-01
0
8.67E-01 1.03E-01
0 0
R103
8.23E-01 2.34E-01 8.99E-01 1.84E-01
0
9.00E-01 1.65E-01
0
8.59E-01 1.46E-01
0
8.57E-01 1.99E-01
0
9.01E-01 1.05E-01
0
R104
6.13E-01 1.39E-01 7.34E-01 1.59E-01
1
7.44E-01 1.99E-01
1
7.55E-01 1.56E-01
1
8.22E-01 1.33E-01
1
7.57E-01 1.54E-01
1
C101
9.26E-01 6.81E-02 9.44E-01 1.02E-01
0
9.26E-01 8.16E-02
0
8.82E-01 8.02E-02
1
1.10E+00 1.14E-01
1
1.11E+00 1.00E-01
1
C102
9.90E-01 1.61E-01 9.85E-01 1.43E-01
0
9.37E-01 8.47E-02
0
8.90E-01 1.24E-01
1
1.15E+00 1.14E-01
1
1.08E+00 2.17E-01
1
C103 1.06E+00 2.47E-01 1.03E+00 1.53E-01
0
9.91E-01 1.25E-01
1
9.26E-01 1.47E-01
1
1.01E+00 3.03E-01
0
1.01E+00 2.24E-01
0
C104
6.69E-01 1.99E-01 8.93E-01 1.92E-01
1
8.33E-01 2.12E-01
1
7.46E-01 1.72E-01
0
7.08E-01 1.80E-01
0
6.48E-01 1.85E-01
0
RC101 7.89E-01 1.49E-01 7.95E-01 2.20E-01
0
8.18E-01 1.37E-01
0
7.30E-01 1.51E-01
0
8.76E-01 1.10E-01
1
8.26E-01 1.05E-01
0
RC102 7.46E-01 2.36E-01 8.04E-01 1.99E-01
0
7.60E-01 1.62E-01
0
7.90E-01 1.27E-01
0
7.43E-01 1.68E-01
0
7.82E-01 1.24E-01
0
RC103 6.08E-01 1.07E-01 7.14E-01 1.45E-01
1
6.79E-01 1.51E-01
0
7.47E-01 1.33E-01
1
8.74E-01 1.35E-01
1
7.10E-01 1.24E-01
1
RC104 7.09E-01 1.74E-01 6.69E-01 1.41E-01
0
7.68E-01 1.32E-01
0
6.08E-01 4.14E-02
0
8.18E-01 1.15E-01
0
7.21E-01 1.25E-01
0
37
660
solutions, but also considers the smallest distance between edge points and boundary solutions. From Spacing index, HMOEA-GL only performed best on 3 problems and ranked second on 4 problems. MOEA-G ranks third on 6 problems. From the Spread index, HMOEA-GL performed best on 3 problems and ranked second on 3 other problems. MOEA-G ranked second on 2 problems
665
and ranked third on 2 problems. The VRPTW discussed in this paper is not a continuous function problem. We found that when the number of vehicles becomes small, the number of feasible solutions existing in the solution space will decrease. As the solution space approaches the direction in which the number of vehicles decreases, the difference in total waiting time between solutions for
670
different vehicle numbers increases. However, HMOEA-GL and MOEA-G are more prefers to look for the solutions which fewer vehicles number, which may be the reason for their general performance on the Spacing index and Spread index. HMOEA-GL does not use special density maintaining strategies (such as calculate crowding distance) to maintain the distribution of the population,
675
and it is acceptable to obtain such distribution performance. Overall, comparing to NSGA-II, SPEA2 and MOEA/D, HMOEA-GL and MOEA-G have better convergence performance in terms of 12 Solomon benchmark problems. The distribution performance of HMOEA-GL and MOEA-G is acceptable. It can be stated that the proposed FSS-GS and the framework of
680
HMOEA-GL are effective. From the experimental results, it can be found that the significance between HMOEA-GL and MOEA-G is not strong, so in order to more accurately judge the ability of RSD-LS, experiment evaluated the C index for HMOEA-GL and MOEA-G. C index is a binary performance index, can directly reflect the convergence performance of the two algorithms. The
685
results of the C index is shown in table 12. It can be seen that HMOEA-GL performs better on 10 problems than MOEA-G on the C index. This means, in most problems, the final solution set obtained by HMOEA-GL dominates more the solution obtained by MOEAG. It can be concluded that RSD-LS is effective.
690
In summary, HMOEA-GL has been verified on 12 Solomon benchmark prob38
Table 12: Mean and standard deviation of C between HMOEA-GL and MOEA-G C
C(HM OEA−GL,M OEA−G) Mean
C(M OEA−G,HM OEA−GL)
S.D.
Mean
S.D.
R101
0.4323
3.75E-01
0.4908
3.91E-01
R102
0.5898
4.07E-01
0.3526
3.89E-01
R103
0.5931
4.39E-01
0.3220
4.18E-01
R104
0.5006
4.68E-01
0.2722
3.93E-01
C101
0.6195
4.52E-01
0.3578
4.30E-01
C102
0.5479
3.68E-01
0.3889
3.79E-01
C103
0.6860
3.39E-01
0.2671
3.35E-01
C104
0.6719
3.90E-01
0.2194
3.35E-01
RC101
0.3506
3.98E-01
0.5322
4.42E-01
RC102
0.4178
4.66E-01
0.4083
4.62E-01
RC103
0.3556
4.79E-01
0.3500
4.76E-01
RC104
0.4278
4.79E-01
0.2389
4.14E-01
lems with four algorithms, MOEA-G, NSGA-II, SPEA2 and MOEA/D, by five performance indices GD, IGD, HV , Spacing and Spread. And conducted an additional evaluation between HMOEA-GL and MOEA-G on the C index. The experimental results show that HMOEA-GL have more excellent performance 695
in convergence, while maintaining a satisfying distribution performance.
6. Conclusions In this paper, HMOEA-GL is proposed for solving VRPTW. Two objectives are mainly considered, one is minimizing the time-wasting due to early arrival, and another is reducing the number of vehicle. In HMOEA-GL, FSS-GS is de700
signed as the global search strategy. FSS-GS mixes the center area-based elite sampling strategy and edge area-based sampling strategy. Center area-based elite sampling strategy can fast determine which individuals are nondominated or dominated but have large domination area, i.e. an individual located in the central area of the Pareto front. Combining the edge area-based sampling strat-
705
egy which prefers the edge region on the Pareto front, FSS-GS increases the attention to individuals which located on the edge and center regions of the Pareto frontier. FSS-GS can guarantee the capability of HMOEA-GL for con39
verging to the multiply directions of Pareto front quickly and its natural uniform distribution performance. RSD-LS is implemented as a local search strategy af710
ter FSS-GS. RSD-LS can further improve the search ability of HMOEA-GL by guiding the individuals with poor performance to move closer to better one. Permutation code is used to represent the solution of VRPTW, and use the appropriate genetic operator to genetic iteration. In the case of using permutation code, exchange order is used to indicate the difference of route sequence between
715
individuals, which support the execution of RSD-LS smoothly. Designing suitable coding method to represent the individual of VRPTW, and proper genetic operators are used to recombine chromosomes. A single objective optimization method SIS is used to reduce the number of vehicles in one solution. Comparing with NSGA-II, SPEA2, MOEA/D on the 12 Solomon benchmark problems,
720
the experimental results show that the HMOEA-GL hybridizing FSS-GS and RSD-LS is effective. As an effectively intelligent algorithm, HMOEA-GL can be the core of an expert and intelligent decision support system, and provide suitable solutions for decision-makers. However, the proposed method might lose the efficiency in
725
solving the hyperdimensional optimization problems, since when the number of objectives is more than four, the sampling strategy of VEGA in FSS-GS will be ineffective and the computation time of calculating the fitness by PDDR-FF will increase considerably. In the future work, it is appropriate to consider some factors of the problems such as multiple types of goods, multiple types of vehicles,
730
multiple depots and green energy etc. In addition, hybridizing HMOEA-GL with other effective algorithms such as ACO, PSO, TS and SA etc., might be an interesting research direction.
Acknowledgements This research work is supported by the Science & Technology Research 735
Project of Henan Province (162102210044), Program for Science & Technology Innovation Talents in Universities of Henan Province (19HASTIT027), Key Re-
40
search Project in Universities of Henan Province (17A520030), and the Grantin-Aid for Scientific Research (C) of Japan Society of Promotion of Science (JSPS) (19K12148).
740
References Abidi, H., Hassine, K., & Mguis, F. (2018). Genetic algorithm for solving a dynamic vehicle routing problem with time windows. In 2018 International Conference on High Performance Computing & Simulation (HPCS) (pp. 782– 788). IEEE.
745
Baˇ nos, R., Ortega, J., Gil, C., Fernndez, A., & Toro, F. D. (2013). A simulated annealing-based parallel multi-objective approach to vehicle routing problems with time windows. Expert Systems with Applications, 40 , 1696–1707. Banzhaf, W. (1990). The “molecular” traveling salesman. Biological Cybernetics, 64 , 7–14.
750
Braekers, K., Ramaekers, K., & Nieuwenhuyse, I. V. (2016). The vehicle routing problem: State of the art classification and review. Computers & Industrial Engineering, 99 , 300–313. Dantzig, G. B., & Ramser, J. H. (1959). The truck dispatching problem. Management science, 6 , 80–91.
755
Davis, L. (1991). Handbook of genetic algorithms, . Desrochers, M., Desrosiers, J., & Solomon, M. (1992). A new optimization algorithm for the vehicle routing problem with time windows. Operations research, 40 , 342–354. Dixit, A., Mishra, A., & Shukla, A. (2019). Vehicle routing problem with time
760
windows using meta-heuristic algorithms: A survey. In Harmony Search and Nature Inspired Optimization Algorithms (pp. 539–546). Springer.
41
Garcia-Najera, A., & Bullinaria, J. A. (2011). An improved multi-objective evolutionary algorithm for the vehicle routing problem with time windows. Computers & Operations Research, 38 , 287–300. 765
Gen, M., & Cheng, R. (2000). Genetic algorithms and engineering optimization. John Wiley & Sons. Gen, M., Cheng, R., & Lin, L. (2008). Network models and optimization: Multiobjective genetic algorithm approach. Springer Science & Business Media. G¨o¸cken, T., Yaktubay, M., & Kili¸c, F. (2017). Improvement of a genetic algo-
770
rithm approach for the solution of vehicle routing problem with time windows. In Artificial Intelligence and Data Processing Symposium (IDAP), 2017 International (pp. 1–8). IEEE. Gonz´ alez, O. M., Segura, C., Pe˜ na, S. I. V., & Le´ on, C. (2017). A memetic algorithm for the capacitated vehicle routing problem with time windows. In
775
Evolutionary Computation (CEC), 2017 IEEE Congress on (pp. 2582–2589). IEEE. Gupta, A., & Saini, S. (2017). An enhanced ant colony optimization algorithm for vehicle routing problem with time windows. In 2017 Ninth International Conference on Advanced Computing (ICoAC) (pp. 267–274). IEEE.
780
Holland, J. H. (1992). Adaptation in natural and artificial systems. Ann Arbor , 6 , 126137. Jiang, X., Zhou, Q., & Ye, Y. (2017). Method of task assignment for uav based on particle swarm optimization in logistics. In Proceedings of the 2017 International Conference on Intelligent Systems, Metaheuristics & Swarm
785
Intelligence (pp. 113–117). ACM. Lenstra, J. K., & Kan, A. H. G. R. (1981). Complexity of vehicle routing and scheduling problems. Networks, 11 , 221–227.
42
Liu, S. C., & Chen, J. R. (2011). A heuristic method for the inventory routing and pricing problem in a supply chain. Expert Systems with Applications, 38 , 790
13223–13231. Long, J., Sun, Z., Pardalos, P. M., Hong, Y., Zhang, S., & Li, C. (2019). A hybrid multi-objective genetic local search algorithm for the prize-collecting vehicle routing problem. Information Sciences, 478 , 40–61. Marinakis, Y., Marinaki, M., & Migdalas, A. (2019). A multi-adaptive parti-
795
cle swarm optimization for the vehicle routing problem with time windows. Information Sciences, 481 , 311–329. Moon, I. K., Lee, J.-H., & Seong, J. (2012). Vehicle routing problem with time windows considering overtime and outsourcing vehicles. Expert Systems with Applications, 39 , 13202–13213.
800
Pratiwi, A., Pratama, A., Sa’diyah, I., & Suprajitno, H. (2018). Vehicle routing problem with time windows using natural inspired algorithms. In Journal of Physics: Conference Series (p. 012025). IOP Publishing volume 974. Pu, E., Wang, F., Yang, Z., Wang, J., Li, Z., & Huang, X. (2017). Hybrid differential evolution optimization for the vehicle routing problem with time
805
windows and driver-specific times. Wireless Personal Communications, 95 , 2345–2357. Salamatbakhsh-Varjovi, A., Tavakkoli-Moghaddam, R., Alinaghian, M., & Najafi, E. (2018). Robust periodic vehicle routing problem with time windows under uncertainty: An efficient algorithm. KSCE Journal of Civil Engineer-
810
ing, 22 , 4626–4634. Schaffer, J. D. (1985). Multiple objective optimization with vector evaluated genetic algorithms. In Proceedings of the First International Conference on Genetic Algorithms and Their Applications, 1985 . Lawrence Erlbaum Associates. Inc., Publishers.
43
815
Schneider, M. (2016). The vehicle-routing problem with time windows and driver-specific times. European Journal of Operational Research, 250 , 101– 119. Solomon, M. M. (1987). Algorithms for the vehicle routing and scheduling problems with time window constraints. Operations research, 35 , 254–265.
820
Storn, R., & Price, K. (1997). Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces. Journal of global optimization, 11 , 341–359. Tang, L., & Wang, X. (2013). A hybrid multiobjective evolutionary algorithm for multiobjective optimization problems. IEEE Transactions on Evolutionary
825
Computation, 17 , 20–45. Tasan, A. S., & Gen, M. (2012). A genetic algorithm based approach to vehicle routing problem with simultaneous pick-up and deliveries. Computers & Industrial Engineering, 62 , 755–761. Wang, Y., & Xing, L. (2018). A multi ant system based hybrid heuristic algo-
830
rithm for vehicle routing problem with service time customization. In International Conference on Bio-Inspired Computing: Theories and Applications (pp. 411–422). Springer. Wang, Z., Li, Y., & Hu, X. (2015). A heuristic approach and a tabu search for the heterogeneous multi-type fleet vehicle routing problem with time windows
835
and an incompatible loading constraint. Computers & Industrial Engineering, 89 , 162–176. Wolpert, D. H., Macready, W. G. et al. (1997). No free lunch theorems for optimization. IEEE transactions on evolutionary computation, 1 , 67–82. Yong, S., Boudouh, T., & Grunder, O. (2016). A hybrid genetic algorithm for
840
a home health care routing problem with time window and fuzzy demand. Expert Systems with Applications, 72 , 160–176.
44
Yu, B., Yang, Z., & Yao, B. (2011). A hybrid algorithm for vehicle routing problem with time windows. Expert Systems with Applications, 38 , 435–441. Yu, X., & Gen, M. (2010). Introduction to evolutionary algorithms. Springer 845
Science & Business Media. Zhang, H., Zhang, Q., Ma, L., Zhang, Z., & Liu, Y. (2019). A hybrid ant colony optimization algorithm for a multi-objective vehicle routing problem with flexible time windows. Information Sciences, 490 , 166–190. Zhang, J., Yang, F., & Weng, X. (2018). An evolutionary scatter search parti-
850
cle swarm optimization algorithm for the vehicle routing problem with time windows. IEEE Access, 6 , 63468–63485. Zhang, W., Gen, M., & Jo, J. (2014). Hybrid sampling strategy-based multiobjective evolutionary algorithm for process planning and scheduling problem. Journal of Intelligent Manufacturing, 25 , 881–897.
855
Zhou, A., Qu, B. Y., Li, H., Zhao, S. Z., Suganthan, P. N., & Zhang, Q. (2011). Multiobjective evolutionary algorithms: A survey of the state of the art. Swarm & Evolutionary Computation, 1 , 32–49.
45
Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:
Author Contribution Statement Wenqiang Zhang: Conceptualization, Methodology, Writing – Review & Editing, Project Administration, Funding Acquisition, Supervision. Diji Yang: Data Curation, Software, Resources, Visualization, Formal Analysis. Guohui Zhang: Validation, Investigation, Writing – Review & Editing. Mitsuo Gen: Visualization, Validation, Writing – Review & Editing, Supervision.
Journal Pre-proof
Hybrid Multiobjective Evolutionary Algorithm with Fast Sampling Strategy-based Global Search and Route Sequence Difference-based Local Search for VRPTW Wenqiang Zhang, Diji Yang, Guohui Zhang, Mitsuo Gen PII: DOI: Reference:
S0957-4174(19)30868-1 https://doi.org/10.1016/j.eswa.2019.113151 ESWA 113151
To appear in:
Expert Systems With Applications
Received date: Revised date: Accepted date:
18 August 2019 3 December 2019 19 December 2019
Please cite this article as: Wenqiang Zhang, Diji Yang, Guohui Zhang, Mitsuo Gen, Hybrid Multiobjective Evolutionary Algorithm with Fast Sampling Strategy-based Global Search and Route Sequence Difference-based Local Search for VRPTW, Expert Systems With Applications (2019), doi: https://doi.org/10.1016/j.eswa.2019.113151
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.
Highlights • HMOEA-GL combines FSS-GS and RSD-LS. • FSS-GS quickly explores the entire solution space. • RSD-LS further enhances the search ability of HMOEA-GL. • Designing suitable coding method and proper genetic operators. • Simple insertion search is used to reduce the number of vehicles.
1
Hybrid Multiobjective Evolutionary Algorithm with Fast Sampling Strategy-based Global Search and Route Sequence Difference-based Local Search for VRPTW Wenqiang Zhanga,∗, Diji Yanga , Guohui Zhangb , Mitsuo Genc a College
of Information Science and Engineering, Henan University of Technology, China of Management Engineering, Zhengzhou University of Aeronautics, China c Fuzzy Logic Systems Institute/Tokyo University of Science, Japan
b School
Abstract The vehicle routing problem with time windows (VRPTW) is an important and widely studied combinatorial optimization problems. This paper aims at VRPTW with the objectives of reducing the number of vehicles and minimizing the time-wasting during the delivery process caused by early arrival. To solve this NP-hard problem, a hybrid multiobjective evolutionary algorithm with fast sampling strategy-based global search and route sequence difference-based local search (HMOEA-GL) is proposed. Firstly, fast sampling strategy-based global search (FSS-GS) of HMOEA-GL extensively explores the entire solution space to quickly guide the search direction towards the center and edge areas of Pareto frontier. Secondly, route sequence difference-based local search (RSD-LS) is executed on the individuals with poor performance in the population obtained by FSS-GS to enhance the search ability of HMOEA-GL. In addition, the suitable coding method and proper genetic operators are designed, especially, a simple insertion search is used to reduce the number of vehicles in VRPTW. Comparing with NSGA-II, SPEA2, and MOEA/D, experimental results on 12 Solomon benchmark test problems indicate that the proposed HMOEA-GL is effective, and more excellent in convergence, while maintaining a satisfying distribution ∗ Corresponding
author Email addresses: [email protected] (Wenqiang Zhang), [email protected] (Diji Yang), [email protected] (Guohui Zhang), [email protected] (Mitsuo Gen)
Preprint submitted to Expert Systems with Applications
December 20, 2019
performance. HMOEA-GL could be an effective intelligent algorithm for expert and intelligent decision support system to help logistics companies to make decisions. Keywords: multiobjective evolutionary algorithm, vehicle routing problem, global search, local search, transportation
1. Introduction Under the environment of economic globalization and information technology, people’s online shopping and long-distance ordering is more and more convenient and popular, and logistics plays an increasingly obvious role in mod5
ern economic activities, at the same time, the logistics problem becomes more complex. It is very necessary to design and use efficient intelligent decision support system(IDSS) to help decision-makers deal with logistics problem. IDSS is a decision support system that makes extensive use of artificial intelligence (AI) techniques, which works like a human consultant to support the decision-
10
makers. In-depth analysis of the logistics problem, it can be found that the core is the vehicle routing problem (VRP). A good solution of VRP can take into account many aspects, not only can complete the customer service on time and on demand, but also can reduce the cost of distribution center consumption, to achieve a win-win. With the progress of the Internet, the type of customers
15
of VRP is more likely to be transformed from large-volume customers like production base and supermarkets to small-volume customers such as family and individuals. Unlike large-volume customers who can receive service 24 hours a day, small-volume customers are generally only allowed to receive service within a few hours. Therefore, the consideration of time window is very necessary. At
20
the same time, small-volume customers generally do not have a parking space, so the waiting cost of vehicles in the city may increase significantly, should not be ignored. With the efforts of researchers, great progress has been made in solving this problem, but due to the complexity of the problem and the differ-
3
ent research background, the research on this problem is still very important 25
and ongoing. VRP is one of the most important and widely studied combinatorial optimization problems, which is relevant to transportation logistics such as post, parcel and distribution services. Its objective is to obtain the lowest-cost set of routes. In the beginning, truck dispatching problem was proposed by Dantzig
30
and Ramser Dantzig & Ramser (1959) as a generalization of the travelling salesman problem (TSP), cost has mostly been associated with the travel distance. In 1964, Clarke and Wright summarized this problem as vehicle routing problem that is well-known: i.e., how to use a fleet of trucks to serve a set of customers which have different locations and needs. The reason VRP can be one of the
35
most widely studied topics in the field of operations research is that VRP have important practical significance and with high complexity when this problem integrated into real-life. However, the current VRP model is very different from the original, because when VRP is combined with real-life, there are many factors that cannot be ignored. Like the number of vehicles, delivery time, parking
40
cost, makespan, and workload balance, etc. are important to consider. Considering these realistic factors, VRP have several common variants that involve different constraints, like time-dependent travel times, with time windows, support delivery while pickup, can update customer information in real time, etc. Among them, the variant with time windows (VRPTW) has vehicles with lim-
45
ited capacity and the specific delivery time windows, and is particularly relevant to practical applications. The main problem studied in this paper is the vehicle routing problem with time windows which objective to minimize the number of vehicles and the waste time caused by the premature arrival. Because the limitations of many constraints of VRP, so it is often to con-
50
sider the VRP as a multiobjective problem. VRP is also considered as an NP-hard problem Lenstra & Kan (1981), when solving large-scale real-world VRPs, heuristics and meta-heuristics are often more suitable, exact algorithms are only efficient for small problem. The same is true for the variant VRPTW, optimal solutions for small instances of the VRPTW can be obtained by exact 4
55
methods, but the computation time required increases considerably for larger instances Desrochers et al. (1992). For this reason, in recent years, various heuristic and meta-heuristics algorithm for solving VRPTW have become the hotspot of research. Algorithms such as simulated annealing (SA), tabu search (TS), ant colony optimization (ACO), genetic algorithm (GA), particle swarm
60
optimization (PSO), differential evolution (DE), have proven to be effective in solving complex multiobjective problems. These algorithms are also used to solve complex VRPTW. Heuristic algorithms can greatly reduce the computational complexity and have many successful applications. However, with the increasing of the complexity and scale of the VRPTW, there may suffer some
65
problems. For example, the speed of the algorithm or the final solution set is still not satisfactory. Therefore, research on VRPTW is still ongoing. According to the No Free Lunch Theorem Wolpert et al. (1997), no algorithm can solve all the problems. Each of the classic algorithms has its own unique advantages, and inevitably there are some drawbacks. Proper combination of
70
different algorithms with different characteristics may further enhance the search effectiveness by adopting the advantages of each algorithm, and consequently may overcome the inherent limitations of single algorithm Tang & Wang (2013). Therefore, it is naturally to mix different techniques and utilize their advantages for solving complicated problems Zhou et al. (2011).
75
As we know, GA is one of the most popular type of evolutionary algorithms (EA), which can easily be adapted to various types of problems Holland (1992)Gen & Cheng (2000)Gen et al. (2008)Yu & Gen (2010)Tasan & Gen (2012). DE is a direction-based stochastic search technique which can improve the convergence of the algorithm by lead the poor solution to close to the good
80
solution Storn & Price (1997). But both of these algorithms have some drawbacks, which is easy to converge to local optima or even arbitrary points rather than the global optimum, resulting in premature convergence. This challenge is particularly hard when the dimension of problems is high and there are a lot of local optima.
85
Therefore, it is an interesting but challenging to hybridize GA and DE re5
garding the evolutionary behavior and potential advantages for VRPTW. Inspired by the successful hybrid strategies and application in previous research, in this paper, hybrid multiobjective evolutionary algorithm with fast sampling strategy-based global search and route sequence difference-based local search 90
(HMOEA-GL) is proposed to solve complicated multiobjective VRPTW. The objectives of the problems one is to reduce the number of vehicles and another is to minimize the time-wasting during the delivery process caused by early arrival. HMOEA-GL combines fast sampling strategy-based global search (FSS-GS) and route sequence difference-based local search (RSD-LS). FSS-GS is
95
designed as the global search strategy which includes center area-based elite sampling strategy and edge area-based sampling strategy. Center area-based elite sampling strategy according to Pareto dominating and dominated relationshipbased fitness function (PDDR-FF) can fast picking out the individuals located in the central area of the Pareto frontier Zhang et al. (2014). Edge area-based
100
sampling strategy based on vector evaluated genetic algorithm (VEGA) which prefers the individuals located in the edge area of the Pareto frontier Schaffer (1985). By mixing these two strategies and executing suitable genetic operators as well as simple insertion search (SIS), FSS-GS can rapidly improve the convergence and distribution performances toward the center and edge regions
105
of Pareto frontier. As a follow-up step, drawing on the evolutionary strategy of DE, RSD-LS is used to further improve the quality of individuals. RSD-LS guides individuals with poor performance to move toward better performing individuals. By hybridizing FSS-GS and RSD-LS, HMOEA-GL can achieve fast convergence and sufficient distribution. As the core of IDSS, HMOEA-GL
110
could be used to solve various scheduling problems such as VRPTW, which can provide multiple appropriate solutions for decision-makers by different weights of the objectives according to decision-makers’ needs. With HMOEA-GL, an IDSS can be established to effectively support decision-makers in a series of scheduling decision-making processes such as vehicle scheduling.
115
The rest of this paper is organized as follows. Section 2 introduces the related work to solve the vehicle routing problem with time windows. Section 6
3 provides definitions of VRPTW. A detailed description of the HMOEA-GL and the main strategies are shown in Section 4. Experimental design and the analysis of results are detailed in Section 5. Conclusions and future work are 120
provided in Section 6.
2. Related Works Over the past decades, the VRP and its variants have grown ever more popular in the academic literature. Braekers et al. present a taxonomic review of the VRP literature published between 2009 and June 2015 Braekers et al. 125
(2016). For VRPTW, Dixit et al. reviewed some of the recent advancements in the VRPTW using meta-heuristic techniques Dixit et al. (2019). SA inspired from annealing in metallurgy is often combined with other algorithms to solve the VRPTW. Pratiwi et al. proposed a natural inspired algorithms for dealing with constraints of VRPTW which involves bat algorithm and
130
cat swarm optimization. Among the proposed algorithms, bat algorithm is hybridized with SA, the worst solution of bat algorithm is replaced by the solution from SA. Algorithm which is based on behavior of cats, cat swarm optimization, is improved using crow search algorithm to make simpler and faster convergence Pratiwi et al. (2018). Gonz´ alez et al. proposed a memetic algorithm (MA) with
135
SA as trajectory-based method for solving the capacitated vehicle routing problem with time windows (CVRPTW). In MA, a novel crossover operator called the single breaking-point sequence based crossover (SBSBX) is introduced. One of the principles behind the design of SBSBX is to reduce the disruptive behavior of SBX, with the aim of providing additional intensification Gonz´ alez et al.
140
(2017). Baˇ nos, Ra´ ul et al. proposed a multiobjective procedure based on SA and the multiple temperature Pareto SA to solve VRPTW of considering the travelled distance and the imbalance of the routes Baˇ nos et al. (2013). TS is a meta-heuristic search method employing local search methods used for mathematical optimization, it was used to solve VRPTW long ago. Schnei-
145
der proposed a tabu search algorithm for the VRP with time windows and
7
driver-specific times, a variant of the classical VRPTW that uses driver-specific travel and service times to model the familiarity of the different drivers with the customers to visit Schneider (2016). Wang et al. focused on the heterogeneous multi-type fleet vehicle routing problem with time windows and an incompatible 150
loading constraint, developed a mathematical model, proposed a ruin-recreate heuristic algorithm, and a threshold tabu search method Wang et al. (2015). Yu et al. proposed a hybrid approach, which consists of TS and ACO to solve VRPTW. In proposed algorithm, a neighborhood search is introduced to improve the performance of ACO, TS is used to maintain the distribution and
155
explore new solutions when close to the convergence Yu et al. (2011). Liu and Lee proposed a mathematical model of inventory routing problem with time windows to solve the VRPTW and the inventory control decision problem simultaneously. The proposed method is divided into two phases, first is to find the initial solution and the second phase is to improve the solution adopting the
160
variable neighborhood TS Liu & Chen (2011). ACO is a probabilistic technique for solving computational problems which can be used to find good paths through graphs, has a good performance in computer science and operations research, it is also widely used in VRPTW. To solve the vehicle routing problem with service time customization (VRPTW-STC),
165
which the aim is to find an optimum solution to the smallest fleet size, the lowest travelling distance as well as the largest total service time of all customers, Wang et al. designed a multi ant system based hybrid heuristic algorithm which inspired by decomposing a multiobjective problem into several single objective ones. ACO algorithms are applied to every single-objective problem. Several
170
local search algorithms are also incorporated into MAS to improve the solution quality Wang & Xing (2018). Gupta and Saini proposed an improved ant colony optimization to solve VRPTW, which uses a new pheromone reset and update function to enhance the route searching and a 2-opt method to improve the optimized path. The experimental results verify the effectiveness of the algorithm
175
and confirms it as the better path optimizing method for VRPTW, in comparison with the traditional ACO and other heuristics Gupta & Saini (2017). Zhang 8
et al. proposed a solution strategy based on ant colony optimization and three mutation operators for multiobjective VRP with flexible time windows. In this special problem, a fleet of vehicles can service a set of customers earlier and later 180
than the required time with a given tolerance, in the hope of enabling a logistics company to save distribution costs at the expense of customer satisfaction Zhang et al. (2019). GA based on recombination, mutation and natural selection also have many related researches, which applied GA on VRPTW. G¨ oc¸ken et al. proposed a
185
hybridized version of genetic algorithm which is a meta-heuristic solution technique with constructive heuristic methods for solving VRPTW. By using sweep algorithm in initial population generation phase of GA, the search starts with high quality solution sets. Comparing with the results of GA based on nearest neighbor, the proposed GA reaches more effective solutions by beginning with
190
sweep based initial population G¨ oc¸ken et al. (2017). Abidi et al. focused on a class of problems known as rich vehicle routing problem (RVRP) with dynamically changing orders and time windows constraints. The proposed GA with a simple heuristic to solve the dynamic vehicle routing problem with time windows Abidi et al. (2018). Shi et al. proposed a hybrid GA integrated with stochas-
195
tic simulation methods to solve VRPTW with uncertain demand. Among, the uncertain demand is considered as a fuzzy variable, a fuzzy chance constraint model is designed Yong et al. (2016). Moon et al. studied VRPTW and extend two cases that allows overtime for drivers and the possibility of using outsourced vehicles. A mixed integer programming model, a hybrid algorithm based on SA
200
and GA is developed to solve this problem Moon et al. (2012). Long et al. proposed a Pareto-based evolutionary algorithm for solving prize-collecting vehicle routing problem (PCVRP). In PCVRP, besides the vehicle assignment and visiting sequencing, customer selection also must be determined. The proposed algorithm combined GA and a local search strategy, a decomposition strategy
205
is designed to decompose PCVRP into multiple sub problem Long et al. (2019). PSO is a population based stochastic optimization technique which inspired by social behavior of bird flocking or fish schooling, can also be used to solve 9
VRPTW. Marinakis et al. proposed a new variant of PSO for the solution of the VRPTW. Three different adaptive strategies are used in the proposed 210
multi-adaptive particle swarm optimization (MAPSO) algorithm. First, adaptive strategy concerns the use of a greedy randomized adaptive search procedure (GRASP). Second, new adaptive strategy concerns the adaptiveness in the movement of the particles from one solution to another. Third, the adaptiveness is in all the parameters of the particle swarm optimization Marinakis et al.
215
(2019). Zhang et al. proposed an evolutionary scatter search particle swarm optimization algorithm (ESS-PSO) to solve the VRPTW. In ESS, a genetic algorithm and a new ’route+/-’ evolutionary operator are introduced in scatter search template. Also Zhang et al. proposed a discrete PSO that sets the route-segment as the velocity of particles and in which the velocity and position
220
updating rules are designed based on the concept of ’ruin and recreate’. A new solution representation called ’auxiliary code’ based on customer allocation is designed to maintain the distribution of the reference set Zhang et al. (2018). Jiang et al. combined fields of logistics and unmanned aerial vehicle (UAV) created a model of task assignment for UAV in logistics. The model takes into
225
account multi-constraints (such as weight coefficients, time windows constraints, the constraints of the UAV and soon), and proposed an improved PSO to solve the task assignment problem with multiple constraints Jiang et al. (2017). Based on the differences between individuals, the DE can guide individuals with poor fitness to close to excellent one, which also has many research
230
works in solving the VRPTW. Pu et al. proposed a hybrid differential evolution (HDE) optimization for the vehicle routing problem with time windows and driver specific times (VRPTW-DST). In this problem, adopts service times and driver-specific travel times to model the degree of familiarity to different drivers with the service for customers. By comparing the results of experiment
235
and simulation, the knowledge for drivers in the routing problem may effectively enhance the efficiency for vehicle routesPu et al. (2017). Salamatbakhsh-Varjovi et al. developed a bi-objective mathematical model to evaluate a periodic vehicle routing problem with time windows under uncertainty for companies, and 10
an improved differential evolution algorithm was used to identify efficient solu240
tionsSalamatbakhsh-Varjovi et al. (2018). After many years of efforts by many researchers, there are already many methods to solve VRP, but still face some problems, such as the speed of the algorithm or the quality of the final solution is not satisfactory. Moreover, different researchers concern the VRP under different environmental backgrounds.
245
The influence of different environments leads to different conditional constraints, so the proposed method may only be applied to certain specific situations. In addition, the combination of different technologies to form a hybrid algorithm has become a new trend. Hybrid algorithm may have unexpected effects for solving VRP, and there are many combination methods of hybrid algorithm
250
still have not been attempted yet. Therefore, research on VRP issues is still very important and is still ongoing.
3. Problem Definitions The VRPTW is one of a variants of VRP which is NP-hard. VRPTW can be seen as a combination of the TSP and the bin packing problem (BPP). In 255
this paper, the objectives of the problem are to find set of routes with minimumwaste time and less number of vehicle to service customers. The customers have specific delivery time windows, and the vehicles are identical and have limited capacity Garcia-Najera & Bullinaria (2011). An instance of VRPTW can be defined as follows. There is a set V =
260
{0, . . . , N } of vertices, which the vertex 0 is the depot, and the other vertices
V 0 = V \{0} are defined as customers. Each customer i ∈ V 0 located at coordinates (xi , yi ), has a demand of goods gi > 0, a delivery time window [bi , ei ], bi is the beginning time of customer allow receive goods, and ei is the ending time of time windows for allow receive goods. Each customer also has a service 265
time si for unloading goods. The depot is positioned at (x0 , y0 ), has demand g0 = 0, and time window [0, e0 ≥ max{ei : i ∈ V 0 }]. The vehicles are the same and with the capacity Q ≥ max{gi : i ∈ V 0 }. 11
The distance dij and the passing time tij between the vertices i and j are generally considered to be the Euclidean distance. The problem is to find a set 270
of k routes R = {r1 , . . . , rk }, each route begins and ends at the depot, and each customer is serviced once. So each vehicle is assigned to a set of customers, and the vehicle capacity Q must be able to satisfy the sum of these customers’ demands. If rk = (u(1, k), . . . , u(Nk , k)) is the sequence of Nk customers supplied in the k-th route, where u(i, k) is the i-th customers to be visited in route k,
275
then Vk0 = u(1, k), . . . , u(Nk , k) is the set of customers serviced. There is no need to display the depot, but when calculating the fitness function, it is necessary to add the depot at the beginning and the end of the route. Dk is total distance of routing rk can be expressed as:
Dk =
Nk X
d(u(i,k),u(i+1,k))
(1)
i=0
Suppose a(u(i, k)) is the arrival time of the vehicle k at vertex i and l(u(i, k)) 280
means the time to leave the vertex. The arrival time can be calculated as follows.
a (u (i, k)) = l (u (i − 1, k)) + t(u(i−1,k),u(i,k))
(2)
Late service is not allowed, that means the route is invalid if the arrival time is later than the end of the customer’s time window. However, arriving early is allowed, but the vehicle has to wait until the beginning of the customer time window, so the vehicle has a waiting time. 0, ifa (u (i, k)) ≥ b u(i,k) w (u (i, k)) = b − a (u (i, k)) , otherwise
(3)
u(i,k)
285
The leaving time from the i-th customer in route k is described as:
l(u(i, k)) = a(u(i, k)) + w(u(i, k)) + su(i,k)
(4)
The total time of route rk can be expressed as:
Tk =
Nk X
(t(u(i,k),u(i+1,k)) + w u (i + 1, k) + su(i+1,k)
i=0
12
(5)
This paper will concentrate on two objectives, one is to reduce the number of vehicles and another is to minimize the time-consuming waste during the delivery process caused by the premature arrival. The first objective number of 290
vehicles can be expressed as:
f1 (R) = |R| = K
(6)
The second objective is total waiting time in each route k. It is necessary to mention that the waiting time of the first vertex which be serviced of each vehicle is set to zero.
f2 (R) =
Nk K X X
w (u (i + 1, k))
(7)
k=1 i=0
The relationship between customer’s demand and vehicle capacity in route 295
rk can be expressed as:
Gk =
X
0 i∈vk
gi ≤ Q, ∀k = 1, . . . , K
(8)
The arrival times need to meet customer’s time window constraints.
a (u (i, k)) ≤ eu(i,k) , ∀k = 1, . . . , K, 1 ≤ i ≤ Nk
(9)
4. Algorithm Description As the VRPTW belongs to NP-hard problem with the high computational complexity and lots of conditional constraints, in this paper HMOEA-GL is pro300
posed for solving VRPTW. In this section, section 4.1 introduces the overview of HMOEA-GL. Section 4.2 describes the details of FSS-GS. Section 4.3 shows the genetic representation of HMOEA-GL, including encoding and decoding. The genetic operators of HMOEA-GL will be described in section 4.4. Section 4.5 introduces the process of SIS and RSD-LS will also be introduced in section
305
4.6.
13
4.1. Overview of the HMOEA-GL The main advantages of HMOEA-GL are hybridizing FSS-GS and RSD-LS, FSS-GS has fast convergence and sufficient distribution, and subsequently using RSD-LS to further improve the performance of the algorithm. Among FSS-GS, 310
center area-based elite sampling strategy according to a fitness function called PDDR-FF Zhang et al. (2014). PDDR-FF is used to evaluate the individuals that can fast determine which individuals are nondominated or dominated but have large domination area. Combining with edge area-based sampling strategy based on VEGA Schaffer (1985), FSS-GS increases the attention to individuals
315
which located on the edge and center regions of the Pareto front. Therefore, FSS-GS can effectively improve the convergence performance of HMOEA-GL while ensuring an acceptable distribution performance. FSS-GS is designed as the global search strategy of HMOEA-GL. In addition, RSD-LS is a local search which can guide the individuals with poor performance move closer to
320
those who perform better. By this way, the search ability of HMOEA-GL can be further enhanced. In addition, in order to reduce the number of vehicles, a single objective optimization method SIS is used in HMOEA-GL, and other suitable genetic operators are also used to recombine chromosomes. HMOEAGL can be divided into four steps roughly, and the framework of HMOEA-GL
325
is shown as follows, and the framework of HMOEA-GL is shown in figure 1, and the pseudo code is shown in algorithm 1. Step 1: Using the edge area-based sampling strategy, individuals are selected from the population Pt by binary tournament selection. Individuals who perform well on objective 1 will be put into SP1 . Similarly, the individuals
330
which perform well on objective 2 will be placed in SP2 . Step 2: Mixing SP1 , SP2 and elite population At as a mating pool, and pick out the parent individuals from it. Performing crossover operator, and mutation operator on parent individuals to generate the offspring. Executing 0 . SIS on offsprings to form temporary population Pt+1
335
Step 3: Performing RSD-LS operations on individuals in the temporary 0 population Pt+1 . Based on the difference between the two individuals, the
14
Start
Initializing Population Pt and Elite Population At
Yes
Stop
Output Pareto solution set
No FSS-GS: edge area-based sampling strategy
Step1
Creating Sub-Population SP1 and SP2 according to VEGA 1. Selecting individuals in Pt who perform well on objective 1 into SP1 2. Selecting individuals in Pt who perform well on objective 2 into SP2
Mixing SP1, SP2 and At as Mating pool Genetic operators: Selection Crossover Mutation
Step2
Simple Insertion Search
Creating temp. Population P’ t+1 1. Creating offSpringSet by perform genetic operators on Mating pool 2. Performing simple insertion search on offSpringSet
Route Sequence Difference-based Local Search
FSS-GS: center area-based elite sampling strategy
Step3
Creating Population Pt+1 1. Creating RSD-LS Population RSD-LSP’ t+1 by performing RSD-LS on P’ t+1 2. Selecting individual from RSD-LSP’ t+1 and P’ t+1 according to PDDR-FF
Step4
Updating Elite Population At to Elite Population At+1 1. Uniting Pt+1 and At as tempEliteSolutionset 2. Selecting individual from tempEliteSolutionset according to PDDR-FF
FSS-GS: center area-based elite sampling strategy
Figure 1: Basic framework of the HMOEA-GL.
15
Algorithm 1 HMOEA-GL Input: data set of VRPTW (for each instance); EA parameters: PopSize, MaxGenerationsSize, ElitePopSize, subPopSize; Output: Pareto SolutionSet; 1: t = 0; 2: Initializing population Pt by encoding routine; 3: Calculating objectives f1 , f2 by decoding routine; 4: Calculating fitness eval by PDDR-FF routine; 5: Initializing elite population At by center area-based elite sampling strategy; 6: while (t ≤ MaxGenerationsSize) do 7:
Creating subPopulation1 SP1 by edge area-based sampling strategy;
8:
Creating subPopulation2 SP2 by edge area-based sampling strategy;
9:
Mixing SP1 , SP2 , and At as Mating pool;
10:
Creating offspringSet by genetic operators on Mating pool;
11:
Performing simple insertion search on offspringSet;
12:
0 Temp population Pt+1 = offspringSet;
13:
0 0 ; Performing RSD-LS on Pt+1 to create RSD-LS Pt+1
14:
0 0 and RSD-LS Pt+1 Creating temp. solutionSet tempSet union of Pt+1 ;
15:
Creating Pt+1 from tempSet by center area-based elite sampling strategy;
16:
Calculating objectives f1 , f2 by decoding routine;
17:
Calculating fitness eval by PDDR-FF routine;
18:
Creating a temp. elite solutionSet tempEliteSet union of Pt+1 and At ;
19:
Creating At+1 from tempEliteSet by center area-based elite sampling strategy;
20:
t = t + 1;
21: end while 22: Creating Pareto SolutionSet by selecting nondominated in At ;
16
individual with poor fitness function can obtain a guiding evolution. PDDR-FF will be used to evaluate the new individual generated by RSD-LS, which judges the new individual whether can be accepted or not. The population optimized 340
by RSD-LS is determined as the new population Pt+1 . Step 4: Using the center area-based elite sampling strategy to update the elite population. The individuals both in the new population Pt+1 and elite population At are evaluated. Best | A | individuals will be selected and formed the new elite population At+1 .
345
HMOEA-GL can be divided into two parts. One is to use the FSS-GS to fast improve the convergence and distribution performance of the HMOEA-GL, and another is to use RSD-LS to guide individuals with poor fitness to further evolve, thus the performance of HMOEA-GL can be improved. The specific details of FSS-GS, RSD-LS, SIS, coding method and various genetic operators
350
in HMOEA-GL will be described later. 4.2. Fast Sampling Strategy-based Global Search FSS-GS can fast converge toward the center and edge regions of Pareto frontier while ensuring sufficient distribution by combining the center area-based elite sampling strategy and the edge area-based sampling strategy. Among,
355
a fitness function called PDDR-FF plays a key role in center area-based elite sampling strategy Zhang et al. (2014). PDDR-FF is used to evaluate individuals according to the dominant relationship among individuals, and each individual can obtain an evaluation value. With the evaluation value of PDDR-FF, FSSGS can fast identify which individuals are nondominated. Moreover, individuals
360
who are dominated but have a larger domination area (have great potential to become nondominated) are also evaluated as good. The specific formula of PDDR-FF for an example of calculating the fitness function of individual Si will be shown in following.
eval (Si ) = q (Si ) +
1 , i = 1, 2, . . . , popSize p (Si ) + 1
17
(10)
where q(Si ) is the number of the individuals that dominate Si , p(Si ) is the 365
number of the individuals who are dominated by Si , and popSize is the size of population. Obviously, the smaller the eval, means the individual is better. The ability of converging to the center area of Pareto frontier can be effectively improved by the center area-based elite sampling strategy. At the same time, in order to maintain the distribution performance of the algorithm,
370
FSS-GS combines center area-based elite sampling strategy with the edge areabased sampling strategy, which is according to the sampling strategy of VEGA Schaffer (1985). The sampling strategy of VEGA can be simply described as decomposing multiobjective problem into multiple single-objective problems, then set a sub-population for each objective. Each sub-population contains individ-
375
uals who prefer one certain objective. The edge area-based sampling strategy prefers to the individuals located in the edge area of Pareto frontier, that can make up for the deficiency of distribution performance of center area-based elite sampling strategy effectively. Taking a bi-objective minimization problem as an example, the main mechanism of FSS-GS is shown in figure 2. min f2
Edge area-based sampling strategy
Individuals prefer objective 1
Center area-based elite sampling strategy
Individuals with a good eval Individuals prefer objective 2 min f1
Figure 2: The main mechanism of Fast Sampling Strategy-based Global Search.
380
According to the center area-based elite sampling strategy, FSS-GS can fast pick out the individuals located around the central area on the Pareto frontier,
18
and the edge area-based sampling strategy prefers the edge area on the Pareto frontier. Mixing two sampling strategies, FSS-GS can effectively improve the convergence and distribution performance of the HMOEA-GL. 385
4.3. Genetic Representation 4.3.1. Encoding and Initialization The encoding of HMOEA-GL is simple and intuitive, which uses permutation code. The length of the chromosome is equal to the number of the customer vertices. Each gene of chromosome represents the one of the customer number.
390
In one chromosome, one customer number will only appear once and recurring is not allowed. In order to easily apply genetic operators on chromosome, there is no depot in the chromosome. But the depot needs to be taken into account when decoding and calculating fitness. When initializing, the population is formed of chromosomes which are randomly generated.
395
4.3.2. Decoding When decoding, the genes of chromosome will be traversed by the order from the first one to the last. Cutting the chromosome when exceeded the conditional limit, then the chromosome is separated into multiple sub-routes. One vehicle is designated to serve customers according to one of the sub-route. It should
400
be noted that each sub-route needs to add the depot at the beginning and the end. An example of decoding is shown in figure 3, and pseudo code of decoding in HMOEA-GL is also be shown in algorithm 2. Taking figure 3 as an example, it can suppose that there is a depot 0, and a group of customers numbers is 1-6. Each customer has a time window indicated
405
by square brackets and the distance between the customers determines the travel time of the vehicle. The service time for each customer is set to 1. For customer 1, the begin time of receiving goods is 1, and the end time is 6. The distance between the depot 0 and customer 1 is set to 1, then distance between the depot 0 and customer 6 is approximately considered to be 1.3. A chromosome {1, 5,
410
4, 2, 3, 6} will be decoded into two sub-routes {0, 1, 5, 4, 2, 0} and {0, 3, 6, 19
[5, 10] [3, 7]
4
1
0.7 [9, 15]
5 Chromosome:
1
5
4
2
3
6
2
1 [2, 6]
Vertex sequence:
1
5
4
2
3
1
6
1
0
1
exceeded conditional limit
Sub-route:
0
1
5
4
Vehicle_1
2
0
0
3
6
0
1
3
1.3
Vehicle_2
[2, 7] 1.5
6 [2, 5]
Figure 3: Decoding method in HMOEA-GL.
0}. When the vehicle finishes serving the customer 2, the end time is 10, it has no enough time to complete the service requirements of customer 3, so the vehicle is returned to the depot and use another vehicle to serve the customer 3. Conditional limit is not only limited by time, but also the load of the vehicle. 415
4.4. Genetic Operators 4.4.1. Selection Operator In HMOEA-GL, the binary tournament selection is used as the selection operator. Since the purpose of using the edge area-based sampling strategy based on VEGA is to improve the distribution of the algorithm, so in the process
420
of forming sub-populations, it is not necessary to divide all the individuals which have excellent performance on one particular objective into one sub-population. Instead, select two individuals randomly and choose the one who performs well on particular objective, and put it into the corresponding sub-population. The center area-based elite sampling strategy according to PDDR-FF is used to
425
improve the convergence of the algorithm. So is often to order the individuals of the population by their fitness values, then picks out the best part of them as the elite population.
20
Algorithm 2 Decoding procedure of Chromosome Input: VerticesData x, y, b, e, s, g; Distance d(i, j)=Euclidean(vi , vj ); VehicleData velocity V =1, capacity Q; Chromosome Sequence C; Output: Routes R = r1 , . . . , rk ; 1: Vn =0; // Vn is the vehicle number 2: cg =0; // cg is the current load of vehicle 3: at = bC(0) ; // at is the arrival time of vertex t 4: lt = bC(0) + sC(0) ; //lt is the leave time of vertex t 5: bdt ; // bdt is the time from vertex t back to depot 0 6: add v0 into r0 ; // add the depot as the beginning of route 7: add C(0) into r0 ; // add the first vertex in C into first route r0 8: cg = cg + gC(0) ; 9: for i = 1 to length of C do 10:
at = lt + d(C(i − 1), C(i))/V ;
11:
cg = cg + gC(i) ;
12:
bd1t = bC(i) + sC(i) + d(C(i), 0)/V ; //for early arrival
13:
bd2t = at + sC(i) + d(C(i), 0)/V ; //for no early arrival
14:
if at < bC(i) and cg ≤ Q and bdt ≤ e0 then
15: 16: 17:
add C(i) into rV n ; //C(i) belongs to the current vehicle lt = bC(i) + sC(i) ; else if bC(i) ≤ at ≤ eC(i) and cg ≤ Q and bd2t ≤ e0 then
18:
add C(i) into rV n ; //C(i) belongs to the current vehicle
19:
lt = at + sC(i) ;
20: 21:
else if at > eC(i) or cg > Q or bd1t > e0 or bd2t > e0 then //C(i) isn’t served by the current vehicle
22:
add v0 into rV n ; //add depot as the ending of route
23:
Vn = Vn + 1; //C(i) is served by the new vehicle
24:
add v0 into rV n ; //add depot as the beginning of route
25:
add C(i) into rV n ;
26:
lt = bC(i) + sC(i) ;
27: 28:
cg = gC(i) ; end if
29: end for
21
4.4.2. Crossover Operator Since the permutation code is adopted, HMOEA-GL uses order crossover 430
(OX) as the crossover operator to perform genetic recombination process Davis (1991). An example of order crossover for permutation code is shown in figure 4. Parent 1 Offspring 1
1 2 3 4 5 6 7 8 9
2 1 3 4 5 6 7 8 9 5 4 6 9 2 1 7 8 3 Parent 2
Figure 4: An example of order crossover.
4.4.3. Mutation Operator Exchange mutation (EM) is the mutation operator of HMOEA-GL Banzhaf 435
(1990), which aims to increase the search ability of the algorithm by causing mutations of the individual. In order to increase the intensity of variation, this paper turns the exchange of gene position from the single into fragments. An example of EM for permutation code is shown in figure 5. 1 2 3 4 5 6 7 8 9
1 7 8 4 5 6 2 3 9
Figure 5: An example of exchange mutation.
4.5. Simple Insertion Search 440
As we all know, the number of vehicles is one of the most important objectives in VRP. Whether in the initial VRP or in variants of various VRP problems, it is necessary to reduce the number of vehicles as much as possible. In reality, adding one vehicle is very likely to make other objectives getting worse such as a significant increase in cost, and it is meaningless that a solution
445
which requires a large number of vehicles. Therefore, this paper uses SIS to 22
reduce the number of vehicles used in one solution. The use of SIS for individuals generated by genetic operation can be seen as a secondary mutation of the individual. Figure 6 shown an example of the main mechanism of SIS. [5, 10] [3, 7]
[3, 7]
5 1.3
[9, 15]
1.5 2
0.7 [9, 15]
2
1 [2, 6]
1
1
1
0
1
4
1
5
2
[2, 6]
1
[5, 10]
4
1
1
3
1.3
0
1
1
3
1.3
[2, 7] 1.5
[2, 7] 1.5
6
6
[2, 5]
[2, 5]
Figure 6: An example of the main mechanism of simple insertion search.
Among them, 0 is the depot, and 1 to 6 are six customers. Each customer 450
has a time window indicated by square brackets, and needs 1 service time. The distance between the customers determines the travel time of the vehicle. Taking figure 6 as an example, figure 7 shows the detailed process of SIS, and the pseudo code is shown in algorithm 3. SIS is used to insert the next customer into a certain route of the current
455
vehicle. Chromosome {1, 2, 4, 5, 3, 6}, which direct decoding will result in a solution that requires 3 vehicles. SIS judges whether the customer 2 can be served before the customer 1, and does not affect the service request of customer 1. Assume that the vehicle first serves customer 2, the end time is 10, and takes 2 units of time to go to customer 1. When the vehicle arriving at customer
460
1, the time is already 12, which has exceeded the time window of customer 1. Obviously, before the customer 2 is placed in the customer 1, the service of the customer 1 will be abnormal, so give up inserting the customer 2 between the depot and the customer 1. Customer 1 is already the last service customer in the current sub-route, then judges if customer 2 can be served after customer
465
1. Similarly, judges in turn whether the customer 4 can be inserted before the customer 1 or before the customer 2. It can be seen when the vehicle completes 23
Chromosome:
1
2
4
5
3
6
1
2
4
5
3
6
0
0
4
5
0
0
Direct decoding
0
1
2
Vehicle_2
Vehicle_1
3
0
6
Vehicle_3
Simple insertion search
x 0
2
1
0
√
1
x 0
1
2
0
1
2
0
1
4
√
2
0
4
1
2
4
2
x
0
1
0
1
4
2 √
5
4
2
x x 3
5
0
0
1
5
4
1
5
4
2
3
6
2
0
0
3
.. . 0
1
5
4
2
0
0
Vehicle_1
0
Vehicle_2
Figure 7: An example of the detailed process of Simple Insertion Search.
servicing customer 1, the end time is 2, after that, takes 2 units of time to go to customer 4, and spends 1 unit time waiting for customers 4 to start receiving the goods, the end time after the service for the customer 4 is 6. At this time, 470
take 0.7 units of time to go to customer 2, and it can be found that the service of the customer 2 is not affected. Therefore, it is feasible to insert the operation of the customer 4 between the customer 1 and the customer 2. In the same way, insert customer 5 before customer 4. In the end, the solution obtained through a simple insertion search requires only 2 vehicles and the number of vehicles has
475
been reduced. 4.6. Route Sequence Difference-based Local Search Differential evolution is an efficient evolutionary algorithm that guides the evolution of individuals based on differences among individuals. HMOEA-GL 24
Algorithm 3 Simple Insertion Search Input: VerticesData x, y, b, e, s, g; Distance d(i, j)=Euclidean(vi , vj ); VehicleData velocity V =1, capacity Q; Chromosome Sequence C; Output: Chromosome after Simple Insertion Search simC;
1: cg =0; // cg is current route load 2: aC(0) = bC(0) ; // at is arrival time of vertex t 3: lC(0) = bC(0) + sC(0) ; //lt is leave time of vertex t 4: 5: 6: 7: 8: 9: 10: 11: 12: 13: 14: 15: 16: 17: 18: 19: 20: 21: 22: 23: 24: 25: 26: 27: 28: 29: 30: 31: 32: 33: 34: 35: 36: 37: 38:
Vn = 0; // Vn is the vehicle number add v0 into r; // add the depot as the beginning of route Vn = Vn + 1;
rVn = v0 ;
add C(0) into r; // add the first vertex in C into first route r Vn = Vn + 1;
rVn = C(0);
cg = cg + gC(0) ;
for i = 1 to length of C do f lag = 0; //flag to judge insert or not,1 is yes,0 is no for j = 1 to Vn − 1 do ctemp = cg + gC(i) ; g
atemp = lrj + d(rj , C(i))/V t
if atemp ≤ bC(i) and ctemp ≤ Q then g t lttemp = bC(i) + sC(i) ;
else if bC(i) < atemp ≤ eC(i) and ctemp ≤ Q then g t lttemp = atemp + sC(i) ; t
else if atemp > eC(i) and ctemp ≤ Q then g t continue; else if ctemp > Q then //This vertex cannot be serviced in r g break; end if atemp2 = lttemp + d(C(i), rj+1 )/V ; t if atemp2 ≤ brj+1 then t Vn = Vn + 1 add C(i) into r between rj and rj+1 ; cg = cg + gC(i) ;
lC(i) = lttemp ;
flag = 1;
break; else if atemp2 > brj+1 then t continue; end if end for if flag = 1 then
continue;
else if flag = 0 then //Determine if C(i) can be serviced at the end of the r if C(i) can be serviced in r then add C(0) into the end of r; else if C(i) can not be serviced in r then record the sequence of r into simC;
39: use a new vehicle to serve customer C(i), 40: end if 41: end if 42: end for
25
borrowed the idea of DE to further improve the performance of the algorithm. 480
With appropriate evolutionary operations, FSS-GS can find effective solutions fast. To further improve the quality of solutions, HMOEA-GL uses a local search called route sequence difference-based local search. RSD-LS is an additional search operation on the population generated by FSS-GS during the iteration of HMOEA-GL. The main idea of RSD-LS can be roughly divided into four
485
steps. • Step 1: Randomly pick an individual from the population. • Step 2: Select another individual. If it is better or worse than the individual selected in Step 1, then execute Step 3, otherwise repeat Step 2.
490
• Step 3: Use exchange order to determine the difference of route sequence between individuals selected in Step 1 and Step 2. • Step 4: Taking a certain proportion of the difference sequence obtained by Step 3, then apply it to the individual who performs poorly. Edge area-based sampling strategy
min f2
Center area-based elite sampling strategy
DLS
(S1 – S2)
S2
F(S1 – S2) S2 DLS
S 2` S 2`
S1
DLS S1 DLS DLS min f1
Figure 8: An example of Route Sequence Difference-based Local Search.
26
Taking figure 8 as an example, according to the above expression, RSD-LS 495
process can be approximated by a formula. S20 = S2 + F (S1 − S2 )
(11)
where S1 and S2 are two individual selected by Step 1 and Step 2. Suppose S2 is the poor one, that is, S2 is dominated by S1 . (S1 -S2 ) is the difference of route sequence between individuals S1 and S2 . F is the proportion, and S20 is the new individual that evolved from S2 by RSD-LS. Each new individual 500
will be evaluated by PDDR-FF, and new individuals which superior to previous evolution will be accepted. Otherwise the original individual has to put into the next generation of the population. By using RSD-LS, the population can be closer to the Pareto front. Because HMOEA-GL uses permutation code, so HMOEA-GL uses exchange
505
sequence to determine the difference of route sequence between individuals. Taking figure 9 as an example to describe the process of how to use exchange order to define the differences of route sequence between two individuals for permutation code. S1 1 2 3 4 5 6 S2 3 2 4 5 1 6
S1
1 2 3 4 5 6
S 2'
3 2 4 5 1 6
S1
1 2 3 4 5 6
S 2'
1 2 4 5 3 6
S1
1 2 3 4 5 6
S 2'
1 2 3 5 4 6
S 2'
1 2 3 4 5 6
{(0,4)}
{(0,4),(2,4)}
{(0,4),(2,4),(3,4)} Over
Figure 9: An example of Exchange Order.
27
Suppose S1 = {1, 2, 3, 4, 5, 6}, S2 = {3, 2, 4, 5, 1, 6}, and calculating process 510
of the difference of route sequence between S1 and S2 is as follows. Starting from the first gene position of the chromosome, compare each gene of the two individuals. Mark the first position subscript (0, x) when there has a difference in S2 . Then find the target gene 1 in S2 , and mark the second position subscript (0, 4). Combining two position subscripts becomes an exchange sequence
515
{(0, 4)}. Subsequently, exchange genes in S2 to obtain a temporary chromosome
S20 . Continue to compare S20 and S1 until the last set of genes of the chromosome
is detected. Sequence {(0, 4), (2, 4), (3, 4)} is the exchange order belong to S2 , represents the difference of route sequence between S1 and S2 . The max length of the exchange sequence for S2 is 3, which can be processed by proportion F . 520
Proportion F is a random number which F ∈ [0, 1]. Suppose F = 0.5, then the length of the exchange sequence executed is 1 (multiply max length of exchange by proportion F and round down). That is, the actual execution sequence is {(0, 4)}, the new individual S20 is generated as S20 = {1, 2, 4, 5, 3, 6}. Based on the route sequence differences between individuals, RSD-LS guides
525
the individuals with poor performance close to the better one, and the search ability of the algorithm can be further improved. 5. Experimental Evaluation In this part, section 5.1 describes the overall setup of the experiment. Section 5.2 introduces the evaluation indicators used in this experiment. Section
530
5.3 discusses the parameter setting of the algorithms. Section 5.4 verifies the effectiveness of SIS. Section 5.5 shows the experimental results and some conclusions. 5.1. Experimental Setup To verify the performance of proposed HMOEA-GL, a publicly available
535
benchmark generated by Solomon is used in this experimental studies Solomon (1987). The problem number used in this experiment are R101-R104, C101C104, RC101-RC104, which is random and with short scheduling horizon. Each 28
benchmark gives information about the depot and 100 customers, including coordinate locations, demand, time window, and service time. In this experiment, 540
indices C, GD, IGD, HV , Spacing and Spread are used to evaluate the performance of the algorithms. Wilcoxon rank-sum test is used for statistical test between proposed algorithm and competing algorithms. In order to make the comparison of experiments more fairer, the parameters of each algorithm are adjusted and optimized before the final experiment. In addition, since the SIS
545
is used for each algorithm to reduce the number of vehicles in this experiment, the effect of the SIS is also verified. Experiments are conducted to compare the effects of the proposed algorithm, NSGA-II, SPEA2 and MOEA/D on 12 Solomon benchmark problems, each algorithm is executed 30 times and record the solutions in the Pareto approximation each time. The population size is 50
550
and the max number of generations is set to 500. Selection operator uses binary tournament selection, crossover operator uses order crossover, mutation operator is exchange mutation, the SIS applied to all algorithms. The experimental results are marked, the darker the cell background color, the better the result. 5.2. Performance Indices
555
Coverage (C). C(S1 , S2 ) is the percent of the individuals in S2 who are weakly dominated by S1 . The larger C(S1 , S2 ) is, the better S1 outperforms S2 in C. Generational Distance (GD). GD evaluates a solution set according to the shortest distance of each solution from the P F ∗ . The smaller GD for one
560
solution set is, the better it is in approaching P F ∗ . In this experiment, the P F ∗ of the problem consists of the best results that all algorithms run multiple times. Inverted Generational Distance (IGD). IGD calculates the distance from P F ∗ to the solution set. IGD mainly evaluates the convergence ability, but
565
when the solution set is approaching P F ∗ , smaller IGD requires the individuals being distributed similarly to P F ∗ . The smaller IGD might mean better convergence and distribution. In IGD, the P F ∗ of the problem is also consists 29
of the best results that all algorithms run multiple times. Hypervolume (HV ). HV is the size of the dominated space of the obtained 570
solution set. Larger HV means better convergence performance. The reference point of HV in this experiment is set as combining the worst fitness values in each objective, which found by all algorithms run multiple times. Spacing (SP ). SP is a distribution performance evaluation indicator which according to the closest distance between each solution and other solutions. A
575
smaller SP might mean better distribution. Spread (Spread). Spread calculates the average of the distance between consecutive solutions. At the same time, consider the smallest distance between edge points and boundary solutions. A smaller Spread means a better distribution and spread.
580
5.3. Parameter Settings This part mainly compares the influence of different crossover rate Cr , mutation rate Mr , and proportion F of RSD-LS on the HMOEA-GL. Table 1 shows the performance indices of the final solution set obtained by running different parameter HMOEA-GL under the Solomon benchmark R101 problem. Based
585
on previous experience, the crossover rate Cr is generally set large, and the mutation rate Mr is generally small. If the proportion F of RSD-LS is too large, it will easily produce the exact same solution, and resulting to fall into local optimum. So the value of Cr , Mr and F is set to: Cr ∈ {0.5, 0.6, 0.7, 0.8, 0.9}, Mr ∈ {0.1, 0.2, 0.3}, and RSD-LS proportion F ∈ {0.5, 0.6, 0.7, 0.8}.
590
It can be seen that the better performing combinations are mainly concentrated between experiments 25-28 and experiments 37-40. In experiments 25-28 the crossover rate Cr is 0.7, and the mutation rate Mr is 0.1. The crossover rate Cr of experiments 37-40 is 0.8, the mutation rate Mr is 0.1. Among them, from the number of times of get top three of evaluation indicator, experiments 37-40
595
is more than experiments 25-28. In experiments 25-28 and experiments 37-40, it can see that when the proportion F of RSD-LS is set to 0.8, HMOEA-GL can perform better. Therefore, the parameters of HMOEA-GL are set as follows: 30
Table 1: The result of HMOEA-GL in different parameter in R101 problem. Problem: R101 Experiment 1 Experiment 2 Experiment 3 Experiment 4 Experiment 5 Experiment 6 Experiment 7 Experiment 8 Experiment 9 Experiment 10 Experiment 11 Experiment 12 Experiment 13 Experiment 14 Experiment 15 Experiment 16 Experiment 17 Experiment 18 Experiment 19 Experiment 20 Experiment 21 Experiment 22 Experiment 23 Experiment 24 Experiment 25 Experiment 26 Experiment 27 Experiment 28 Experiment 29 Experiment 30 Experiment 31 Experiment 32 Experiment 33 Experiment 34 Experiment 35 Experiment 36 Experiment 37 Experiment 38 Experiment 39 Experiment 40 Experiment 41 Experiment 42 Experiment 43 Experiment 44 Experiment 45 Experiment 46 Experiment 47 Experiment 48 Experiment 49 Experiment 50 Experiment 51 Experiment 52 Experiment 53 Experiment 54 Experiment 55 Experiment 56 Experiment 57 Experiment 58 Experiment 59 Experiment 60
Cr
Mr 0.1
0.5
0.2
0.3
0.1
0.6
0.2
0.3
0.1
0.7
0.2
0.3
0.1
0.8
0.2
0.3
0.1
0.9
0.2
0.3
F 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8 0.5 0.6 0.7 0.8
GD 21.3388 23.8032 16.6762 20.0018 25.4418 23.7405 29.8896 23.0882 30.9756 46.1087 36.6538 33.4476 16.3872 16.8135 18.4615 20.8670 29.9787 23.9819 19.7299 23.5509 29.0993 28.0303 30.8049 30.9219 12.9247 14.1581 13.9804 14.2389 16.4505 20.6195 19.6137 20.1798 29.3523 20.9034 25.5161 26.9648 12.1302 19.7369 12.9911 14.5651 16.7877 16.8810 23.2771 18.0375 22.3219 23.4938 26.1277 25.1435 15.6213 16.1051 16.3387 15.6077 21.6154 20.6990 18.8295 15.9019 22.0410 30.1262 19.4764 24.5973
31
IGD 13.3350 13.7668 12.8990 14.6975 16.1033 15.1378 15.7111 16.7063 19.0228 19.9408 19.2888 17.5324 11.0455 10.9044 11.7491 12.2983 13.8058 14.1500 12.0492 13.5544 16.4089 15.7600 16.9685 17.3720 9.7634 9.9506 10.5774 9.0088 12.5619 12.3818 12.0495 12.9600 17.5648 15.3338 16.1236 14.0746 9.9663 10.2646 8.9186 9.3273 11.3039 12.0827 12.3263 11.3698 14.6250 15.6468 15.0392 15.1042 10.8904 9.8199 10.1662 10.2055 13.2986 12.9988 12.1309 12.4055 15.5344 16.6866 14.8417 16.0166
HV 19887.0763 20119.8256 20162.2998 19898.8391 19658.3691 19838.2241 19880.7286 19988.1011 19254.2900 19086.5119 19414.6001 19657.7081 20471.0208 20438.4427 20623.4639 20461.0787 20169.9439 20300.1074 20568.4782 20205.0562 19662.4597 19937.2317 19589.2903 19822.9547 20626.8089 20773.6158 20637.2858 21121.5829 20367.7308 20265.9197 20582.3587 20539.8800 19733.4465 20125.2448 19838.5651 20196.2071 20758.1035 20726.7729 20928.8921 21046.1259 20503.8451 20537.7070 20433.4750 20580.5499 20228.5421 19907.7994 20195.4369 20370.9136 20834.4784 20818.3227 20704.2913 21045.9968 20397.1851 20390.4750 20397.9926 20550.8149 19961.9319 19946.6018 20063.2120 19966.5271
Spacing 17.9338 19.3704 17.1619 18.7967 23.1657 19.4627 25.3967 16.8518 17.9079 34.9161 31.6266 31.1964 17.3784 18.4281 22.7160 17.6778 24.8959 17.3863 18.9225 18.1104 20.8584 22.6456 25.3131 24.0197 17.8313 14.6415 13.1678 21.6098 12.5120 18.8373 19.8368 17.8454 22.2331 17.4680 21.9535 26.0599 13.7202 24.9500 15.8211 16.2298 20.1938 16.5157 28.3139 17.3354 18.3644 13.0944 19.4031 21.9715 16.5970 20.0627 17.6052 20.0819 19.2909 18.6517 24.8851 14.6071 21.0580 23.0481 19.8246 22.2405
Spread 0.9233 0.9530 0.9434 0.9659 0.9802 0.9228 1.0051 0.9214 0.9589 1.0237 1.0553 1.0098 0.9672 0.9902 1.0284 0.9682 0.9934 0.9368 0.9630 0.9191 0.9543 0.9839 0.9893 1.0090 1.0008 0.9678 0.9599 0.9994 0.9471 0.9599 0.9531 0.9441 0.9623 0.9326 0.9525 1.0178 0.9659 0.9946 0.9666 0.9884 1.0280 0.9791 1.0006 0.9751 0.9876 0.9165 0.9588 0.9461 0.9838 0.9753 0.9754 1.0121 0.9761 0.9818 0.9969 0.9025 0.9306 0.9196 0.9389 0.9628
crossover rate Cr =0.8, mutation rate Mr =0.1, RSD-LS proportion F =0.8. Table 2: The result of NSGA-II in different parameter in R101 problem. Problem: R101
Cr
Mr
GD
IGD
HV
Spacing
Spread
0.1
45.8989
24.9881
16632.5620
28.0686
0.9587
0.2
49.4310
28.7195
16619.2598
25.9404
0.9556
Experiment 3
0.3
54.5880
28.5689
16229.2180
24.7489
0.9457
Experiment 4
0.1
46.9625
28.5377
16873.7573
20.3073
0.9208
0.2
55.2580
30.1363
16406.4515
20.9027
0.9321
0.3
74.7008
34.9182
15672.5058
27.7712
0.9915
0.1
69.6310
37.5529
16079.3407
28.8825
0.9664
0.2
71.8652
39.6551
15995.1530
23.8254
0.9201
0.3
91.9884
47.0863
15185.3113
17.9976
0.8885
Experiment 1 Experiment 2
Experiment 5
0.7
0.8
Experiment 6 Experiment 7 Experiment 8
0.9
Experiment 9
Table 3: The result of SPEA2 in different parameter in R101 problem. Problem: R101
Cr
Mr
GD
IGD
HV
Spacing
Spread
0.1
22.3897
13.5210
18759.3233
14.4591
0.9351
0.2
29.4389
16.8471
18268.8967
17.8439
0.9539
Experiment 3
0.3
41.8934
18.5128
18010.6933
22.3111
0.9996
Experiment 4
0.1
27.0167
14.6118
18541.1259
11.2576
0.8937 0.9195
Experiment 1 Experiment 2
Experiment 5
0.7
0.2
27.6127
16.4009
18586.4625
15.7125
Experiment 6
0.3
46.3519
19.9279
17565.6125
14.3577
0.9243
Experiment 7
0.1
32.4969
15.4946
18548.8857
14.3954
0.9058
0.2
36.3992
16.8760
18342.2598
17.4023
0.9537
0.3
42.5555
17.7649
17957.3059
12.8516
0.9428
Experiment 8 Experiment 9
0.8
0.9
Similarly, for the competing algorithms, table 2 to 4 shows the performance 600
of NSGA-II, SPEA2 and MOEA/D with different parameters under Solomon benchmark R101 problem. It can be seen that for NSGA-II and SPEA2, the better parameters are also set to crossover rate Cr =0.8 and mutation rate Mr =0.1. For MOEA/D, can be found that MOEA/D performs better in the case of crossover rate Cr =0.9 and mutation rate Mr =0.1. In order to make a more fair
605
comparison, MOEA/D under two parameters MOEA/D-a(Cr =0.8, Mr =0.1) and MOEA/D-b(Cr =0.9, Mr =0.1) is set in the following experiments. 5.4. Effect of the SIS Component Section 4.5 refers to a local search method called SIS that can reduce vehicle number. In this part, experiments are conducted between HMOEA-GL with SIS 32
Table 4: The result of MOEA/D in different parameter in R101 problem. Problem: R101
Cr
Mr
GD
IGD
HV
Spacing
Spread
0.1
15.4365
15.2818
14697.7947
14.0840
0.9604
0.2
13.4091
14.2141
15218.6083
15.0232
0.9909
0.3
11.6031
12.6066
16047.6508
12.0886
0.9466
0.1
12.6256
14.8540
15002.5816
12.3173
0.9770
0.2
10.4925
12.1173
16308.3994
10.6865
0.9208
Experiment 6
0.3
9.5854
11.9201
16612.0227
11.3378
0.9655
Experiment 7
0.1
7.8974
12.9815
16285.5243
5.8513
0.9284
0.2
9.1806
12.9264
16228.4371
10.3583
0.9744
0.3
9.3730
12.1153
16516.0961
11.8013
0.9653
Experiment 1 Experiment 2
0.7
Experiment 3 Experiment 4 Experiment 5
Experiment 8
0.8
0.9
Experiment 9
610
and without SIS to verify the effectiveness of SIS. The parameters of this part are set as follows: Cr =0.8 , Mr =0.1, F =0.8. Table 5 shows the average value and best value of the vehicle number of the solution obtained by HMOEA-GL with SIS and without SIS. Table 5: The vehicle number of HMOEA-GL with SIS and without SIS HMOEA-GL with SIS
HMOEA-GL without SIS
Vehicle Nmuber Mean
Best
Mean
Best
R101
38.6498
31
39.1064
32
R102
34.3476
29
33.6790
29
R103
27.2037
23
27.0000
24
R104
21.3469
19
22.0000
20
C101
21.1860
16
23.7532
18
C102
20.8333
14
21.4593
15
C103
17.9293
13
18.1917
15
C104
13.3077
12
13.7805
12
RC101
33.3000
30
32.8857
30
RC102
30.3056
27
29.5366
27
RC103
26.1190
24
26.0000
24
RC104
22.9444
21
22.9744
21
It can be seen from table 5 that, for the average of vehicle number, in the 615
half of experimental problems HMOEA-GL with SIS is better than HMOEA-GL without SIS. In six problems, HMOEA-GL with SIS is able to find a solution with fewer vehicle number, and on the remaining problems it is the same as with HMOEA-GL without SIS. It can be shown that the addition of SIS can 33
effectively reduce the number of vehicles in the solution. In the following exper620
iments, SIS is added to all the algorithms participating in the comparison. 5.5. Comparisons with other Multiobjective Algorithms Experiment will compare HMOEA-GL, NSGA-II, SPEA2, MOEA/D-a, MOEA/Db and multiobjective evolutionary algorithm with FSS-GS but no RSD-LS (MOEAG). According to the results discussed in section 5.3, the parameters of algo-
625
rithms are set as table 6. Table 6: The parameters of algorithms Parameters
HMOEA-GL
MOEA-G
NSGA-II
SPEA2
MOEA/D-a
MOEA/D-b 50
Population Size
50
50
50
50
50
Elite Population Size
25
25
50
50
/
/
Sub-population Size
25
25
/
/
/
/
Crossover Rate
0.8
0.8
0.8
0.8
0.8
0.9
Mutation Rate
0.1
0.1
0.1
0.1
0.1
0.1
DLS proportion F
0.8
/
/
/
/
/
Number of Neighborhood T
/
/
/
/
5
5
Number of Sub-problems N
/
/
/
/
50
50
The average and standard deviation of convergence indices, distribution indices, and spread indices is shown in table 7 to 11, the result of the hypothesis test between HMOEA-GL and competing algorithms are expressed in the H.T. column. 630
From the GD index as shown in table 7, compared with the other algorithms, among the 12 questions tested, HMOEA-GL performed best on 8 problems and ranked second on 2 problems. MOEA-G performed best on 2 problems and ranked second on 4 problems. This shows that the final solution set obtained by HMOEA-GL and MOEA-G are closer to P F ∗ than the solution set obtained
635
by NSGA-II, SPEA2 and MOEA/D, which means that the convergence performance is better. The result of Wilcoxon rank-sum test can also support these conclusions. From the IGD index and the result of Wilcoxon rank-sum test as shown in table 8, HMOEA-GL performed best on 8 problems again and ranked second 34
Table 7: The experimental results of GD index HMOEA-GL
MOEA-G
NSGA-II
SPEA2
MOEA/D-a
MOEA/D-b
GD Mean
S.D.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
R101 1.46E+01 1.10E+01 1.30E+01 9.64E+00
0
4.70E+01 3.03E+01
1
2.70E+01 2.12E+01
1
1.26E+01 7.14E+00
0
1.09E+01 6.27E+00
0
R102 1.59E+01 1.44E+01 1.55E+01 1.15E+01
0
4.74E+01 2.66E+01
1
2.47E+01 1.55E+01
1
1.01E+01 4.79E+00
0
1.06E+01 6.85E+00
0
R103 1.18E+01 1.85E+01 1.56E+01 1.96E+01
0
4.01E+01 1.98E+01
1
1.66E+01 1.54E+01
1
1.87E+01 1.33E+01
1
1.51E+01 1.37E+01
1
R104 4.21E+00 5.31E+00 7.03E+00 7.28E+00
0
1.64E+01 1.34E+01
1
8.16E+00 5.77E+00
1
1.12E+01 9.40E+00
1
9.79E+00 8.43E+00
1
C101 1.75E+02 1.35E+02 2.21E+02 1.77E+02
0
1.17E+03 3.26E+02
1
4.52E+02 2.41E+02
1
1.96E+03 4.89E+02
1
2.13E+03 5.02E+02
1
C102 8.27E+01 6.79E+01 8.40E+01 5.79E+01
0
5.40E+02 2.29E+02
1
1.73E+02 1.43E+02
1
1.01E+03 3.15E+02
1
9.15E+02 4.99E+02
1
C103 4.51E+01 4.93E+01 8.31E+01 7.23E+01
1
1.68E+02 9.60E+01
1
7.71E+01 6.50E+01
1
3.52E+02 3.41E+02
1
3.74E+02 3.48E+02
1
C104 3.08E+01 4.21E+01 4.02E+01 2.78E+01
1
9.74E+01 7.20E+01
1
3.63E+01 3.94E+01
0
1.12E+02 1.72E+02
0
1.18E+02 1.66E+02
1
RC101 7.66E+00 5.55E+00 6.74E+00 7.10E+00
0
1.49E+01 7.37E+00
1
1.23E+01 5.42E+00
1
1.03E+01 3.14E+00
1
1.02E+01 2.66E+00
1
RC102 6.70E+00 7.13E+00 9.91E+00 1.16E+01
0
1.26E+01 1.52E+01
1
7.84E+00 4.41E+00
1
7.86E+00 3.81E+00
1
8.16E+00 3.08E+00
1
RC103 4.39E+00 5.32E+00 3.82E+00 3.13E+00
0
1.26E+01 1.24E+01
1
8.63E+00 4.17E+00
1
8.17E+00 5.54E+00
1
1.04E+01 7.24E+00
1
RC104 2.60E+00 2.23E+00 5.96E+00 9.09E+00
0
6.74E+00 4.50E+00
1
7.43E+00 3.89E+00
1
6.63E+00 1.53E+00
1
7.17E+00 2.84E+00
1
Table 8: The experimental results of IGD index HMOEA-GL
MOEA-G
NSGA-II
SPEA2
MOEA/D-a
MOEA/D-b
IGD Mean
S.D.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
R101 9.33E+00 4.37E+00 9.24E+00 3.29E+00
0
2.85E+01 1.01E+01
1
1.46E+01 5.10E+00
1
1.49E+01 2.16E+00
1
1.25E+01 1.79E+00
1
R102 5.03E+00 1.95E+00 6.56E+00 3.00E+00
1
1.51E+01 3.16E+00
1
9.63E+00 3.18E+00
1
9.96E+00 2.21E+00
1
9.06E+00 1.39E+00
1
R103 2.75E+00 1.50E+00 3.84E+00 2.09E+00
1
8.93E+00 2.37E+00
1
6.57E+00 1.36E+00
1
7.98E+00 1.57E+00
1
7.62E+00 1.87E+00
1
R104 1.51E+00 1.02E+00 2.22E+00 1.33E+00
1
4.97E+00 1.23E+00
1
5.48E+00 1.17E+00
1
5.89E+00 1.26E+00
1
5.76E+00 1.27E+00
1
C101 2.04E+02 8.98E+01 2.77E+02 1.41E+02
1
9.50E+02 2.40E+02
1
3.80E+02 1.92E+02
1
6.49E+01 2.43E+01
1
6.78E+01 2.89E+01
1
C102 3.99E+01 2.03E+01 4.55E+01 1.86E+01
0
3.23E+02 1.37E+02
1
9.80E+01 1.01E+02
1
5.29E+01 1.83E+01
1
5.54E+01 2.30E+01
1
C103 2.93E+01 1.38E+01 3.34E+01 1.53E+01
0
4.04E+01 1.49E+01
1
2.51E+01 1.18E+01
0
4.71E+01 1.00E+01
1
4.39E+01 1.22E+01
1
C104 1.65E+00 1.05E+00 3.07E+00 1.25E+00
1
4.03E+00 1.20E+00
1
3.27E+00 1.11E+00
1
4.64E+00 1.04E+00
1
4.62E+00 9.88E-01
1
RC101 4.68E+00 1.06E+00 4.24E+00 1.81E+00
0
1.06E+01 1.78E+00
1
8.79E+00 1.42E+00
1
9.69E+00 1.77E+00
1
9.85E+00 1.89E+00
1
RC102 2.60E+00 1.35E+00 2.99E+00 1.43E+00
0
6.49E+00 1.31E+00
1
6.37E+00 1.44E+00
1
7.24E+00 1.66E+00
1
7.05E+00 1.46E+00
1
RC103 2.17E+00 1.04E+00 2.19E+00 9.68E-01
0
5.87E+00 1.21E+00
1
7.38E+00 1.49E+00
1
6.90E+00 1.45E+00
1
7.40E+00 1.34E+00
1
RC104 1.99E+00 9.73E-01 2.35E+00 9.94E-01
0
5.57E+00 1.16E+00
1
6.60E+00 1.19E+00
1
6.45E+00 1.39E+00
1
6.47E+00 1.33E+00
1
35
640
on 4 problems. MOEA-G also performed best on 2 problems and ranked second on 7 problems. IGD index can not only judge the convergence performance of an algorithm, but also consider the performance of distribution performance. The results of IGD experiments show that the final solution set obtained by HMOEA-GL and MOEA-G is not only closer to P F ∗ , but more consistent with
645
the entire P F ∗ . Table 9: The experimental results of HV index
HMOEA-GL
MOEA-G
NSGA-II
SPEA2
MOEA/D-a
MOEA/D-b
HV Mean
S.D.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
R101 2.10E+04 9.38E+02 2.09E+04 6.80E+02
0
1.69E+04 9.81E+02
1
1.85E+04 9.77E+02
1
1.50E+04 9.71E+02
1
1.62E+04 1.04E+03
1
R102 1.52E+04 6.71E+02 1.55E+04 9.02E+02
0
1.22E+04 1.03E+03
1
1.30E+04 1.04E+03
1
1.08E+04 1.29E+03
1
1.15E+04 6.78E+02
1
R103 6.06E+03 4.22E+02 5.99E+03 4.30E+02
0
4.61E+03 3.95E+02
1
4.82E+03 4.16E+02
1
4.02E+03 5.10E+02
1
4.12E+03 5.90E+02
1
R104 1.00E+03 1.12E+02 9.80E+02 9.91E+01
0
6.48E+02 1.19E+02
1
5.66E+02 1.01E+02
1
5.28E+02 1.20E+02
1
5.25E+02 1.30E+02
1
C101 2.80E+05 7.63E+03 2.74E+05 8.59E+03
1
2.18E+05 1.20E+04
1
2.61E+05 1.02E+04
1
1.37E+05 2.23E+04
1
1.45E+05 2.00E+04
1
C102 1.65E+05 5.57E+03 1.64E+05 6.38E+03
0
1.34E+05 7.94E+03
1
1.56E+05 6.26E+03
1
9.25E+04 1.13E+04
1
9.66E+04 1.11E+04
1
C103 6.99E+04 2.57E+03 6.92E+04 1.77E+03
0
5.98E+04 2.81E+03
1
6.62E+04 2.43E+03
1
4.45E+04 3.73E+03
1
4.61E+04 4.36E+03
1
C104 6.91E+03 5.50E+02 7.03E+03 5.00E+02
0
5.43E+03 6.69E+02
1
6.20E+03 5.71E+02
1
3.15E+03 8.74E+02
1
3.29E+03 1.06E+03
1
RC101 1.63E+03 1.11E+02 1.67E+03 1.59E+02
0
1.09E+03 1.34E+02
1
1.22E+03 1.37E+02
1
1.10E+03 1.79E+02
1
1.08E+03 1.88E+02
1
RC102 2.56E+03 2.55E+02 2.52E+03 2.49E+02
0
1.66E+03 2.32E+02
1
1.64E+03 2.84E+02
1
1.41E+03 3.62E+02
1
1.49E+03 3.45E+02
1
RC103 1.45E+03 1.29E+02 1.48E+03 1.39E+02
0
9.86E+02 1.48E+02
1
7.81E+02 1.80E+02
1
8.09E+02 1.81E+02
1
7.81E+02 1.73E+02
1
RC104 6.38E+02 8.91E+01 6.16E+02 7.76E+01
0
3.30E+02 8.79E+01
1
2.28E+02 1.15E+02
1
2.35E+02 1.21E+02
1
2.43E+02 1.07E+02
1
From the HV index and the result of statistical test as shown in table 9, HMOEA-GL still ranked first on 8 problems and ranked second on 4 other problems. MOEA-G achieved the first place on 4 problems and ranked second of other 8 problems. The larger the HV value, the larger the dominance area, 650
which means that the final solution set obtained by HMOEA-GL and MOEA-G have more larger dominance area than the solution set obtained by NSGAII, SPEA2 and MOEA/D. That is, HMOEA-GL and MOEA-G have better performance in convergence and distribution. It can be seen that HMOEA-GL and MOEA-G are not outstanding in the
655
Spacing index and Spread index as shown in table 10 and 11. Spacing index calculates the sum of the variance between the distance of successive solutions and the shortest distance in the entire solution set. It is mainly used to judge the distribution performance of the solution set, whether the distribution of the solution is uniform. Spread not only calculates the distance between successive
36
Table 10: The experimental results of Spacing index HMOEA-GL
MOEA-G
NSGA-II
SPEA2
MOEA/D-a
MOEA/D-b
Spacing Mean
S.D.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
R101 1.62E+01 1.59E+01 1.64E+01 2.08E+01
0
2.03E+01 1.94E+01
0
1.13E+01 7.63E+00
0
1.23E+01 1.38E+01
0
1.26E+01 1.19E+01
0
R102 1.72E+01 2.20E+01 1.29E+01 1.42E+01
0
1.55E+01 1.34E+01
0
1.11E+01 1.06E+01
0
5.58E+00 6.30E+00
1
6.88E+00 1.18E+01
1
R103 1.47E+01 3.05E+01 1.77E+01 3.60E+01
0
2.40E+01 2.45E+01
1
1.03E+01 1.34E+01
0
1.66E+01 2.37E+01
0
1.49E+01 1.70E+01
0
R104 1.68E+00 5.76E+00 3.63E+00 6.31E+00
1
7.83E+00 1.43E+01
1
3.29E+00 9.64E+00
0
6.11E+00 1.15E+01
1
5.21E+00 1.30E+01
0
C101 1.13E+02 8.51E+01 1.18E+02 8.62E+01
0
1.80E+02 1.79E+02
0
1.01E+02 5.28E+01
0
3.39E+02 3.93E+02
1
2.52E+02 1.33E+02
1
C102 1.12E+02 1.05E+02 1.19E+02 9.68E+01
0
1.46E+02 9.59E+01
0
8.30E+01 6.15E+01
0
4.92E+02 5.79E+02
1
4.25E+02 5.45E+02
1
C103 1.01E+02 9.23E+01 1.12E+02 1.18E+02
0
1.45E+02 1.32E+02
0
6.60E+01 5.71E+01
0
3.95E+02 4.58E+02
0
4.39E+02 4.16E+02
1
C104 1.78E+01 2.96E+01 4.38E+01 4.47E+01
1
9.52E+01 1.09E+02
1
2.31E+01 3.67E+01
0
9.41E+01 1.61E+02
0
6.46E+01 1.18E+02
0
RC101 3.89E+00 5.87E+00 3.78E+00 6.08E+00
0
5.63E+00 1.10E+01
0
1.42E+00 4.46E+00
1
1.31E+00 3.19E+00
0
3.11E-01 7.65E-01
1
RC102 5.48E+00 1.03E+01 8.36E+00 1.22E+01
0
9.05E+00 2.91E+01
0
2.44E+00 4.94E+00
0
3.92E+00 6.09E+00
0
1.08E+00 2.64E+00
0
RC103 2.55E-02 8.47E-02 2.18E+00 4.53E+00
0
4.62E+00 1.75E+01
0
1.03E+00 2.06E+00
0
1.32E+00 3.97E+00
0
4.41E-01 1.76E+00
0
RC104 3.45E+00 7.71E+00 1.03E+00 3.67E+00
0
4.06E+00 1.10E+01
0
4.22E-02 8.44E-02
0
0.00E+00 0.00E+00
0
0.00E+00 0.00E+00
0
Table 11: The experimental results of Spread index HMOEA-GL
MOEA-G
NSGA-II
SPEA2
MOEA/D-a
MOEA/D-b
Spread Mean
S.D.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
Mean
S.D.
H.T.
R101
9.88E-01 1.66E-01 9.55E-01 2.02E-01
0
9.21E-01 1.17E-01
0
8.94E-01 1.14E-01
1
9.77E-01 1.41E-01
0
9.34E-01 1.30E-01
R102
9.31E-01 2.01E-01 9.20E-01 2.03E-01
0
9.53E-01 1.14E-01
0
8.59E-01 1.35E-01
0
8.81E-01 1.03E-01
0
8.67E-01 1.03E-01
0 0
R103
8.23E-01 2.34E-01 8.99E-01 1.84E-01
0
9.00E-01 1.65E-01
0
8.59E-01 1.46E-01
0
8.57E-01 1.99E-01
0
9.01E-01 1.05E-01
0
R104
6.13E-01 1.39E-01 7.34E-01 1.59E-01
1
7.44E-01 1.99E-01
1
7.55E-01 1.56E-01
1
8.22E-01 1.33E-01
1
7.57E-01 1.54E-01
1
C101
9.26E-01 6.81E-02 9.44E-01 1.02E-01
0
9.26E-01 8.16E-02
0
8.82E-01 8.02E-02
1
1.10E+00 1.14E-01
1
1.11E+00 1.00E-01
1
C102
9.90E-01 1.61E-01 9.85E-01 1.43E-01
0
9.37E-01 8.47E-02
0
8.90E-01 1.24E-01
1
1.15E+00 1.14E-01
1
1.08E+00 2.17E-01
1
C103 1.06E+00 2.47E-01 1.03E+00 1.53E-01
0
9.91E-01 1.25E-01
1
9.26E-01 1.47E-01
1
1.01E+00 3.03E-01
0
1.01E+00 2.24E-01
0
C104
6.69E-01 1.99E-01 8.93E-01 1.92E-01
1
8.33E-01 2.12E-01
1
7.46E-01 1.72E-01
0
7.08E-01 1.80E-01
0
6.48E-01 1.85E-01
0
RC101 7.89E-01 1.49E-01 7.95E-01 2.20E-01
0
8.18E-01 1.37E-01
0
7.30E-01 1.51E-01
0
8.76E-01 1.10E-01
1
8.26E-01 1.05E-01
0
RC102 7.46E-01 2.36E-01 8.04E-01 1.99E-01
0
7.60E-01 1.62E-01
0
7.90E-01 1.27E-01
0
7.43E-01 1.68E-01
0
7.82E-01 1.24E-01
0
RC103 6.08E-01 1.07E-01 7.14E-01 1.45E-01
1
6.79E-01 1.51E-01
0
7.47E-01 1.33E-01
1
8.74E-01 1.35E-01
1
7.10E-01 1.24E-01
1
RC104 7.09E-01 1.74E-01 6.69E-01 1.41E-01
0
7.68E-01 1.32E-01
0
6.08E-01 4.14E-02
0
8.18E-01 1.15E-01
0
7.21E-01 1.25E-01
0
37
660
solutions, but also considers the smallest distance between edge points and boundary solutions. From Spacing index, HMOEA-GL only performed best on 3 problems and ranked second on 4 problems. MOEA-G ranks third on 6 problems. From the Spread index, HMOEA-GL performed best on 3 problems and ranked second on 3 other problems. MOEA-G ranked second on 2 problems
665
and ranked third on 2 problems. The VRPTW discussed in this paper is not a continuous function problem. We found that when the number of vehicles becomes small, the number of feasible solutions existing in the solution space will decrease. As the solution space approaches the direction in which the number of vehicles decreases, the difference in total waiting time between solutions for
670
different vehicle numbers increases. However, HMOEA-GL and MOEA-G are more prefers to look for the solutions which fewer vehicles number, which may be the reason for their general performance on the Spacing index and Spread index. HMOEA-GL does not use special density maintaining strategies (such as calculate crowding distance) to maintain the distribution of the population,
675
and it is acceptable to obtain such distribution performance. Overall, comparing to NSGA-II, SPEA2 and MOEA/D, HMOEA-GL and MOEA-G have better convergence performance in terms of 12 Solomon benchmark problems. The distribution performance of HMOEA-GL and MOEA-G is acceptable. It can be stated that the proposed FSS-GS and the framework of
680
HMOEA-GL are effective. From the experimental results, it can be found that the significance between HMOEA-GL and MOEA-G is not strong, so in order to more accurately judge the ability of RSD-LS, experiment evaluated the C index for HMOEA-GL and MOEA-G. C index is a binary performance index, can directly reflect the convergence performance of the two algorithms. The
685
results of the C index is shown in table 12. It can be seen that HMOEA-GL performs better on 10 problems than MOEA-G on the C index. This means, in most problems, the final solution set obtained by HMOEA-GL dominates more the solution obtained by MOEAG. It can be concluded that RSD-LS is effective.
690
In summary, HMOEA-GL has been verified on 12 Solomon benchmark prob38
Table 12: Mean and standard deviation of C between HMOEA-GL and MOEA-G C
C(HM OEA−GL,M OEA−G) Mean
C(M OEA−G,HM OEA−GL)
S.D.
Mean
S.D.
R101
0.4323
3.75E-01
0.4908
3.91E-01
R102
0.5898
4.07E-01
0.3526
3.89E-01
R103
0.5931
4.39E-01
0.3220
4.18E-01
R104
0.5006
4.68E-01
0.2722
3.93E-01
C101
0.6195
4.52E-01
0.3578
4.30E-01
C102
0.5479
3.68E-01
0.3889
3.79E-01
C103
0.6860
3.39E-01
0.2671
3.35E-01
C104
0.6719
3.90E-01
0.2194
3.35E-01
RC101
0.3506
3.98E-01
0.5322
4.42E-01
RC102
0.4178
4.66E-01
0.4083
4.62E-01
RC103
0.3556
4.79E-01
0.3500
4.76E-01
RC104
0.4278
4.79E-01
0.2389
4.14E-01
lems with four algorithms, MOEA-G, NSGA-II, SPEA2 and MOEA/D, by five performance indices GD, IGD, HV , Spacing and Spread. And conducted an additional evaluation between HMOEA-GL and MOEA-G on the C index. The experimental results show that HMOEA-GL have more excellent performance 695
in convergence, while maintaining a satisfying distribution performance.
6. Conclusions In this paper, HMOEA-GL is proposed for solving VRPTW. Two objectives are mainly considered, one is minimizing the time-wasting due to early arrival, and another is reducing the number of vehicle. In HMOEA-GL, FSS-GS is de700
signed as the global search strategy. FSS-GS mixes the center area-based elite sampling strategy and edge area-based sampling strategy. Center area-based elite sampling strategy can fast determine which individuals are nondominated or dominated but have large domination area, i.e. an individual located in the central area of the Pareto front. Combining the edge area-based sampling strat-
705
egy which prefers the edge region on the Pareto front, FSS-GS increases the attention to individuals which located on the edge and center regions of the Pareto frontier. FSS-GS can guarantee the capability of HMOEA-GL for con39
verging to the multiply directions of Pareto front quickly and its natural uniform distribution performance. RSD-LS is implemented as a local search strategy af710
ter FSS-GS. RSD-LS can further improve the search ability of HMOEA-GL by guiding the individuals with poor performance to move closer to better one. Permutation code is used to represent the solution of VRPTW, and use the appropriate genetic operator to genetic iteration. In the case of using permutation code, exchange order is used to indicate the difference of route sequence between
715
individuals, which support the execution of RSD-LS smoothly. Designing suitable coding method to represent the individual of VRPTW, and proper genetic operators are used to recombine chromosomes. A single objective optimization method SIS is used to reduce the number of vehicles in one solution. Comparing with NSGA-II, SPEA2, MOEA/D on the 12 Solomon benchmark problems,
720
the experimental results show that the HMOEA-GL hybridizing FSS-GS and RSD-LS is effective. As an effectively intelligent algorithm, HMOEA-GL can be the core of an expert and intelligent decision support system, and provide suitable solutions for decision-makers. However, the proposed method might lose the efficiency in
725
solving the hyperdimensional optimization problems, since when the number of objectives is more than four, the sampling strategy of VEGA in FSS-GS will be ineffective and the computation time of calculating the fitness by PDDR-FF will increase considerably. In the future work, it is appropriate to consider some factors of the problems such as multiple types of goods, multiple types of vehicles,
730
multiple depots and green energy etc. In addition, hybridizing HMOEA-GL with other effective algorithms such as ACO, PSO, TS and SA etc., might be an interesting research direction.
Acknowledgements This research work is supported by the Science & Technology Research 735
Project of Henan Province (162102210044), Program for Science & Technology Innovation Talents in Universities of Henan Province (19HASTIT027), Key Re-
40
search Project in Universities of Henan Province (17A520030), and the Grantin-Aid for Scientific Research (C) of Japan Society of Promotion of Science (JSPS) (19K12148).
740
References Abidi, H., Hassine, K., & Mguis, F. (2018). Genetic algorithm for solving a dynamic vehicle routing problem with time windows. In 2018 International Conference on High Performance Computing & Simulation (HPCS) (pp. 782– 788). IEEE.
745
Baˇ nos, R., Ortega, J., Gil, C., Fernndez, A., & Toro, F. D. (2013). A simulated annealing-based parallel multi-objective approach to vehicle routing problems with time windows. Expert Systems with Applications, 40 , 1696–1707. Banzhaf, W. (1990). The “molecular” traveling salesman. Biological Cybernetics, 64 , 7–14.
750
Braekers, K., Ramaekers, K., & Nieuwenhuyse, I. V. (2016). The vehicle routing problem: State of the art classification and review. Computers & Industrial Engineering, 99 , 300–313. Dantzig, G. B., & Ramser, J. H. (1959). The truck dispatching problem. Management science, 6 , 80–91.
755
Davis, L. (1991). Handbook of genetic algorithms, . Desrochers, M., Desrosiers, J., & Solomon, M. (1992). A new optimization algorithm for the vehicle routing problem with time windows. Operations research, 40 , 342–354. Dixit, A., Mishra, A., & Shukla, A. (2019). Vehicle routing problem with time
760
windows using meta-heuristic algorithms: A survey. In Harmony Search and Nature Inspired Optimization Algorithms (pp. 539–546). Springer.
41
Garcia-Najera, A., & Bullinaria, J. A. (2011). An improved multi-objective evolutionary algorithm for the vehicle routing problem with time windows. Computers & Operations Research, 38 , 287–300. 765
Gen, M., & Cheng, R. (2000). Genetic algorithms and engineering optimization. John Wiley & Sons. Gen, M., Cheng, R., & Lin, L. (2008). Network models and optimization: Multiobjective genetic algorithm approach. Springer Science & Business Media. G¨o¸cken, T., Yaktubay, M., & Kili¸c, F. (2017). Improvement of a genetic algo-
770
rithm approach for the solution of vehicle routing problem with time windows. In Artificial Intelligence and Data Processing Symposium (IDAP), 2017 International (pp. 1–8). IEEE. Gonz´ alez, O. M., Segura, C., Pe˜ na, S. I. V., & Le´ on, C. (2017). A memetic algorithm for the capacitated vehicle routing problem with time windows. In
775
Evolutionary Computation (CEC), 2017 IEEE Congress on (pp. 2582–2589). IEEE. Gupta, A., & Saini, S. (2017). An enhanced ant colony optimization algorithm for vehicle routing problem with time windows. In 2017 Ninth International Conference on Advanced Computing (ICoAC) (pp. 267–274). IEEE.
780
Holland, J. H. (1992). Adaptation in natural and artificial systems. Ann Arbor , 6 , 126137. Jiang, X., Zhou, Q., & Ye, Y. (2017). Method of task assignment for uav based on particle swarm optimization in logistics. In Proceedings of the 2017 International Conference on Intelligent Systems, Metaheuristics & Swarm
785
Intelligence (pp. 113–117). ACM. Lenstra, J. K., & Kan, A. H. G. R. (1981). Complexity of vehicle routing and scheduling problems. Networks, 11 , 221–227.
42
Liu, S. C., & Chen, J. R. (2011). A heuristic method for the inventory routing and pricing problem in a supply chain. Expert Systems with Applications, 38 , 790
13223–13231. Long, J., Sun, Z., Pardalos, P. M., Hong, Y., Zhang, S., & Li, C. (2019). A hybrid multi-objective genetic local search algorithm for the prize-collecting vehicle routing problem. Information Sciences, 478 , 40–61. Marinakis, Y., Marinaki, M., & Migdalas, A. (2019). A multi-adaptive parti-
795
cle swarm optimization for the vehicle routing problem with time windows. Information Sciences, 481 , 311–329. Moon, I. K., Lee, J.-H., & Seong, J. (2012). Vehicle routing problem with time windows considering overtime and outsourcing vehicles. Expert Systems with Applications, 39 , 13202–13213.
800
Pratiwi, A., Pratama, A., Sa’diyah, I., & Suprajitno, H. (2018). Vehicle routing problem with time windows using natural inspired algorithms. In Journal of Physics: Conference Series (p. 012025). IOP Publishing volume 974. Pu, E., Wang, F., Yang, Z., Wang, J., Li, Z., & Huang, X. (2017). Hybrid differential evolution optimization for the vehicle routing problem with time
805
windows and driver-specific times. Wireless Personal Communications, 95 , 2345–2357. Salamatbakhsh-Varjovi, A., Tavakkoli-Moghaddam, R., Alinaghian, M., & Najafi, E. (2018). Robust periodic vehicle routing problem with time windows under uncertainty: An efficient algorithm. KSCE Journal of Civil Engineer-
810
ing, 22 , 4626–4634. Schaffer, J. D. (1985). Multiple objective optimization with vector evaluated genetic algorithms. In Proceedings of the First International Conference on Genetic Algorithms and Their Applications, 1985 . Lawrence Erlbaum Associates. Inc., Publishers.
43
815
Schneider, M. (2016). The vehicle-routing problem with time windows and driver-specific times. European Journal of Operational Research, 250 , 101– 119. Solomon, M. M. (1987). Algorithms for the vehicle routing and scheduling problems with time window constraints. Operations research, 35 , 254–265.
820
Storn, R., & Price, K. (1997). Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces. Journal of global optimization, 11 , 341–359. Tang, L., & Wang, X. (2013). A hybrid multiobjective evolutionary algorithm for multiobjective optimization problems. IEEE Transactions on Evolutionary
825
Computation, 17 , 20–45. Tasan, A. S., & Gen, M. (2012). A genetic algorithm based approach to vehicle routing problem with simultaneous pick-up and deliveries. Computers & Industrial Engineering, 62 , 755–761. Wang, Y., & Xing, L. (2018). A multi ant system based hybrid heuristic algo-
830
rithm for vehicle routing problem with service time customization. In International Conference on Bio-Inspired Computing: Theories and Applications (pp. 411–422). Springer. Wang, Z., Li, Y., & Hu, X. (2015). A heuristic approach and a tabu search for the heterogeneous multi-type fleet vehicle routing problem with time windows
835
and an incompatible loading constraint. Computers & Industrial Engineering, 89 , 162–176. Wolpert, D. H., Macready, W. G. et al. (1997). No free lunch theorems for optimization. IEEE transactions on evolutionary computation, 1 , 67–82. Yong, S., Boudouh, T., & Grunder, O. (2016). A hybrid genetic algorithm for
840
a home health care routing problem with time window and fuzzy demand. Expert Systems with Applications, 72 , 160–176.
44
Yu, B., Yang, Z., & Yao, B. (2011). A hybrid algorithm for vehicle routing problem with time windows. Expert Systems with Applications, 38 , 435–441. Yu, X., & Gen, M. (2010). Introduction to evolutionary algorithms. Springer 845
Science & Business Media. Zhang, H., Zhang, Q., Ma, L., Zhang, Z., & Liu, Y. (2019). A hybrid ant colony optimization algorithm for a multi-objective vehicle routing problem with flexible time windows. Information Sciences, 490 , 166–190. Zhang, J., Yang, F., & Weng, X. (2018). An evolutionary scatter search parti-
850
cle swarm optimization algorithm for the vehicle routing problem with time windows. IEEE Access, 6 , 63468–63485. Zhang, W., Gen, M., & Jo, J. (2014). Hybrid sampling strategy-based multiobjective evolutionary algorithm for process planning and scheduling problem. Journal of Intelligent Manufacturing, 25 , 881–897.
855
Zhou, A., Qu, B. Y., Li, H., Zhao, S. Z., Suganthan, P. N., & Zhang, Q. (2011). Multiobjective evolutionary algorithms: A survey of the state of the art. Swarm & Evolutionary Computation, 1 , 32–49.
45
Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:
Author Contribution Statement Wenqiang Zhang: Conceptualization, Methodology, Writing – Review & Editing, Project Administration, Funding Acquisition, Supervision. Diji Yang: Data Curation, Software, Resources, Visualization, Formal Analysis. Guohui Zhang: Validation, Investigation, Writing – Review & Editing. Mitsuo Gen: Visualization, Validation, Writing – Review & Editing, Supervision.