Three-stage multi-modal multi-objective differential evolution algorithm for vehicle routing problem with time windows

Abstract

In this paper, the mathematical model of Vehicle Routing Problem with Time Windows (VRPTW) is established based on the directed graph, and a 3-stage multi-modal multi-objective differential evolution algorithm (3S-MMDEA) is proposed. In the first stage, in order to expand the range of individuals to be selected, a generalized opposition-based learning (GOBL) strategy is used to generate a reverse population. In the second stage, a search strategy of reachable distribution area is proposed, which divides the population with the selected individual as the center point to improve the convergence of the solution set. In the third stage, an improved individual variation strategy is proposed to legalize the mutant individuals, so that the individual after variation still falls within the range of the population, further improving the diversity of individuals to ensure the diversity of the solution set. Based on the synergy of the above three stages of strategies, the diversity of individuals is ensured, so as to improve the diversity of solution sets, and multiple equivalent optimal paths are obtained to meet the planning needs of different decision-makers. Finally, the performance of the proposed method is evaluated on the standard benchmark datasets of the problem. The experimental results show that the proposed 3S-MMDEA can improve the efficiency of logistics distribution and obtain multiple equivalent optimal paths. The method achieves good performance, superior to the most advanced VRPTW solution methods, and has great potential in practical projects.

Keywords

Logistics distribution vehicle routing problem multi-modal multi-objective optimization three-stage strategy differential evolution algorithm

1. Introduction

With the development of online shopping, logistics distribution, as a key component of online shopping, has received more and more attention. The efficiency of distribution affects customers’ satisfaction directly. Vehicle routing problem (VRP) and its variants [1, 2] have been widely popularized because they can simulate practical applications in various fields. One of the classic variants is VRPTW. The goal of VRP is to design a set of optimal distribution routes for vehicles of a certain size, so as to provide services for customers in logistics distribution. It represents the essence of vehicle allocation and path planning under the lowest cost in logistics distribution. Therefore, it is a key problem in logistics distribution and one of the most widely studied problems in the field of combinatorial optimization. In VRPTW, each customer is associated with a time window, which limits the idle time interval for obtaining services.

Based on the constraints and problem structure of VRPTW, the improvement of one goal may lead to the deterioration of other goals. VRPTW is a combination optimization problem with multiple constraints and multi-objective decisions, and the importance of its objective function varies from field to field. For example, for the food distribution and medical industries, delay time is critical. The freight transport industry can consider the total journey as the key objective to minimize compared with other objectives, because the fuel consumption is proportional to the driving distance. Therefore, from an economic point of view, it is important to minimize the total distance traveled by all vehicles. For small industries, the minimization of the number of vehicles may be the highest priority compared to other goals. When planning the vehicle path, the traditional solution attempts to find a single optimal path without considering the diversity of the path. However, the decision-maker hopes to obtain multiple solutions (i.e. multiple equivalent optimal paths) that meet the target requirements at the same time. On the one hand, decision-makers can better understand the problems to be optimized. On the other hand, if a specific solution becomes infeasible due to changes in the environment, the decision-maker can easily change to another equivalent solution to ensure the stability of the decision, that is, at least two equivalent global Pareto optimal solutions correspond to the same point on the Pareto Front (PF) [3, 4]. Therefore, for some problems, it is important to obtain all equivalent Pareto Sets (PS). When solving VRPTW, a multi-objective optimization problem, there may be two or more global or local PSs, some of which may correspond to the same PF. If the decision-maker hopes to get multiple solutions under the constraint conditions, VRPTW can be regarded as a multi-modal multi-objective optimization problem (MMOP), that is, multi-modal multi-objective VRPTW (MMVRPTW).

Figure 1 shows a simple multi-mode multi-objective vehicle path planning scheme [5, 6], which includes cost and time objectives. In Fig. 1, the triangle and circle in the left path diagram represent different forms of transfer stations, and the dotted line in the right figure represents the PF of the problem. Route 1 and Route 4, Route 2 and Route 3 have the same time and cost respectively. Obviously, for passengers who need to complete the trip in a shorter time, there are two options: Route 2 and Route 3, and for passengers who need to complete the trip at less cost, there are also two options: Route 1 and Route 4. When the transit station of a route scheme is affected by weather and other factors, passengers can choose alternative schemes. If the final result can only provide one PS for each passenger, it can not meet the needs of different passengers. There are at least four equivalent Pareto optimal solution sets for the problem: {Route 1, Route 2} {Route 1, Route 3} {Route 4, Route 2} {Route 4, Route 3}. In this problem, a set of Pareto optimal solution sets can not well meet the needs of decision makers. Therefore, it is of great practical significance to study how to obtain more PSs at the same time, so as to provide more choices for decision makers.

Figure 1.

MMOP in path planning.

Since Solomon et al. [7] first proposed the VRPTW and gave a classical example of this kind of problem, VRPTW has become a research hotspot, which is considered as a typical variant of VRP. In recent years, researchers have proposed a large number of multi-objective optimization algorithms to solve VRPTW and its related variants. However, due to the complexity of problem modeling, the difficulty of solving, and the multi-modality of the problem, the research results are relatively few. Xu et al. [8] proposed hybrid genetic algorithm and particle swarm optimization algorithm for vehicle routing problem with time window, and used particle real number coding method to decode the path, thus reducing the computational burden. At the same time, it is combined with the crossover operator of genetic algorithm to avoid falling into local optimum. Jakub Nalepa et al. [9] solved the VRPTW by using the adaptive parameter modulus algorithm. Jose et al. [10] solved the heterogeneous vehicle routing problem with time windows (HVRPTW) using the ACS based memetic algorithm. Wang et al. [11] proposed a multi-objective multi station vehicle routing problem with time windows, and developed a two-stage multi-objective evolutionary algorithm (TS-MOEA) to solve the problem. Lu et al. [12] considered the VRP with multiple stops and multiple trips with time windows and delivery times, established a mixed integer programming model, and proposed a hybrid particle swarm optimization algorithm and a hybrid genetic algorithm to solve the model. Cueto et al. [13]combined the branch cutting algorithm with the heuristic algorithm to solve VRPTW with multiple parking lots and multiple trips. Cheng et al. [14] proposed a solution framework that combines a multi task optimization framework with multi-objective evolutionary algorithms for solving multi-objective VRPTW. Shu et al. [15] presented a two-stage multi-objective evolutionary algorithm based on classified population (TSCEA) to solve a three-objective VRPTW. In the first stage, a population is explored using the proposed algorithm and then classified according to the number of vehicles, we call this process population classification; In the second stage, Pareto solution set of tri-objective VRPTW is obtained by optimizing the classified population again. Hou et al. [16] proposed a multi-objective differential evolution algorithm for solution evaluation (differential mutation strategy based on feasible solutions) for multi-objective VRPTW. Cai et al. [17] proposed a hybrid evolutionary multitask algorithm, termed HEMT, to address MOVRPTWs under the framework of evolutionary multitasking, where multiple MOVRPTWs are optimized simultaneously by leveraging the similarity between them. Srivastava et al. [18] proposed two evolutionary approaches, viz., a steady-state grouping genetic algorithm and a discrete differential evolution algorithm, to solve a three-objective VRPTW. According to the analysis of previous studies, although some research achievements have been made on VRPTW, there are still some problems: (1) The solution objective is single. At present, the research on VRPTW mainly solves the objective from a certain angle, such as the minimum total cost, the minimum travel time, the minimum average waiting time of customers, and the minimum travel time. (2) The final solution is single. At present, the research on VRPTW optimizes multiple objectives in the form of a single objective, so only one final solution can be returned to the decision-makers.

The research on the methods of solving the optimal vehicle routing is endless. Most methods need to establish mathematical models and complete the optimization of vehicle routing by defining different types of variables, constraint functions and objective functions. Common methods mainly include precise search algorithm and heuristic search algorithm [19, 20, 21, 22, 23]. Among them, the precise search algorithm is to establish a corresponding mathematical model for a specific problem, and then use mathematical methods to solve it. It must be able to find the optimal solution of the problem, mainly focusing on branch and bound method [24], branch and cut method [25], etc. Since VRPTW is a NP hard problem, strict mathematical methods must be introduced when using accurate algorithms to solve it. Therefore, too much computing power and storage space will be consumed during calculation, which will limit the accuracy of the optimal vehicle path. This method can only be applied to small-scale VRPTW solutions. With the gradual increase of the problem scale, some scholars propose to use heuristic search algorithm to solve the vehicle routing problem. The heuristic algorithm is proposed based on the optimization algorithm, and its basic idea is to give a feasible solution of the combination optimization problem to be solved within an acceptable range. The heuristic search algorithms mainly used in VRPTW include ant colony algorithm [10], genetic algorithm [12], etc. Compared with the precise search algorithm, the heuristic search algorithm has better robustness and feasibility when dealing with large-scale VPRTW. Based on the existing research, it can be seen that precise search algorithm and intelligent heuristic search algorithm can be used to solve VRPTW and related problems. The exact algorithm can find the optimal solution of the problem, but it is highly dependent on the solution space, the number of constraints and the number of decision variables in the problem model, and cannot provide a general solution strategy for different types of variables, goals and constraints. By designing a heuristic function, the heuristic algorithm can obtain the optimal solution of the search problem in a very short time. The heuristic algorithm can further improve the accuracy of vehicle routing.

In the vehicle path planning problem, the decision-makers hope to get multiple paths with the same target value. In addition, some decision-makers can accept other paths that are slightly worse than the best solution when possible. Therefore, it is of great significance to improve the diversity of solutions obtained. In order to solve this problem, this paper first explains the importance of solving VRPTW using multi-modal multi-objective algorithm. Then, a 3-stage multi-modal multi-objective differential evolution algorithm (3S-MMDEA) is proposed to maintain the diversity of solutions in the decision space. Finally, the effectiveness of the algorithm is verified by experiments. The main contributions of this paper are as follows:

In order to improve the efficiency of logistics distribution and customer satisfaction, and meet the planning needs of different decision-makers for different vehicle routes, a multi-modal multi-objective vehicle routing problem model with time window constraints, namely MMVRPTW, is established.

Aiming at MMVRPTW, a 3S-MMDEA is proposed to solve the problem. In the first stage strategy, the GOBL strategy is used to generate a reverse population, expand the search scope of subsequent individuals, and improve the diversity of individuals.

In the second stage, a search strategy of reachable distribution area is proposed to divide the population to reduce the scale of logistics distribution and improve the efficiency of logistics distribution.

In the third stage strategy, an improved individual variation strategy is proposed to make individuals fall within the range of the population after mutation, further improving the diversity of individuals, so as to ensure the diversity of the solution set. Based on the above strategies, the vehicle paths are optimized through forward output and feedback, and the optimization of MMVRPTW is completed cooperatively.

2. Relevant theories and model establishment

2.1 Relevant theories

The multi-objective optimization problem (MOP) is defined as follows [26]:

$\displaystyle\min\,f(x)=[f_{1}(x),f_{2}(x),\ldots,f_{m}(x)].$ $\displaystyle\textit{s.t. }g_{i}(x)\leqslant 0,i=1,2,\ldots,k;$ (1) $\displaystyle h_{j}(x)=0,j=1,2,\ldots,p.$

Where, $m$ is the number of targets to be optimized, $x=(x_{1},x_{2},\ldots,x_{n})$ is an $n$ -dimensional decision variable, $g_{i}(x)\leqslant 0(i=1,2,\ldots,k)$ is $k$ inequality constraints, $h_{j}(x)=0(j=1,2,\ldots,p)$ is $p$ equality constraints. The space $R^{n}$ satisfying the above equality and inequality constraints is called the decision space, and the space $R^{m}$ mapped by the objective function of the decision space is called the objective space. Suppose $x_{1}$ and $x_{2}$ are two feasible solutions to the optimization problem, if and only if $\forall i=1,2,\ldots m$ , $f_{i}(x_{1})\leqslant f_{i}(x_{2})\wedge\exists j\in[1,m]$ , $f_{j}(x_{1})<f_{j}(x_{2})$ , $x_{1}$ dominates $x_{2}$ [27]. If a solution is not dominated by any other solution, it is called a non dominated solution. All non dominated solutions form a Pareto optimal set (PS). The set of PS mapped to the objective space is called Pareto front (PF).

Definition (multi-modal multi-objective optimization, MMOP) [28, 29]: If a MOP satisfies one of the following two conditions, the problem is called multi-modal multi-objective optimization problem:

The problem has at least one local Pareto optimal solution set.

The problem has at least two equivalent global Pareto optimal solutions, which correspond to the same point on PF.

The objective function values corresponding to the Pareto optimal solutions For the problem of two objectives, the Pareto optimal frontier is usually a line. For multiple targets, the Pareto optimal frontier is usually a hyper-surface. The local Pareto optimal solution refers to the solution that is not dominated by any neighborhood solution, and the global Pareto optimal solution refers to the solution that is not dominated by any solution in the feasible region.

MMOP is a special kind of MOP, which is mainly manifested in two cases: (1) each solution in the decision space has multiple equivalent solutions. (2) There are several equivalent solutions in the middle decomposition of decision space.

The goal of MMOP is not only to obtain the PF approximation with good approximation and diversity in the objective space, but also to obtain enough equivalent solutions in the decision space. In order to improve the search ability of the algorithm and ensure the diversity of decision space and objective space at the same time, researchers have proposed different solutions, mainly including methods to improve the diversity of decision space and methods to improve the search ability [30, 31, 32, 33, 34, 35].

In the process of solving MMOP, designers need to pay more attention to the performance in the decision space, that is, they hope to find the equivalent PS of multiple groups of uniform distribution corresponding to a group of PFs. Multi-modal multi-objective optimization (MMO) can provide more elegant solutions, and provide diversified decisions for decision-makers in real optimization problems. How to balance the convergence, diversity and feasibility in the search process is the key and main difficulty in solving MMOP.

2.2 Basic theory of differential evolution algorithm

Differential evolution (DE) algorithm [36] was proposed by Storn and Price in 1997. As a new intelligent algorithm, it has simple principle, fewer controlled parameters, good robustness and easy implementation. Its essence is a multi-objective (continuous variable) optimization algorithm, which is mainly used to solve the overall optimal solution in multi-dimensional space.

The basic idea of DE algorithm is [37, 38]: The difference component of two individual vectors randomly selected from the population is used as the disturbance of the third random benchmark vector to obtain the variation vector, and then the variation vector is hybridized with the benchmark vector (or target vector) to generate the trial vector. Finally, the benchmark vector competes with the test vector, and the better one is saved in the next generation group. In this way, the differential evolution algorithm improves the population quality generation by generation and guides the population to focus on the location of the optimal solution.

Similar to other population-based stochastic algorithms, the process of DE algorithm mainly includes four steps: population initialization, variation, crossover and selection. As shown in Fig. 2, initialization is a one-time process, and the other three mechanisms are repeated in the DE searching process in the D-dimensional solution space until the termination conditions are met.

Figure 2.

Consecutive phase of DE algorithm.

(1) Initialization

Let D be the dimension of the individual, NP be the population size, and $t$ be the evolutionary algebra, $X(t)$ is the population of generation $t$ . Firstly, the population of generation 0 is randomly generated in the decision space of the problem:

$\displaystyle X(0)=\{x_{1}(0),x_{2}(0),\ldots,x_{NP}(0)\}$ (2)

Where, ${\rm{\bf x}}_{i}(t)=(x_{i,1}(t),x_{i,2}(t),\ldots,x_{i,D}(t))$ is used to represent the $i^{\text{th}}$ individual in the $t$ generation population.

The value of each dimension of an individual can be generated according to the following formula:

$\displaystyle x_{i,j}(0)=L_{j}+\textit{rand}_{i,j}[0,1](U_{j}-L_{j})$ (3)

Where, $1\leqslant i\leqslant NP$ , $1\leqslant j\leqslant D$ , $[L_{j},U_{j}]$ is the value range in the $j^{\text{th}}$ dimension, $\textit{rand}_{i,j}[0,1]$ is a uniformly distributed random number between 0 and 1.

(2) Variation

From a biological point of view, variation means a change in the sequence of genes in the chromosome. In the field of evolutionary computing, variation is regarded as the change of an element. In the differential evolution algorithm, take the simplest variation operation (DE/RAND/1) as an example, the variation operation for the $i^{\text{th}}$ benchmark vector ${\rm{\bf x}}_{i}^{t}$ is: Randomly select three vectors $x_{r1}(t),x_{r2}(t),x_{r3}(t)$ from the current population. $r_{1},r_{2},r_{3}$ is required to be mutually different integers randomly selected from set $\{1,2,\ldots,NP\}\backslash\{i\}$ , and the difference between the two vectors is scaled and added with the third vector to obtain a variation vector ${\rm{\bf V}}_{i}(t)$ :

$\displaystyle{\rm{\bf V}}_{i}(t)=x_{r1}(t)+F(x_{r2}(t)-x_{r3}(t))$ (4)

Where $F$ is the scaling factor of DE, and the value range is $[0,1]$ .

Obviously, the smaller the difference vector between $x_{r2}(t)$ and $x_{r3}(t)$ , the smaller the disturbance. This means that in the initial stage of the algorithm, due to the large difference between individuals and large disturbance, the algorithm searches in a relatively large range. In the later stage of the algorithm, because the groups are close to the optimal individuals and the disturbance value is small, the algorithm searches in a small range.

At present, the other four widely used variation methods [38] are as follows:

DE/rand/2: ${\rm{\bf V}}_{i}(t)=x_{r1}(t)+F[(x_{r2}(t)-x_{r3}(t)+(x_{r4}(t)-x_{r5}(t))]$ DE/best/1: ${\rm{\bf V}}_{i}(t)=x_{\textit{best}}(t)+F((x_{r2}(t)-x_{r3}(t))$ DE/best/2: ${\rm{\bf V}}_{i}(t)=x_{\textit{best}}(t)+F[(x_{r2}(t)-x_{r3}(t)+(x_{r4}(t)-x_{% r5}(t))]$ DE/current-to-best/1: ${\rm{\bf V}}_{i}(t)=x_{i}(t)+F[(x_{\textit{best}}(t)-x_{i}(t)+(x_{r1}(t)-x_{r2% }(t))]$

Where, $r1,r2,r3,r4,r5$ are mutually different integers randomly selected from set $\{1,2,\ldots,NP\}\backslash\linebreak\{i\}$ , $x_{\textit{best}}(t)$ is the best individual in the current group.

(3) Crossover

In order to improve the diversity of evolutionary population, the discrete crossover operator is introduced. Unlike other evolutionary algorithms, which exchange genes based on multiple reference vectors from the parent, the crossover operator in differential evolution algorithm uses reference vector and variation vector to operate.

(4) Selection

The selection operation of DE is a greedy selection mechanism [38], which is to retain the best in both the target vector and its corresponding test vector, so that the fitness value of the offspring is always better than that of the parent, resulting in the population always evolving to the position of the optimal solution and gradually focusing on the position of the optimal solution or satisfactory solution.

2.3 Model of MMVRPTW

VRP is a typical MMOP. When solving this problem, the goal is to find one solution that may be sufficient to obtain an acceptable solution. If multiple shortest paths are not determined, the decision-maker may not be able to consider solutions that may improve performance. Since the preference of decision-makers is a priori unknown, it is necessary to provide a variety of excellent solutions. Generally, in MMVRPTW, each vehicle starts from a central depot, serves customers in the distribution area, and returns to the same central depot after the service. Among them, each car has its own capacity, and each customer has its own demand and time window to obtain services. A soft time window is used here to better optimize costs. The MMVRPTW is related to many factors. The mathematical model established is very complex and has many constraints. In order to facilitate modeling, the following assumptions are made in this study:

Only consider the logistics distribution of a single logistics distribution center.

Vehicles responsible for logistics distribution must start from the distribution center and return to the distribution center after completing all customer order distribution tasks.

Each vehicle only completes the distribution of one line.

The demand and location coordinates of each customer are known and fixed.

The customer nodes are interconnected, that is, the distribution vehicles can be transferred from customer node $i$ to customer node $j$ , or from customer node $j$ to customer node $i$ .

The arc formed between customer nodes has two-way weight, representing distance and time cost.

In the process of vehicle distribution, the impact of temporary vehicle failure or goods distribution error and other factors will not be considered for the time being.

Based on the above assumptions, MMVRPTW can be formulated as an optimization problem in a complete directed graph. The goal of MMVRPTW is to find a solution with multiple equivalent optimal path solutions to serve all customers, meet all constraints and minimize the following five goals. In order to serve all customer points, the following constraints need to be followed: maximum delay time constraint, that is, each customer should be served before the maximum allowable delay time, and the service request should be fully satisfied at one time. Vehicle loading constraint, that is, the total demand of customers in the route shall not exceed the vehicle capacity.

$f_{1}$ number of vehicles:

$\displaystyle f_{1}=R=\left\lceil{\sum\limits_{i=1}^{n}{\frac{q_{i}}{w}}}\right\rceil$ (5)

Wherein, $R$ represents the number of vehicles required to complete the distribution task, $n$ represents the number of customer points, $q_{i}$ represents the distribution volume (or demand volume) of customers $i$ , $w$ represents the maximum carrying capacity of vehicles, and $\lceil\,\rceil$ represents rounding up.

$f_{2}$ total travel distance:

$\displaystyle f_{2}=\min\sum\limits_{k\in R}{R_{k}}$ (6)

Wherein, $R_{k}$ refers to the $k^{\text{th}}$ route, i.e. the route in the charge of the $k^{\text{th}}$ vehicle.

$f_{3}$ makespan:

$\displaystyle f_{3}=\max_{k\in R}T(R_{k})$ (7)

Wherein, $T(R_{k})$ represents the travel time of path $R_{k}$ . $T(R_{k})$ consists of travel time, waiting time and service time of all sides of path $R_{k}$ .

$f_{4}$ total waiting time:

$\displaystyle f_{4}=\sum\limits_{k\in R}{WT(R_{k})}$ (8)

Wherein, $WT(R_{k})$ represents the waiting time of the route $R_{k}$ . If the vehicle arrives at the customer’s point before the time window, the vehicle shall wait for work. $WT(R_{k})$ consists of the waiting time of all customer points in the path $R_{k}$ .

$f_{5}$ total delay time:

$\displaystyle f_{5}=\sum\limits_{k\in R}{DT(R_{k})}$ (9)

Wherein, $DT(R_{k})$ represents the delay time of the route. If the vehicle arrives at the customer point after its time window, the delay time will be generated. $DT(R_{k})$ consists of the delay time of all customer points in the path $R_{k}$ .

Figure 3.

Three stage framework and population distribution.

Figure 4.

Flow chart of 3S-MMDEA.

3. 3S-MMDEA for MMVRPTW

3.1 Basic framework of 3S-MMDEA

Traditional MOP solutions will be more likely to choose solutions with better convergence to form a new group. However, for problems with multiple different solutions corresponding to the same target value, convergence first will lead to premature emergence of the optimal solution, and the lack of diversity maintenance mechanism will lead to poor quality in finding all equivalent solutions.

MMDEA pursues two basic but often conflicting goals, namely convergence and diversity. This conflict has a negative impact on the optimization process of the algorithm and may be aggravated in many objective optimization. Therefore, in the process of MMDEA optimization, it is necessary to deal with the balance between convergence and diversity. In order to handle the balance well, the 3S-MMDEA shown in Fig. 3 is used to solve MMVRPTW. In 3S-MMDEA, a three-stage strategy is proposed to improve the convergence and diversity of MMDE algorithm.

3S-MMDEA consists of the following parts: population initialization, population expansion, individual selection, population division, individual variation, crossover and environmental selection. In the whole differential evolution process, in the first stage, a reverse population is generated through GOBL strategy (that is, multiple customer points are generated reversely), which increases the range of individual screening and improves the diversity of solution set. In the second stage, the reachable distribution area search strategy is used to divide the population (all customer points) with the selected individual (that is, the starting customer point) as the center point to improve the convergence of the solution set. In the third stage, through the improvement of individual variation strategy, the individual is legalized to make the individual still fall within the range of the population, further improving the diversity of the solution set. Based on the synergy of the above three stages of strategies, the diversity of individuals is ensured, so as to improve the diversity of solution sets, and multiple equivalent optimal paths are obtained to meet the planning needs of different decision-makers. The basic flow of 3S-MMDEA is shown in Fig. 4. It should be noted that the existing studies [11, 39, 40] and 3S-MMDEA use different staged methods to solve different types of problems based on different motivations.

3.2 Stage I: GOBL for population initialization

The initialization of the population is a key step in solving the problem, because it can control the quality of the final solution set, and the initialization process is an irreversible process. It is very important to conduct the operation of the non repeatable process, because the subsequent individual selection is based on the initial population, and the individual selection affects the overall quality of the initial solution obtained, thus affecting the overall process of the solution algorithm.

Generalized opposition-based learning (GOBL) [41] can effectively use the resources of the reverse solution, expand the search direction and improve the search efficiency.

Suppose $P=(z_{1},z_{2},\ldots,z_{D})$ is a point in a $D$ -dimensional space (assuming $P$ is a candidate solution), where $z_{1},z_{2},\ldots,z_{D}\in R$ and $z_{j}\in[L_{j},U_{j}]$ , $L$ and $U$ represent the upper and lower bounds, $\forall j\in\{1,2,\ldots,D\}$ . $f(\bullet)$ is the fitness value of the objective function of the candidate solution, then the reverse point of $P$ is $\bar{P}=(\bar{z}_{1},\bar{z}_{2},\ldots,\bar{z}_{D})$ , where $\bar{z}_{j}=k(L_{j}+U_{j})-z_{j}$ , $k=\textit{random}(0,1)$ . If $\bar{z}_{j}<L_{j}$ or $\bar{z}_{j}>U_{j}$ , then $\bar{z}_{j}=\textit{random}(L_{j},U_{j})$ . If $f(\bar{P})<f(P)$ , it means that $\bar{P}$ has a better fitness value than $P$ . At this time, select $\bar{P}$ instead of $P$ , otherwise it remains unchanged.

In the process of evolution, GOBL can restrict the search space of the population, accelerate the convergence speed and efficiency of the algorithm, and quickly find the optimal solution. The pseudo code in the initialization phase is as follows:

Algorithm 1: GOBL strategy in population initialization stage
1 Randomly generate initialization population $x(0)$ # The population number is NP
2 for $i=$ 1:NP
3 for $j=$ 1:D
4 $X(1)_{i,j}=k\times(L_{j}+U_{j})-X(0)_{i,j}$ ; $0<k<1$
5 if $X(1)_{i,j}<L_{j}$ or $X(1)_{i,j}>U_{j}$ then # If it exceeds the boundary, it will be regenerated
6 $X(1)_{i,j}=\textit{rand}(L_{j},U_{j})$
end if
end for
end for

Using GOBL strategy to obtain the original population and the reverse new population, at this time, the population number becomes two, greatly expanding the range of individual selection. The individual selection needs to be conducted in two populations at the same time. The traditional DE algorithm compares the fitness of each individual of the parent vector in the population with the corresponding individual in the new population one by one, and saves the individuals with better fitness. Under this screening mechanism, there may be two individuals with better fitness or worse fitness for comparison. At this time, an individual with poor fitness will be retained, which will undoubtedly affect the evolution speed and easily fall into local optimization. In this paper, by using the generation method of difference vector and the calculation method of crowding distance in [42], the diversity of individuals is considered in both the current population and the changing population to achieve the optimal selection of individuals.

3.3 Stage II: Search strategy of reachable distribution area for population division

When using MMDE algorithm to optimize vehicle routing, customers need to be searched. That is, after individual selection, it is necessary to consider the access order of subsequent individuals to complete traversal search. Common search strategies include blind search strategy and heuristic search strategy [43]. The blind search strategy follows the specified route and does not use heuristic information related to the problem. The heuristic search strategy needs to use heuristic information related to the problem and guide the search process with these heuristic information. Combined with the research in this paper, heuristic search strategy should be adopted when searching individuals (i.e. customer points). In addition, we still need to consider the search range of individuals. On the one hand, because the experimental results are path length oriented, we need to consider the distance between different customer nodes. On the other hand, we need to consider the distribution range, that is, the maximum transportation distance of vehicle distribution. Therefore, this paper proposes a search strategy based on the reachable distribution region. The effectiveness of this method has been verified in our published paper [40].

Based on directed graph $G(V,E)$ , $V$ represents customer node set (i.e. individual set), $E$ represents path set, Dijkstra algorithm [40] can be used to calculate the shortest path of any two nodes (individual). Assuming $S$ is a partitioned set, initially empty, the distribution center $v_{0}$ is placed in the set, and the distance $\textit{dis}[j]$ from $v_{0}$ to any customer node $j$ is calculated. The shortest corresponding customer node $\textit{dis}[j]$ and the current vehicle load are recorded, and the node is included in the set $S$ . Then, the distance from the current node to the other nodes is calculated based on the current node. The above process is repeated until the load in the set $S$ is greater than the vehicle load, and the current set $S$ is output as a partition. Based on the vehicle’s carrying parameters (mainly including carrying capacity and maximum traveling distance), the node set of the reachable distribution area of the next object node corresponding to each individual can be determined, so as to achieve population division. When optimizing the search, this search strategy can ensure that VRPTW decisions conform to vehicle constraints, and effectively improve the search efficiency and speed. The pseudo code of the strategy in this phase is shown in Algorithm 2.

Algorithm 2: Population division based on search strategy of reachable distribution area
1 Initial node $S\leftarrow{\{}v_{0}{\}}$
2 $\textit{dist}\leftarrow\textit{Edge}[0],j=1,2,\ldots,n-1$ ;
3 //Find the node $k$ corresponding to the minimum path:
4 $\textit{dist}[k]==\min{\{}\textit{dist}[i]{\}},i\in G-S$ ;
5 $S\leftarrow S\cup{\{}k{\}}$ ; //Mark point $k$ to enter set $S$
6 //For each $i\in G-S$ modification:
7 $\textit{dist}[i]\leftarrow\min{\{}\textit{dist}[i],\textit{dist}[k]+\textit{% Edge}[k][i]{\}}$
8 Judge: If $\sum{q_{i}}$ in $S$ is greater than the vehicle load, output the node set to realize a division
9 Otherwise turn to Line 3

3.4 Stage III: Improved individual variation strategy

In the variation step, it is better to have strong exploration ability because multiple solution sets need to be obtained. Therefore, DE/rand/2 [38] is adopted in this paper. The difference vector is generated into:

$\displaystyle v_{i}=x_{r_{1}}+F\cdot[(x_{r_{2}}-x_{r_{3}})+(x_{r_{4}}-x_{r_{5}% })]$ (10)

Wherein, $v_{i}$ represents the difference vector and $r_{1},r_{2},r_{3},r_{4},r_{5}$ are an integer unequal to each other. $F$ is the scale factor used to scale the difference vector.

Due to the particularity of the problem, the mutated vector individual may not be feasible, that is, the individual generated by the variation does not meet the boundary conditions and directly exceeds the scope of the solution space, that is, falls outside the search space. Repairing the infeasible solution [41, 42] is a common means of evolutionary algorithm. In this paper, a new repair operator is designed to deal with the individuals crossing the boundary. The processing methods are as follows:

Firstly, the average Euclidean distance $d_{\textit{avg}}$ between the first cross-border individual and other variation completed individuals is calculated.

Secondly, the individual $x^{\prime}$ whose Euclidean distance from the first out-of-boundary individual is the closest distance to $d_{\textit{avg}}$ .

Then, for each dimension of the individual $x^{\prime}$ which is out of the allowed boundaries, the value of that dimension in that individual is mutated for the second time as follows:

$\displaystyle{\rm{\bf V}}_{i,j}=x_{r1,j}-F[(x_{r2,j}-x_{r3,j})+(x_{r4,j}-x_{r5% ,j})]$ (11)

If the just mutated dimension is still out of the allowed boundaries, repair according to the following formula:

$\displaystyle v_{i,j}=\left\{{\begin{array}[]{l}\min\{U_{j},2L_{j}-v_{i,j}\},v% _{i,j}<L_{j};\\ \max\{L_{j},2U_{j}-v_{i,j}\},v_{i,j}>U_{j}.\\ \end{array}}\right.$ (12)

Where, $v_{i,j}$ represents the value of the $i^{\text{th}}$ individual in the $j^{\text{th}}$ dimension, $U_{j}$ and $L_{j}$ represent the upper and lower bounds of the decision space.

The reason for using the above repair operator is that when the better NP individuals are selected from the original population and reverse new species population obtained by GOBL strategy for generation, there may be randomly generated individuals or individuals at the population boundary in the population, which may exceed the range of solution space after individual variation, resulting in fewer individuals, reducing individual diversity, and finally reducing the diversity of optimal solutions in PS. Figure 5 illustrates the legalization of mutant individuals. Where, black solid dots represent five individuals generated by difference, and hollow dots represent individuals generated by the first variation. Figure 5a illustrates that individual $A_{i}$ exceeds the solution space after variation, and calculates the average Euclidean distance between $A_{1}$ and other individuals after mutation $d_{\textit{avg}}=(d_{A_{1}B_{1}}+d_{A_{1}C_{1}}+d_{A_{1}D_{1}}+d_{A_{1}E_{1}})% /4=19/2$ , which less than the solution space boundary, then take the individual with the Euclidean distance closest to $19/2$ from $A_{1}$ in the solution space as $A^{\prime}$ , $A^{\prime}$ falls between $B_{1}$ and $C_{1}$ , and conduct the secondary variation of $A^{\prime}$ according to Eq. (11) to obtain the individual with secondary variation $A_{1}^{\prime}$ , $A_{1}^{\prime}$ falls in the solution space, then ${\{}B_{1},C_{1},D_{1},E_{1},A_{1}^{\prime}{\}}$ constitutes the set of individuals after variation in the understanding space to ensure the diversity of individuals in the understanding space. Figure 5b shows that if individual $A$ is out of solution space after variation, and the average Euclidean distance between the mutated individual $A_{1}$ and other mutated individuals $d_{\textit{avg}}=(d_{A_{1}B_{1}}+d_{A_{1}C_{1}}+d_{A_{1}D_{1}}+d_{A_{1}E_{1}})% /4=21/2$ . Similarly, take the individual whose Euclidean distance from $A_{1}$ is closest to $21/2$ as $A^{\prime}$ , $A^{\prime}$ still falls outside the solution space, and make a secondary variation of $A^{\prime}$ according to Eq. (11) to obtain $A_{1}^{\prime}$ , but $A_{1}^{\prime}$ still falls outside the solution space. At this time, repair $A_{1}^{\prime}$ according to Eq. (12) to obtain $A_{1}^{\prime\prime}$ in the solution space, then ${\{}B_{1},C_{1},D_{1},E_{1},A_{1}^{\prime\prime}{\}}$ constitutes the individual after variation in the understanding space to ensure the diversity of individuals in the understanding space.

Figure 5.

Schematic diagram of variation legalization.

3.5 Crossover and environmental selection

In order to improve the diversity of population, discrete crossover operator is introduced into DE algorithm. Different from the crossover operators in other evolutionary algorithms that exchange genes based on multiple benchmark vectors from the parent, the crossover operators in differential evolutionary algorithms use benchmark vectors and variation vectors to operate.

The selection operation of DE algorithm is a selection mechanism based on greed, which is to retain the optimal value between the target vector and its corresponding test vector, so that the fitness value of the child individual is always better than the fitness value of the parent individual, thus causing the population to always move towards the location of the optimal solution, and gradually focus on the location of the optimal solution or the location of the satisfactory solution. The commonly used selection mechanism may have the situation that the desirable number of individuals in the first level and the second level in the non dominated ranking is opposite. Therefore, this paper uses the environmental selection mechanism in [42] to screen individuals, so that the obtained PF becomes more complete with the increase of iterations.

4. Experimental results and analysis

4.1 Experimental settings

Parameters: In GOBL strategy, $k$ is set as 0.2, the initial value of shrinkage factor $F$ is 0.5, and the crossover probability $C R$ is 0.5. In order to verify the effectiveness of the proposed method, we use MATLAB language for experimental simulation. The computer is configured as Intel Core i7-3630QM 2.40GHz, 8GB RAM, and executed in Windows 7 system.

Data-set: In this study, the data set in the same actual situation as that used by Zhou and Wang [46] and Castro Gutierrez et al. [47] is selected. This data set is composed of 45 real VRP data instances. These data are composed of 3 customers of different sizes, 5 time window configurations and 3 vehicle models of different capacities. Each instance is named “ $\textit{num}_{1}-\textit{num}_{2}-\textit{num}_{3}$ ”, where, $\textit{num}_{1}$ represents the number of customers (50, 100, 250), $\textit{num}_{2}$ represents the index of different types of vehicle capacity (large, medium, small), and $\textit{num}_{3}$ represents the index of time window configuration (0, 1, 2, 3, 4). This data set and configuration file can be downloaded from [48].

Comparison algorithms: Since MMVRPTW is a new problem and no existing results can be directly used for comparison, we compared INSGA-II [49] (Improved NSGA-II) (equivalent to the third stage of the 3S-MMDEA) and HSA-HGBS [39] (equivalent to the second and third stages of the 3S-MMDEA). All algorithms are carried out under equal conditions. Equal conditions mean using the same starting and termination criterion, where the same starting criterion means the same primitive population, in the first stage of the three-stage algorithm, the population is initialized by using GOBL, while other comparative algorithms directly use their respective initialization methods to process the original population, equal number of starting search points, the same data set, the same hardware running the algorithms.

Evaluation indicators: A single performance indicator cannot comprehensively measure the performance of the multi-objective optimization algorithm. Therefore, we used four metrics, namely, Inverted generation distance (IGD) [50], 1/HV (HV is super capacity, hyper volume [51]), 1/PSP (PSP is Pareto Sets Proximity [34], it reflects the similarity between the obtained PSs and the true PSs.) and Number of Solution (NOS) [52]. Among them, 1/HV and 1/PSP, the most commonly used performance indicators for IGD multi-objective optimization problems, can measure performance in decision space and objective space respectively. The smaller the value of 1/HV and 1/PSP, the better the performance. They are commonly used evaluation indicators for MMOP.

In addition, the purpose of solving MMVRPTW is to find all Pareto optimal solutions that may have the same target value. Therefore, the number of different optimal solutions (NOS) is used to evaluate the results: NOS represents the number of solution sets $x^{*}$ that satisfy both $x^{*}\in S^{*}$ and $x^{*}\in S$ , where $S^{*}$ represents the real Pareto optimal solution set, and $S$ represents the obtained solution set.

For IGD, 1/HV and 1/PSP, we used Wilcoxon signed rank test at 5% significance level. ‘B/S/W’ indicates that the effect of the proposed algorithm is significantly better/basically similar/significantly inferior to the current algorithm. In addition, the final ranking of all algorithms on the instance set is given using Friedman test [11]. The statistical result of the instance set is summarized as ‘B/S/W ’, which means that the effectiveness of the proposed algorithm is significantly better/basically similar/significantly inferior to the current comparison algorithm on B/S/W instances.

4.2 Experimental results and analysis

Table 1 shows the results of three algorithms on IGD, 1/HV and 1/PSP based on actual examples. Table 2 shows the performance comparison statistics results of all algorithms under the actual example. Table 3 shows the sorting results of Friedman test of the running results of the actual example. All the above results are the average of 30 times of algorithm execution. The B/S/W results in Table 2 show that for 1/HV and IGD, 3S-MMDEA is significantly better than HSA-HGBS in about half of the instances and INSGA-II in most instances. In Wilcoxon signed rank test, the R+value obtained by 3S-MMDEA is higher than the R-value in all cases. This means that 3S-MMDEA is superior to other comparison algorithms in all cases. In addition, the $p$ value is less than 0.05 in all cases, indicating that 3S-MMDEA is obviously superior to other comparison algorithms. It can be seen from Table 3 that under the three evaluation indicators, 3S-MMDEA ranks first, followed by HSA-HGBS and INSGA-II.

Table 1
Average values of IGD, 1/HV, and 1/PSP of 3S-MMDEA INSGA-II, HSA-HGBS

Instance	IGD			1/HV			1/PSP
	INSGA- II	HSA- HGBS	3S- MMDEA	INSGA- II	HSA- HGBS	3S- MMDEA	INSGA- II	HSA- HGBS	3S- MMDEA
50-0-0	0.001424	0.001305	0.001337	1.253734	1.248022	1.245421	0.051083	0.051125	0.052085
50-0-1	0.004022	0.004017	0.004019	2.305209	2.432974	2.243053	0.051562	0.052531	0.053652
50-0-2	0.005609	0.004839	0.005203	2.681102	2.546863	2.546051	0.050431	0.047782	0.033601
50-0-3	0.002023	0.002362	0.002517	1.607124	1.615102	1.616053	0.050342	0.047027	0.034695
50-0-4	0.005032	0.005304	0.004825	1.065537	1.065969	1.066763	0.048761	0.043652	0.037127
50-1-0	0.001981	0.001364	0.000964	1.220176	1.215252	1.173075	0.094013	0.085508	0.085371
50-1-1	0.004552	0.003843	0.003016	2.304409	2.216015	2.213072	0.084637	0.077642	0.068023
50-1-2	0.005776	0.005503	0.004532	2.733521	2.590561	2.497621	0.091201	0.086521	0.066292
50-1-3	0.002728	0.002237	0.001342	1.611206	1.586207	1.524507	0.089052	0.085562	0.076525
50-1-4	0.005637	0.005261	0.003739	1.073185	1.077122	1.079531	0.086012	0.083082	0.085072
50-2-0	0.006247	0.006116	0.006218	1.114207	1.109421	1.109431	0.097316	0.087716	0.052876
50-2-1	0.004982	0.004026	0.003732	1.409902	1.385012	1.360542	0.097001	0.096521	0.096087
50-2-2	0.004235	0.003721	0.004153	1.623906	1.568215	1.125895	0.098862	0.086642	0.073695
50-2-3	0.004127	0.003016	0.003148	1.372018	1.351548	1.420634	0.096201	0.093082	0.094116
50-2-4	0.009391	0.007623	0.006854	1.146036	1.137208	1.121108	0.091341	0.091762	0.091802
150-0-0	0.001902	0.001003	0.001237	1.315102	1.300894	1.300658	0.063651	0.059902	0.061498
150-0-1	0.006106	0.005319	0.005896	3.334686	3.348148	3.349852	0.062983	0.057652	0.057232
150-0-2	0.007208	0.005902	0.005885	3.906072	3.784509	3.339864	0.063019	0.062838	0.063082
150-0-3	0.002902	0.001874	0.001772	2.816015	2.745207	2.810984	0.061084	0.059987	0.060411
150-0-4	0.008025	0.007852	0.007906	2.513307	2.493121	2.466726	0.066391	0.059027	0.053114
150-1-0	0.001292	0.001201	0.001126	1.323201	1.298532	1.253861	0.078826	0.074872	0.065792
150-1-1	0.006403	0.005276	0.004865	3.401895	3.397461	3.401086	0.078211	0.067762	0.055748
150-1-2	0.007336	0.006017	0.005529	3.816265	3.751283	3.374027	0.078621	0.071124	0.067085
150-1-3	0.002936	0.002015	0.001437	2.816114	2.779921	2.811072	0.079982	0.072308	0.070127
150-1-4	0.008942	0.007725	0.007428	2.512085	2.493428	2.390524	0.078643	0.072487	0.071648
150-2-0	0.002484	0.001663	0.001238	1.217989	1.190409	1.180103	0.087652	0.078752	0.075815
150-2-1	0.005803	0.005797	0.005772	2.672107	2.614605	2.496652	0.086431	0.076673	0.075932
150-2-2	0.007538	0.006328	0.005986	4.055982	4.006213	3.986417	0.088735	0.078647	0.076561
150-2-3	0.002922	0.001985	0.001325	2.764018	2.596109	2.409721	0.088027	0.081098	0.081391
150-2-4	0.008936	0.007982	0.007341	2.025998	1.999218	1.846542	0.083772	0.082863	0.082752
250-0-0	0.002721	0.001997	0.001583	1.360986	1.354523	1.320015	0.102673	0.099076	0.098562
250-0-1	0.006834	0.005969	0.005875	4.082013	3.590917	3.604087	0.111092	0.107208	0.107285
250-0-2	0.007132	0.006762	0.007021	4.582458	4.271026	3.743086	0.111175	0.106725	0.093025
250-0-3	0.005423	0.004382	0.004216	1.775679	1.712576	1.698084	0.113872	0.108826	0.093531
250-0-4	0.010639	0.009218	0.008542	2.108001	2.079207	2.076207	0.116018	0.109117	0.094072
250-1-0	0.002775	0.002108	0.001108	1.360862	1.354421	1.327632	0.141102	0.126603	0.109615
250-1-1	0.006861	0.006765	0.006015	4.087982	3.999254	3.865321	0.143652	0.129062	0.110245
250-1-2	0.007903	0.006662	0.006006	4.552961	4.518743	4.550972	0.146291	0.132086	0.109082
250-1-3	0.005228	0.004019	0.003548	1.776011	1.763329	1.769806	0.145309	0.132001	0.109627
250-1-4	0.010628	0.009017	0.008542	2.119821	2.045202	2.024095	0.144421	0.133197	0.092082
250-2-0	0.002427	0.001782	0.001531	1.276972	1.266318	1.243806	0.133872	0.129252	0.091453
250-2-1	0.007207	0.005721	0.006012	4.312631	4.148561	3.843371	0.142798	0.131187	0.090112
250-2-2	0.007569	0.006527	0.006365	4.804902	4.702237	4.606027	0.135672	0.127321	0.092013
250-2-3	0.005204	0.004118	0.003652	1.774872	1.766523	1.695863	0.142088	0.131129	0.092107
250-2-4	0.012404	0.010107	0.009652	2.329787	2.318302	2.280753	0.137872	0.126635	0.092208

Table 2

Statistics performance comparisons on practical instances

	Indicators	B/S/W	R $+$	R $-$	$p$ -value
3S-MMDEA VS INSGA-II	IGD	35/9/1	3098.1	872.6	8.08121e-010
	1/HV	39/4/2	3109.4	896.5	6.79529e-008
	1/PSP	39/3/3	2899.7	1008.6	7.95421e-009
3S-MMDEA VS HSA-HGBS	IGD	27/11/7	3682.2	289.7	2.12937e-012
	1/HV	29/12/4	3597.5	424.6	2.54417e-010
	1/PSP	30/10/5	3318.6	398.5	3.09125e-011

It can be seen from Tables 1 and 2 that the difficulty of the problem increases with the increase of the number of customers and the decrease of the vehicle capacity. The reason is that the problem of having more customers and smaller vehicles will have more path planning solutions, so it will be more difficult to converge on all goals.

Table 3

Average ranking of the algorithms by Friedman test on practical instances

Indicators		3S-MMDEA	HSA-HGBS	INSGA_II
IGD	Average ranking value	1.20643	1.45298	2.54235
	Final rank	1	2	3
1/HV	Average ranking value	1.18762	1.45872	2.20871
	Final rank	1	2	3
1/PSP	Average ranking value	1.07663	1.54214	2.09764
	Final rank	1	2	3

In the real world, preference for one target (objective) may be higher than that for other targets. In the MMVRPTW with five targets considered in this paper, from a certain point of view, the target $f_{2}$ (total driving distance) may be more important than other targets, because the driving distance is proportional to fuel consumption, so it has a direct impact on environmental pollution. In addition, for logistics companies, the target $f_{5}$ (delay time) may be more important than other targets, because service delay will affect customer satisfaction. Therefore, it is very important for MMVRPTW to find the optimal value among all objectives by one method, and the preference of decision-makers is a priori unknown, so it is also very important to obtain multiple equivalent optimal paths.

Table 4 shows the average number of NOS obtained by all algorithms. It can be seen from Table 4 that with an increase in the number of customers (from 50 to 150 to 250), the number of solutions increases. In all instances, the 3S-MMDEA can stably obtain all different solutions. In addition, compared with the other two algorithms, the number of solution sets obtained by the three-stage strategy is obviously better.

In order to further verify the efficiency of the multi-mode multi-target algorithm proposed in this paper, we compare the method proposed in this paper with the MMOP algorithm DN-NSGAII [30], Omni optimizer algorithm [32], MO_PSO_MM [35], MMODEA_ICD [42] made comparison on “50-0-3”, “50-1-2”, “150-0-1”, “150-2-4”, “250-1-1” and “250-2-3” (randomly selected). Similarly, all algorithms are carried out under equal conditions. Equal conditions mean using the same starting and termination criterion, equal number of starting search points, the same data set, the same hardware running the algorithms. The comparison results between different MMO algorithms are shown in Table 5. It can be seen from Table 5 that 3S-MMDEA shows good effect on all evaluation indicators.

Table 4

Average number of solutions obtained by 3S-MMDEA and other algorithms in 30 experiments for each problem instance

Instance	INSGA-II	HSA-HGBS	3S-MMDEA	Instance	INSGA-II	HSA-HGBS	3S-MMDEA
50-0-0	1.2	1.2	3.2	150-1-3	2.5	3.6	6.8
50-0-1	1.2	1.4	3.4	150-1-4	2.6	3.7	6.7
50-0-2	1.1	1	3.5	150-2-0	2.4	3.5	7.1
50-0-3	1	1	3	150-2-1	2.5	3.8	6.6
50-0-4	1	1.3	3.3	150-2-2	2.9	3.9	6.8
50-1-0	1	1	3.2	150-2-3	2.8	3.6	7.2
50-1-1	1	1	3.4	150-2-4	2.7	3.6	7.2
50-1-2	1.2	1	3	250-0-0	4.6	5.4	9.5
50-1-3	1	1	3.5	250-0-1	4.8	5.5	9.6
50-1-4	1	1.3	3.1	250-0-2	4.5	5.3	9.7
50-2-0	1	1	3.4	250-0-3	4.7	5.2	9.4
50-2-1	1.3	1	3	250-0-4	4.5	5.4	9.7
50-2-2	1	1.1	3.2	250-1-0	4.7	5.6	9.4
50-2-3	1.2	1	3.4	250-1-1	4.8	5.2	9.6
50-2-4	1	1	3	250-1-2	4.7	5.3	9.5
150-0-0	2.4	3.8	6.9	250-1-3	4.6	5.2	9.4
150-0-1	2.5	3.7	7.2	250-1-4	4.5	5.3	9.6
150-0-2	2.4	3.7	6.8	250-2-0	4.6	5.4	9.6
150-0-3	2.6	3.9	7.2	250-2-1	4.7	5.3	9.7
150-0-4	2.5	3.6	7.0	250-2-2	4.7	5.2	9.5
150-1-0	2.4	3.8	7.1	250-2-3	4.9	5.5	9.5
150-1-1	2.5	3.9	7.2	250-2-4	4.8	5.4	9.4
150-1-2	2.7	3.9	7.1	$\backslash$	$\backslash$	$\backslash$	$\backslash$

Table 5

Comparison results of different multi-modal multi-objective algorithms

Indicators	Algorithms	50-0-3	50-1-2	150-0-1	150-2-4	250-1-1	250-2-3
IGD	DN-NSGAII	0.004765	0.008625	0.009852	0.010026	0.009754	0.007652
	Omni-optimizer	0.005673	0.009317	0.010265	0.009821	0.010652	0.008725
	MO_PSO_MM	0.004216	0.007625	0.008524	0.009255	0.008639	0.006251
	MMODEA_ICD	0.003529	0.005428	0.006552	0.008728	0.007126	0.004237
	3S-MMDEA	0.002517	0.004532	0.005896	0.007341	0.006015	0.003652
1/HV	DN-NSGAII	2.014216	3.766583	4.220867	2.762296	4.773815	2.963926
	Omni-optimizer	1.997622	3.976215	4.390115	2.983148	4.826419	3.002641
	MO_PSO_MM	1.872534	3.524192	4.012972	2.659721	4.562164	2.873165
	MMODEA_ICD	1.762538	3.029721	3.875214	2.012862	4.109273	2.39175
	3S-MMDEA	1.616053	2.497621	3.349852	1.846542	3.865321	1.695863
1/PSP	DN-NSGAII	0.065213	0.089216	0.082017	0.091752	0.162962	0.102681
	Omni-optimizer	0.702645	0.100286	0.083865	0.098621	0.163786	0.120172
	MO_PSO_MM	0.629721	0.078271	0.072011	0.091027	0.156522	0.098621
	MMODEA_ICD	0.047652	0.072019	0.062917	0.086512	0.130972	0.093826
	3S-MMDEA	0.034695	0.066292	0.057232	0.082752	0.110245	0.092107
NOS	DN-NSGAII	1.2	1.4	4.2	4.5	7.2	7.4
	Omni-optimizer	1.4	1.4	4.6	4.5	7.5	7.7
	MO_PSO_MM	2.1	2.2	4.6	4.6	7.7	7.8
	MMODEA_ICD	2.7	2.8	6.5	6.2	8.2	8.1
	3S-MMDEA	3	3	7.2	7.2	9.6	9.5

To better illustrate the performance of the proposed 3S-MMDEA, considering factors such as the actual logistics distribution constraints and the applicable solution scale of the compared algorithms, taking “50-1-2” as an example, which has a certain representativeness, we give the distribution of 3S-MMDEA and HSA-HGBS, MMODEA_ICD and INSGA-II in the final stage objective space. In multi-modal multi-objective optimization, the solution obtained when the objective function exceeds three is a hypersurface, which cannot be reflected in the coordinate system.

This article considers the solution of $f_{1}-f_{5}$ objective functions, namely the minimum number of vehicles, the shortest driving distance, the minimum makespan, the minimum total waiting time, and the minimum total delay time. When mapping these five objective functions to the solution space, due to the limitations of experimental simulation, only three-dimensional coordinates can be displayed. Among the five objective functions, $f_{3}$ , $f_{4}$ and $f_{5}$ have a certain correlation, and when compared with [39, 49], $f_{5}$ can represent $f_{3}$ and $f_{4}$ to a certain extent. Therefore, take $f_{1}$ , $f_{2}$ and $f_{5}$ as the three dimensions displayed in the solution space, and the weight is set at a ratio of 1:1:1, that is, process $f_{1}$ , $f_{2}$ and $f_{5}$ in a weighted sum manner. The distribution of solutions obtained by the three objective functions in the objective space is shown in Fig. 6.

Figure 6.

Distribution of solution sets in the objective space obtained by different algorithms.

It can be seen from Fig. 6 that for the 3S-MMDEA, the optimal solution set is uniformly distributed in the objective space, showing a great diversity of solution sets, and the final solution set converges to a small region. The convergence and diversity of the other three algorithms are poor, and even some algorithms have lost some dominant solutions.

To sum up all the experiments and results, it can be seen that the 3S-MMDEA can maintain all the optimal solutions even if it has the same target value, while ensuring the diversity and convergence of the final solution set. For other MMO algorithms, it is difficult to ensure the diversity and convergence of objective space and decision space at the same time. For the traditional multi-objective optimization algorithm, although it is easy to find the Pareto front, it is difficult to find all equivalent solutions.

5. Conclusion

Vehicle routing planning is a practical problem aimed at finding the best route and realizing efficient logistics distribution. Since the preference of decision-makers is a priori unknown, it is necessary to provide multiple alternative solutions, that is, to maintain the diversity of solutions in the decision space and the convergence quality in the target space. This paper regards VRPTW as an MMOP, to solve this problem, a three-stage multi-mode multi-objective differential evolution algorithm is proposed to solve VRPTW. In the first stage, GOBL strategy is used to generate the reverse population to expand the selection range of individuals. In the second stage, a reachable distribution area division strategy is proposed to realize the population division with the selected individuals as the center point to improve the convergence of the solution set. In the third stage, the diversity of the solution set is further improved by improving the individual variation strategy to legalize the mutation individuals.

In this paper, the effectiveness of 3S-MMDEA have verified through logistics distribution of different scales in the standard data-set. By comparing with INSGA-II and HSA-HGBS, the experimental results show that the 3S-MMDEA proposed in this paper can efficiently solve VRPTW of different scales, and can obtain more equivalent optimal paths. The algorithm performs better in convergence and diversity. In addition, by comparing with DN-NSGAII, Omni-optimizer, MO_PSO_MM, MMODEA_ICD and other multi-modal multi-objective algorithms on the standard data-set, the experimental results show that 3S-MMDEA can obtain multiple equivalent optimal solutions while ensuring convergence.

Although the three-stage algorithm proposed in this paper can better optimize the vehicle routing problem with time window constraints, there are still some limitations. 3S-MMDEA is not an accurate algorithm, and in future research, we should consider the certainty of probability in optimization, and combine the characteristics of intelligent algorithm probability optimization to improve runtime and optimization efficiency. What’s more, in the subsequent research, it is still necessary to consider the local search in line with specific objectives and improve the performance of the algorithm.

Footnotes

Acknowledgments

This work is supported by the National Natural Science Foundation of China under Grant 61806006, China Postdoctoral Science Foundation under Grant No. 2019M660149, Graduate Innovation Foundation of Jiangsu Province under Grant No. KYLX16_0781, the 111 Project under Grants No. B12018, and PAPD of Jiangsu Higher Education Institutions.

References

Deng

et al., Multi-objective Path Optimization Method in Terminal Building Based on Improved Genetic Algorithm, in: 2020 Chinese Automation Congress (CAC), 2020, pp. 3181–318.

Miao

et al., Path planning optimization of indoor mobile robot based on adaptive ant colony algorithm, Computers & Industrial Engineering 156 (2021), 107230.

Yue

C.T.

et al., A novel scalable test problem suite for multi-modal multi-objective optimization, Swarm and Evolutionary Computation 48 (2019), 62–71.

Liu

Yen

G.G.

and Gong

, A multi-modal multi-objective evolutionary algorithm using two-archive and recombination strategies, IEEE Transactions on Evolutionary Computation 23(4) (2019), 660–674.

Liang

J.J.

et al., Problem Definitions and Evaluation Criteria for the CEC 2021 on Multi-modal Multi-objective Path Planning Optimization, Computational Intelligence Laboratory, Zhengzhou University, Zhengzhou China and Technical Report, Nanyang Technological University, Singapore, Krakow, Poland, June, 2021.

Wang

et al., A novel multi-objective competitive swarm optimization algorithm for multi-modal multi objective problems, in: 2019 IEEE Congress on Evolutionary Computation (CEC), IEEE, 2019, pp. 271–278.

Solomon

M.M.

, Algorithms for the vehicle routing and scheduling problem with time window constraints, Operations Research 35(2) (1987), 254–265.

S.H.

et al., A combination of genetic algorithm and particle swarm optimization for vehicle routing problem with time windows, Sensors 15(9) (2015), 21033–21053.

Nalepa

and Blocho

, Adaptive memetic algorithm for minimizing distance in the vehicle routing problem with time windows, Soft Computing 20(6) (2016), 2309–2327.

10.

Jose

et al., An ACS-based memetic algorithm for the heterogeneous vehicle routing problem with time windows, Expert Systems With Applications 157 (2020).

11.

Wang

Weng

and Zhang

, A Two-Stage Multi-objective Evolutionary Algorithm for Multi-objective Multi-depot Vehicle Routing Problem With Time Windows, IEEE Transactions on Cybernetics, 2018; 1–12.

12.

et al., Multi-depot multi-trip vehicle routing problem with time windows and release dates, Transportation Research Part E 135 (2020), 1–21.

13.

Cueto

P.F.

et al., A solution approach for multi-trip vehicle routing problems with time windows, fleet sizing, and depot location, Networks 78(4) (2021), 503–522.

14.

Cheng

Cai

and Fu

, A Preliminary Study of Evolutionary Multitasking for Multi-objective Vehicle Routing Problem With Time Windows, in: The 2021 6th International Conference on Computational Intelligence and Applications (ICCIA), IEEE, 2021.

15.

Shu

et al., Two-stage multi-objective evolutionary algorithm based on classified population for tri-objective VRPTW, International Journal of Unconventional Computing 16 (2021).

16.

Hou

and Han

, Solution Evaluation-Oriented Multi-objective Differential Evolution Algorithm for MOVRPTW, in: The 2021 8th International Conference on Information, Cybernetics, and Computational Social Systems (ICCSS), Piscataway: IEEE, 2021, pp. 50–55.

17.

Cai

et al., A hybrid evolutionary multitask algorithm for the multi-objective vehicle routing problem with time windows, Information Sciences, 2022.

18.

Srivastava

and Singh

, Two evolutionary approaches with objective-specific variation operators for vehicle routing problem with time windows and quality of service objectives, Applied Soft Computing, 2023.

19.

Zhang

H.F.

et al., Review of vehicle routing problems: Models, classification and solving algorithms, Archives Computational Methods Engineering 29(1) (2022), 195–221.

20.

Akyo

and Alatas

, Plant intelligence based meta heuristic optimization algorithms, Artificial Intelligence Review 47(4) (2016).

21.

Alatas

and Bingol

, Comparative assessment of light-based intelligent search and optimization algorithms, Light & Engineering 3 (2020), 51–59.

22.

Alatas

and Bingol

, A physics based novel approach for travelling tournament problem: Optics inspired optimization, Information Technology and Control 3 (2019).

23.

Bingol

and Alatas

, Chaos based optics inspired optimization algorithms as global solution search approach, Chaos, Solitons & Fractals 141 (2020).

24.

Heßler

and Irnich

, Abranch-and-cut algorithm for the soft-clustered vehicle-routing problem, Discrete Applied Mathematics 288 (2021), 218–234.

25.

Lahyani

Khemakhem

and Semet

, Rich vehicle routing problems: From a taxonomy to a definition European, Journal of Operational Research 241(1) (2015), 1–14.

26.

Zhou

et al., Multi-objective evolutionary algorithms: A survey of the state of the art, Swarm and Evolutionary Computation 1(1) (2011), 32–49.

27.

Deb

Pratap

Agarwal

et al., A fast and elitist multi-objective genetic algorithm: NSGA-II, IEEE Transactions on Evolutionary Computation 6(2) (2002), 182–197.

28.

Yue

C.T.

et al., A novel scalable test problem suite for multi-modal multi-objective optimization, Swarm and Evolutionary Computation 48 (2019), 62–71.

29.

Liu

Yen

G.G.

and Gong

, A multi-modal multi-objective evolutionary algorithm using two-archive and recombination strategies, IEEE Transactions on Evolutionary Computation 23(4) (2019), 660–674.

30.

Liang

J.J.

Yue

C.T.

and Qu

B.Y.

, Multi-modal multi-objective optimization: A preliminary study, in: Proceedings of the IEEE Congress on Evolutionary Computation, Vancouver: IEEE, 2016, pp. 2454–2461.

31.

et al., A self-organizing multi-modal multi-objective pigeon-inspired optimization algorithm, Science China Information Sciences 62(7) (2019), 1–17.

32.

Deb

and Tiwari

, Omni-optimizer: A procedure for single and multi-objective optimization, in: Proceedings of the International Conference on Evolutionary Multi-Criterion Optimization, Berlin: Springer, 2005, pp. 47–61.

33.

Yan

et al., A performance enhanced niching multi-objective bat algorithm for multi-modal multi-objective problems, in: Proceedings of the IEEE Congress on Evolutionary Computation, Wellington: IEEE, 2019, pp. 1275–1282.

34.

Yue

C.T.

B.Y.

and Liang

J.J.

, A multi-objective particle swarm optimizer using ring topology for solving multi-modal multi-objective problems, IEEE Transactions on Evolutionary Computation 22(5) (2018), 805–817.

35.

Liang

J.J.

et al., A self-organizing multi-objective particle swarm optimization algorithm for multi-modal multi-objective problems, in: Proceedings of the International Conference on Swarm Intelligence, Cham: Springer, 2018, pp. 550–560.

36.

Storn

and Price

, Differential evolution: A simple and efficient adaptive scheme for global optimization over continuous spaces, Journal of Global Optimization 11(4) (1997), 341–359.

37.

Gupta

and Su

, Multiple individual guided differential evolution with time varying and feedback information-based control parameters, Knowledge-based systems, 2023.

38.

Petrovi

et al., Controlling highway toll stations using deep learning, queuing theory, and differential evolution, Engineering Applications of Artificial Intelligence 119 (2023), 105683.

39.

Sun

Zhang

H.F.

and Dang

J.W.

, Two-stage vehicle routing optimization for logistics distribution based on HSA-HGBS algorithm, IEEE Access 10 (2022), 99646–99660.

40.

Zhang

H.F.

et al., Combining Affinity Propagation with Differential Evolution for Three-echelon logistics distribution optimization, Applied Soft Computing 131C(109878) (2022).

41.

Wang

et al., Space transformation search: a new evolutionary technique, in: Genetic & Evolutionary Computation Conference, DBLP, 2009, pp. 537–544.

42.

Yue

C.T.

et al., Differential evolution using improved crowding distance for multi-modal multi-objective optimization, Swarm and Evolutionary Computation 62(9) (2021), 100849.

43.

Zitzler

et al., Performance assessment of multi-objective optimizers: An analysis and review, IEEE Transactions on Evolutionary Computation 7(2) (2003), 117–132.

44.

Liang

J.J.

et al., Multi-modal multi-objective optimization with differential evolution, Swarm and Evolutionary Computation 44 (2019), 1028–1059.

45.

Yong

Cai

and Zhang

, Differential evolution with composite trial vector generation strategies and control parameters, IEEE Transactions on Evolutionary Computation 15(1) (2011), 55–66.

46.

Zhou

and Wang

, A local search-based multi-objective optimization algorithm for multi-objective vehicle routing problem with time windows, IEEE Systems Journal 9(3) (2017), 1100–1113.

47.

Gutiérrez

J.P.C.

Landa-Silva

and Moreno-Pérez

J.A.

, Nature of real-world multi-objective vehicle routing with evolutionary algorithms, in: IEEE International Conference on Systems, 2011, pp. 257–264.

48.

Castro-Gutierrez

Landa-Silva

and Moreno Perez

J.A.

, MOVRPTW data set, 2010. https://github.com/psxjpc/.

49.

Srivastava

Singh

and Mallipeddi

, NSGA-II with objective-specific variation operators for multi-objective vehicle routing problem with time windows, Expert Systems with Applications 176(1) (2021), 114779.

50.

Zhang

and Li

, MOEA/D: A multi-objective evolutionary algorithm based on decomposition, IEEE Transactions on Evolutionary Computation 11(6) (2008), 712–731.

51.

Zitzler

et al., Performance assessment of multi-objective optimizers: an analysis and review, IEEE Transactions on Evolutionary Computation 7(2) (2003), 117–132.

52.

Yao

et al., Multi-modal multi-objective evolutionary algorithm for multiple path planning, Computers & Industrial Engineering 169(Pt2) (2022).

Three-stage multi-modal multi-objective differential evolution algorithm for vehicle routing problem with time windows

Abstract

Keywords

1. Introduction

2.1 Relevant theories

(1) Initialization

(2) Variation

(3) Crossover

(4) Selection

3.1 Basic framework of 3S-MMDEA

3.2 Stage I: GOBL for population initialization

3.3 Stage II: Search strategy of reachable distribution area for population division

3.4 Stage III: Improved individual variation strategy

4. Experimental results and analysis

4.1 Experimental settings

4.2 Experimental results and analysis

Table 1 Average values of IGD, 1/HV, and 1/PSP of 3S-MMDEA INSGA-II, HSA-HGBS

Footnotes

Acknowledgments

References

Table 1
Average values of IGD, 1/HV, and 1/PSP of 3S-MMDEA INSGA-II, HSA-HGBS