Nondominated Sorting Differential Evolution Algorithm to Solve the Biobjective Multi-AGV Routing Problem in Hazardous Chemicals Warehouse

For the multiple automated guided vehicle (multi-AGV) routing problems in the warehousing link of logistics, where the optimization objective is to minimize both the number of AGVs used and the maximum pickup time simultaneously, a nondominant sorting differential evolution (NSDE) algorithm is proposed. In the encoding and decoding stages, the pickup point area is divided. AGVs are allocated to each region according to the proposed rule based on avoiding duplicate paths. Meanwhile, the pickup points within the region can be adjusted to optimize the pickup paths and improve the pickup efficiency. The fast nondominated sorting method and elitist selection strategy in the nondominated sorting genetic algorithm II (NSGA-II) are introduced into the differential evolution algorithm, which sorts all the regions to obtain the best Pareto solution set. Lastly, the domination of the proposed NSDE algorithm in Pareto frontier evaluation indicators is verified by some numerical experiments.


Introduction
Dangerous chemical accidents have occurred frequently in many countries. Dangerous chemicals have the properties of corrosion, toxicity, and harmfulness, and are easy to explode and combust, so to do a good job in the safety management of the storage, loading and unloading of dangerous chemicals is an important measure to avoid safety accidents. To achieve safety and e ciency in material handling, automated guided vehicle (AGV) system is widely applied in warehouses, manufacturing systems, container terminals, and service industries including hospitals [1]. As a exible and intelligent warehousing logistics equipment integrating a variety of advanced technologies, AGV has the characteristics of automation, high exibility, high e ciency, high reliability, and parallel operation [2], which can well meet this logistics demand of dangerous chemicals. erefore, more and more automatic warehouses of dangerous chemicals use AGV for transportation. is type of transportation problem can be classi ed as multi-AGV routing problem in an automated warehouse: multiple AGVs start from their respective entrances, drive along passageways pickup the goods at the pickup point mentioned in the order, and nally carry the goods to the exit. AGV routing involves the operation sequence and driving path arrangement of AGV [3]. At present, metaheuristic algorithm is still one of the major ways to solve this kind of complex combinatorial optimization problem, and many scholars have carried out relevant research under various warehouse transportation backgrounds.
AGV handling system has been an essential component of modern manufacturing systems. Kabir and Suzuki [4] conducted a comparative analysis of four AGV routing heuristics for battery management: select the nearest battery station, select minimum delay battery station, select the nearest battery station from the initial point, and select the nearest battery station to drop-off point. In terms of system productivity, selecting minimum delay battery station performs the best and selecting the nearest battery station from the initial point performs the worst. Considering energy consumption or its environmental impact indicators, Zhang et al. [5] formulated an energy-efficient path planning model for a single-load AGV and selected two optimization objectives, transport distance, and energy consumption. ey put forward a particle swarm optimization-based method where an indirect priority-based particle encoding scheme is proposed. For the multi-compartment AGV scheduling problem in matrix workshop, Zou et al. [6] addressed an effective iterated greedy algorithm with some advanced techniques including accelerations for evaluating neighborhood solutions, two improved constructive heuristics based on nearest-neighbour and sweep, an improved destruction procedure, and simulated annealing type of acceptance criterion. Zou et al. [7] then investigated a new AGV scheduling problem with pickup and delivery that simultaneously optimized two objectives of maximizing the overall customers' satisfaction and minimizing the total distribution costs. To solve the problem, an effective multiobjective evolutionary algorithm is proposed, containing a constructive heuristic in population initialization, a multiobjective local search based on an ideal point to improve the capability of exploitation, a two-point crossover operation to exploit helpful information collected in the nondominated solutions, and a restart strategy to prevent trapping into a local optimum. For the problem of collision and deadlock, Zając and Małopolski [8] presented an algorithm introducing a structural online control strategy, which has a considerably high efficiency in three types of road systems: unidirectional, bidirectional, and mixed. For the cell formation and intercellular integrated scheduling problems in the cellular manufacturing system where AGVs are used for transferring exceptional parts, Goli et al. [9] developed a fuzzy mixed-integer linear programming, designing an efficient whale Optimization Algorithm. Drótos [10] designed a centralized approach to maintain a central schedule to guarantee conflict-free operation of a large fleet of vehicles without deadlock or live lock, proposing a method based on local search to support the optimization of plans. For the multi-AGVs dispatching problem in a matrix manufacturing workshop, Zhang et al. [11] provided an improved iterated greedy algorithm, where an AGV route merging strategy and a workshop partition strategy are proposed to reduce travelling distance and the cost of AGVs, two rules to identify infeasible solutions are presented to save the operation time, four effective operators in the local search stage are applied to improve the solution quality, and a repair strategy is used to avoid the local optima. For the fleet sizing and routing problem with synchronization for AGVs with dynamic demands, Aziez et al. [12] used mathematical programming and a powerful matheuristic algorithm to cope with large instances and real-time operations, handling dynamic demand by replanning the routes immediately with the arrival of new requests.
In the automated container terminal, reasonable AGV routing is beneficial to the smoothness of the whole logistics transport. For the integrated scheduling problem of handling equipment and AGV, Yang et al. [13] set up a bilevel programming model to minimize the makespan, where AGV path planning is the lower level, and the integrated scheduling of quay cranes, AGVs, and automatic railmounted gantries is the upper level. ey designed an effective bilevel general algorithm based on the preventive congestion rule. Aiming at the goal of minimizing multi-AGV delays, Zhong et al. [14] considered conflict-free path planning in simultaneous scheduling and presented a hybrid genetic algorithm-particle swarm optimization algorithm with fuzzy logic controller to adaptive auto-tuning which is more reliable than the genetic algorithm, especially on largescale problems. Ji et al. [15] investigated the combinatorial optimization of two problems in the synchronous loading and unloading operation pattern, and established a programming model aiming at the layout and operation flow of most container terminals to minimize the completion time of all containers. ey proposed two bilevel genetic algorithms based on conflict resolution strategy, introducing elitism preservation strategy, tabu list, and catastrophe operator. Zhou et al. [16] developed an improved anisotropic Q-learning routing algorithm to find the shortest-time paths in the guide-path network of crosslane type, considering AGV real-time status and designing selecting-action strategy based on AGV waiting time estimation.
Similar to manufacturing systems and container terminals, the application of AGVs in automated warehouses has greatly improved the efficiency in handling goods. Xing et al. [17] designed two conflict elimination methods and presented a novel tabu search algorithm. e algorithm adopts two novel neighborhood operators including relocation and exchange operations, which greatly optimizes the order of pickup points in the automatic warehouse and shortens the total travel distance. Fransen et al. [18] proposed a dynamic path planning approach to solve the realtime control problem of grid-based AGV systems in parcel sorting, baggage handling, and semiconductor fabrication and adopted a graph-representation of the system with constantly updating vertex weights, which can recover from deadlock situations. To solve the integration difficulty between scheduling and routing aspects of the multi-AGV problem, Lee et al. [19] established a model of combinatorial auction-based winner determination problem for multiple AGVs and decomposed the model into a two-phase problem, proposing a new genetic algorithm with knowledgebased operators to obtain a better solution. Chen et al. [20] set up an Ant-agent optimized by repulsive potential field for multi-AGV routing, designing a scheme of varying velocity to adjust the AGV velocity through decentralized control, and presenting a novel transition rule to determine the AGV motion state by combining centralized and decentralized control.
Most AGV path optimization problems only optimize one objective. However, there are often multiple conflicting optimization objectives existing in practice. When the number of AGVs used is small, the pickup time will be long; when the requirements for pickup speed are high, more AGVs will need to be started. erefore, AGV quantity and maximum pickup time are two conflicting goals, which are also the focus of this paper. A large number of studies have found that the population-based intelligent optimization algorithm can maintain a population of solutions in the iterative process, and has greater advantages in solving multi-objective optimization problems. Differential evolution (DE) algorithm has the advantages of simple structure, easy implementation, fast convergence, and strong robustness and are widely used in engineering practice, such as computer network, mechanical design and robotics, image processing, industrial control, biology, and other fields [21][22][23].
When the classical differential evolution algorithm is applied to solve optimization problems in some complex environments, it exists the phenomenon of low solution accuracy and poor stability. e improvement for the classical differential evolution algorithm mainly includes: adaptive parameter adjustment mechanism, mutation and crossover strategy adjustment technology, population structure design, and hybrid algorithm evolution mechanism. Zhao et al. [24] added the competitive evolutionary mechanism to mutation operation where the strategy to produce superior individuals are employed in population with larger sizes. Yu et al. [25] put forward an adaptive selection mutation method where individuals are selected by their fitness values and constraint violations, improving the exploitation, and maintaining the exploration. Ahmadianfar et al. [26] created a novel adaptive mutation mechanism, introducing the operators of particle swarm optimization algorithm to improve global search capability. Using the generalized Lehmer mean and linear bias reduction, Stanovov et al. [27] controlled the parameter adaptation bias for fitness improvement-based and distance-based DE algorithms. To exploit the evolutionary state of the population, Gupta and Su [28] proposed a mutation strategy that uses a dynamic number of fitness-based leading individuals, creating a novel way to set control parameters based on each individual's evolutionary state during trial generation. A new mapping mechanism of continuous and discrete variables was proposed by Ali et al. [29]. In addition, a repair method based on k-means clustering is designed to enhance the initial population, and an ensemble of mutation strategies is presented to improve the exploration ability. Wang et al. [30] used the accompanying population whose individuals are composed of suboptimal solutions to optimize the mutation strategy and control parameters, which realized the adaptation of the strategy and parameters of the main population. Besides, they utilized reverse individuals to enhanced population diversity in evolution and radial spatial projection technology to optimize the evolution direction. Tan et al. [31] presented an adaptive mutation operator based on fitness landscape where the relationship between three mutation strategies and fitness landscape features are trained by random forest offline, designing an adaptable parameter mechanism based on historical memory and a linear reduction method of population size. Based on five mutation strategies, Deng et al. [32] developed an optimal mutation strategy to improve local search ability while ensuring global search ability. To ensure diversity of solutions and accelerate convergence, they employed the wavelet basis function and normal distribution function to control the scaling factor and crossover rate, respectively. e crossover rate sorting mechanism introduced by Li et al. [33] enables each individual to be assigned a crossover rate value according to their fitness value, allowing good elements to be more inherited. A dynamic population reduction strategy is employed to enhance convergence speed and balance exploration and exploitation. A selection operator with three candidate vectors was proposed by Zeng et al. [34] for the escaping of the local optimal value when the individual is stuck in a state of stagnation. Sun et al. [35] designed the reverse learning mechanism for generating the initial subpopulations to improve the convergence velocity and maintain population diversity, developing a new multipopulation parallel control strategy to maintain the search efficiency in subpopulations. Deng et al. [36] applied the quantum chromosome encoding to enhance the population diversity and quantum rotation to speed up the convergence speed, using the cooperative coevolution framework to change the large-scale and high-dimensional complex optimization problem to multiple low-dimensional subproblems to improve the solution efficiency. Sun et al. [37] applied a finite-horizon Markov decision process to adaptively control the parameter, avoiding the problem of poor convergence caused by random values or empirical measurements. Houssein et al. [38] proposed a hybrid algorithm of slime mould algorithm and DE Algorithm where an adaptive guided mutation method improve the swarm agents' local search, enhance the diversity of population, and prevent premature convergence. e performance improvement of the traditional differential evolution algorithm can be implemented by considering the characteristics of complex problems and designing strategies for new coding structure and decoding scheme, mutation operator, crossover operator, and evolution mechanism. To our knowledge, there are no known studies that apply DE algorithms to solve such biobjective AGV routing problem, therefore this paper presents a nondominant sorting differential evolution (NSDE) algorithm with a well-designed encoding and decoding method, introducing an elitist selection strategy. By discussing the characteristics of automated warehouse layout and AGV driving mode, a new encoding and decoding method based on pickup location is designed, which can effectively improve the decoding quality. Nondominated sorting genetic algorithm II (NSGA-II) is one of the most influential and widely applied multi-objective genetic algorithms, based on controlled elitism concepts. In the iteration process, the elite strategy of NSGA-II is introduced to accelerate population convergence and maintain the quality and Pareto diversity of the population. e remaining sections are organized as follows: First, the mathematical formulation of biobjective multi-AGV routing problem is provided. en, the NSDE algorithm is detailed. Afterwards, the paper introduces the numerical experiments and results. Lastly, a summary is given.

Problem Statement
is section introduces the working environment and process of AGV and the basic concept of multi-objective problem, leading to the proposed biobjective multi-AGV routing problem and mathematical model. Figure 1 shows the layout of a common automated warehouse which is represented by the grid method [39], divided into 10 rows and 15 columns, a total of 150 square grids, and numbered in order. To maximize the use of warehouse space, the passageway is only designed as the width of one grid, and each shelf occupies two columns. e white area and gray area are driving passageways (transportation passageways) of AGVs and shelves, respectively. e entrances of AGVs are concentrated in the first row in the figure. When AGVs enter the warehouse, they are arranged from left to right, and the exit is located in the grid of the lower right corner.

Workflow Description for AGV.
When AGV works, it starts from its own entrance and runs in the passageway between shelves. According to the requirements of the order, AGV takes out the corresponding goods from the shelves and delivers them to the exit. e goods involved in the order are called the pickup point. When picking up goods, AGV requires driving to the grid next to the pickup point. e grid like this is called the pickup position of the pickup point (for example, the pickup position of pickup point 51 is grid 50), and the passageway located is the pickup passageway of the pickup point.

Multi-Objective Optimization Problem.
In multi-objective optimization problem (MOP), multiple objectives conflict with each other, and the performance of one objective is often improved at the expense of the performance of other objectives. In fact, it is impossible for MOP to have a solution that can achieve the optimal state of all objects, and the solution can usually be obtained from a set of non-inferior solutions, i.e., Pareto solution set [40].
For MOP, suppose there are m optimization objectives f i (x), i � 1, 2, . . ., m, and each optimization objective is to minimize. If any two feasible solutions a and b satisfy the following relationship, then solution a dominates solution b, which can be expressed as a > b. Otherwise, there is no dominating relationship between the two solutions as follows: When there are multiple Pareto optimal solutions, if more information about the problem is unknown, it is difficult to decide which solution is more desirable, i.e., all Pareto optimal solutions are considered to be equally important. erefore, for MOP, the most important task is to find as many Pareto optimal solutions as possible, mainly completing the following two tasks: (1) Find a set of solutions as close to Pareto optimal front as possible.
(2) Find a set of solutions that are as different as possible.

Problem Description and Mathematical
Formulation. e biobjective multi-AGV routing problem is described as follows: enough AGVs are available and n goods need to be taken out according to a batch of customer orders. AGVs start from their respective entrances, drive along the passageway, take out the goods involved in the order, and finally transport them to the exit. Collision is not allowed during AGV driving, i.e., any grid can only be occupied by one AGV at the same time. e goal of the problem is to simultaneously minimize two objects, i.e., the number of AGVs used and the maximum pickup time of AGVs, so as to obtain the conflict-free driving path of AGVs.
is problem is based on the following assumptions: (1) e time when the goods is taken out of the shelf and put into the AGV is not considered. (2) Only focus the case that the quality or volume of goods is small enough and will not exceed the AGV capacity, regardless of capacity limitation of AGV. e decision variable is: e MIP model of the problem is provided as follows: n j�1 x 0j1 � 1, k∈M i∈S j∈S e two objectives of optimization are represented by (2): simultaneous minimization of the number of AGVs used and the maximum pickup time; (3) is the objective function of the number of AGV used; (4) is the objective function of maximum pickup time; constraints (4)- (5) ensure that at least one AGV is working; constraint (7) means that entrance and exit need to be arranged for the working AGV, otherwise they will not be arranged; constraints (7)-(10), respectively, illustrate that the following operations are not allowed in the pickup sequence of AGV: return the entrance from the pickup point, return the pickup point from the exit, go directly to the exit from the entrance, and go directly to the entrance from the exit; constraints (11)- (12) indicate that the AGV should leave a node after accessing the node; constraint (14) eliminates subtours; formula (15) is the calculation formula of the pickup time of each AGV; constraints (16) provides decision variable constraint.

Nondominated Sorting Differential
Evolution Algorithm e differential evolution algorithm proposed by Storn and Price [41] is a population-based optimization algorithm, one of the global intelligent optimization algorithms to solve continuous optimization problems. For the biobjective problem proposed, a NSDE algorithm is presented. Firstly, an effective encoding and decoding method is designed to eliminate the influence of repeated paths, and then the fast nondominated sorting method and elite selection strategy of NSGA-II are introduced to improve the selection operation. Excellent individuals are selected from the two generations of the population for iteration, so that the Pareto solution set is continuously updated. e procedure of the NSDE algorithm is as follows: 3.1. Encoding and Decoding. Since the driving passageways in an automated warehouse are limited by the position of shelves, and AGV has no capacity limitation and does not need to deal with the situation of exiting, unloading, and reentering if it is full, a relatively superior solution should satisfy the following two rules to reduce duplicate paths: Rule 1: e driving direction on the horizontal passageway should be consistent. Rule 2: Reversals on the vertical passageway should be as little as possible.
Regardless of conflicts, for instance, the pickup sequence (123, 126, 129, and 132) is better than (123, 129, 126, and 132) and the pickup sequence (55, 70, 85, 100) is superior to (55, 85, 70, 100) in Figure 1. Figures 2 and 3 clearly display the path arrangement and driving directions on the horizontal and vertical passageways for these two examples, respectively. erefore, applying these two rules can reduce the number of times through the same grid, i.e., decrease repeated paths and reduce the time wasted to some extent. A novel coding and decoding method with the rule of avoiding duplicate paths is designed.

Encoding Method.
e pickup points involved in each vertical pickup passageway should be taken out by two AGVs at most. For example, the goods in the upper and lower sections of the vertical passageway can be picked up by two AGVs, respectively. If another AGV is added to work in this passageway, its path would be identical to that of one of the previous two AGVs. For applying in the discrete, the individual in the NSDE algorithm is expressed as a region-based sequence. erefore, the coding method proposed firstly divides the pickup points related to each vertical pickup passageway into upper and lower areas, which are represented by two numbers, respectively. e pickup points corresponding to each number are arranged in increasing order, so the individual length is equal to the number of vertical passageways multiplied by 2.
ere are several exceptional cases that need to be explained: (1) If there is only one pickup point for a certain vertical passageway, it will be assigned to the corresponding number according to the area it is located, and the other number code related to this passageway has no meaning, but only exists for the convenience of calculation.
(2) If the pickup points of a certain vertical passageway are only distributed in the upper half or the lower half, and there are two or more pickup points, in order to make the two digital codes related to this passageway meaningful, the regional number without pickup points contains the nearest pickup point to this region. (3) If the row number of the warehouse is odd and there are some pickup points in the middle row, they can be distributed to the area that these pickup points are closer to, and if the distances are the same or their pickup passageway only involves these pickup points, they can be distributed to any area. (4) If not every vertical passageway contains pickup points, the number corresponding to the passageway without pickup points does not represent any pickup point, but it can play the role of dividing AGV in the decoding stage to be introduced below. Figure 4 indicates the goods numbers of a batch of orders only, and divides the areas. It is assumed that the pickup points are distributed in five vertical passageways, where 1 to 5 represents the pickup points in the top area of each passageway, and 6 to 10 presents the pickup points in the area below. According to the above coding mode, the meaning of each coded number, i.e., the corresponding pickup points can be obtained in Table 1.
e current meaning of coded numbers is the initial state before decoding, which will be adjusted during decoding to optimize the path.

Decoding Method.
In the decoding phase, the region number in an individual should be converted into the pickup sequence of one or more AGVs, and the picking points should be adjusted across regions to make decoded paths better. In order to reduce the repeated paths on the horizontal passageway, AGV will not turn around on the horizontal passageway when arranging the pickup sequence, i.e., rule 1 should be observed under the background of the warehouse layout in this paper. In the following process, the pickup points represented by each digital code are sorted from small to large. When calculating the pickup time, connect the pickup point sequences represented by the numbers. If the flip sequence can make the pickup sequence closer to the previous pickup point, flip it. e following is a detailed introduction to the decoding process.   If the max{C 2 , C 3 } becomes smaller after adjustment, the above adjustment will be adopted; otherwise, it will be restored to the original codes, i.e., the scheme with the minimum difference between the pick times of two AGVs. (4) Adjust entrance: when there is no adjustable scheme in (3), find out the AGV with the longest pickup time: if its entrance is the nearest, stop adjustment and enter (5); otherwise, try to swap other AGV entrance for it: if the maximum pickup time of the two AGVs can be shortened, the shortest solution is selected, and then return (3); otherwise, stop the adjustment and go to (5). An instance is provided as follows: suppose the following results are obtained after step (3) (only the numbers of goods needed for calculation are given): e pickup time of AGV 1 is assumed to the maximum pickup time, i.e., C 1 � max{C 1 , C 2 , C 3 }, then the nearest entrance 2 is tried to arrange for the route of AGV 1, i.e., AGV 1: 2 up (19, 49, and 66)-2 down -3 down , AGV 2: 1 up (16, 18, and 48)-1 down -5 up -5 down . If the max{C 1 , C 2 } of the two AGVs becomes smaller after adjustment, this adjustment should be taken and all entrances will be further checked similarly; otherwise, the original route will be restored.

Initial Population Generation Method. e target population composes of Λ target individuals. e hth target individual is
where τ is the current iteration number, h � 1, 2, . . ., Λ, and Λ is the population size. Two initial individuals that are easy to calculate are given in the initial population.
In the first individual X 0 1 , AGV number is 1. First, sort the pickup points involved in each vertical passageway in increased order and this sequence will make the AGVs drive in one direction when picking up the goods from the same passageway. en connect these sequences by following the order of the pickup passageway from left to right. e sequence can be reversed if reverse make the sequence closer to the previous pickup point. Finally, the individual corresponding to the obtained AGV pickup sequence is taken as the initial individual X 0 1 . In the second individual X 0 2 , AGV quantity is the number of vertical passageways, i.e., each AGV is responsible for only pickup points of one passageway respectively, and takes the goods in increasing order of the pickup point.
e individual corresponding to the obtained AGV pickup sequences is taken as the initial individual X 0 2 . e remaining individuals of the initial population are randomly generated. ey are divided into two halves at first, and then generated from the two initial individuals X 0 1 and X 0 2 , respectively. Repeatedly execute the operation of exchanging randomly two gene locations in the new obtained individual until the initial population is generated.

Mutation, Crossover, and Selection Operations.
In the iterative process, individuals in the parent population are successively mutated, crossed, and selected [42]. e detailed processes for these operations are provided in the following subsections:

Mutation Operation.
e mutation process of NSDE resembles that of the classic DE algorithm. Four parent individuals and one Pareto individual are randomly selected for mutation operation to generate one mutation individual that is expressed as For each target individual in parent population X τ−1 h , the relevant mutation individual is generated as by (17): where X τ−1 Pareto is a Pareto individual chosen at random from the current Pareto frontier individual set, and X τ−1 b 2 are four different parent individuals at iteration (τ-1), a 1 , b 1 , a 2, and b 2 are random integers in [1, Λ], and Z∈ [0, 1] is a mutation coefficient. e operator ⊗ is calculated by equation (18): where G τ h � [g τ h,1 , g τ h,2 , . . . , g τ h,n ] is a temporary individual. e operator ⊕ is calculated by (19): where operator mod is modulus. Individual V τ h maybe not valid because some encoding elements may be repeated or loss. An instance is provided in Tables 2-4 to explain the mutation process.

Crossover Operation. A mutant individual can increase the perturbation of the target individual and jump out of local optimality. A valid trial individual
can be obtained by the crossover operation after mutation operation [43]. In the crossover operation process, select parts of the mutant individual V τ h and add it into the parent individual X τ−1 h to obtain a trial individual U τ h . A three-point insertion method is applied to enhance the perturbation. e insertion points p1, p2, and p3 are three different locations of the target individual, where three parts of the mutant individual are inserted, and they are different random integers in [1, n] Table 4, the crossover operation of X τ−1 h is displayed in Table 5. e crossover coefficient Y � 0.6 and the insertion points are p1 � 2, p2 � 4 and p3 � 8.

Non Dominated Sorting Selection Operation.
e selection process of NSDE introduced the fast nondominated sorting method, crowding comparison, and elite selection strategy in NSGA-II. After mutation and crossover, parent population V τ−1 and progeny population U τ are merged into the set C τ whose size is Λ × 2.
e main process of fast nondominated sorting of individuals in the two generation population set C τ is as follows: (1) First, judge the dominate relationship among all individuals in C τ . For any individual X in C τ , two parameters S X and n X need to be calculated, where S X is a set that stores all individuals dominated by individual X and its initial value is null; n X is a variable that records the quantity of all individuals in C τ that can dominate X and its initialization value is 0. (2) Perform fast nondominated sorting, calculate the nondominated rank (layer) rank X of each individual, and use the set F i to store the result of each rank after sorting. When n X � 0, no individual in C τ can dominate X and all individuals conforming to n X � 0 in C τ are taken out and put into set F 1 that is used to collect all individuals with the highest nondominated rank in C τ . en assign each individual in F 1 the nondominated rank rank X � 1. For the individual set S X dominated by each individual X in F 1 , the corresponding parameter n X' of each individual X′ in S X minus 1. If n X' − 1 � 0, X′ is the individual with the highest nondominated rank in current individual set C τ . en delete individual X′ from C τ and add it into set F 2 . Similarly, continue the above operation until C τ is empty to stop [44]. Tables 6 and 7 is an example of fast nondominated sorting, in which Table 6 represents the objective function value and dominance relation of individuals in the two generations of population, and Table 7 represents the calculation results of nondominated sorting rank.
After the fast nondominated sorting of individuals in C τ , each individual is divided into different levels. In the selecting process of the new target population, it is necessary to inherit as many genes of excellent individuals as possible to the next generation to ensure the good distribution of the population. To judge the environment density of individual, the crowding distance is used to calculate the Euclidean distance between the individual and two adjacent individuals. e greater the crowding distance, the better the individual distribution. e calculation method of crowding distance is as follows: e individuals are sorted in ascending order according to each objective function value of them. e crowding distance of the first and last individuals are all set to infinity, and the crowding distance cd i of individual X i (2 ≤ i ≤ Λ × 2-1) is calculated by the equation below.      e selection operation of NSDE algorithm adopts the elite strategy, i.e., integrates the two generations of population set C τ after fast nondominated sorting and crowding distance calculation, and half of them are reserved as the next parent population. Priority is given to individuals in the higher level, and those with large crowding distance are preferred in the process of selecting individuals in the same level to ensure the diversity of population. For the example in Table 6, four individuals need to be selected. e crowding distance in first layer F 1 are calculated as cd 1 � ∞, cd 2 � 0.45, cd 4 � 0.475, cd 6 � ∞, and cd 8 � 0.575. e four individuals with the largest crowding distance are selected as X 1 , X 4 , X 6 , X 8 .

Algorithm Flowchart.
e specific process of NSDE algorithm is as follows: Step 1: Set Algorithm parameters: population size Λ, mutation coefficient Z, crossover coefficient Y and maximum iteration τ max and set the Pareto solution set to be null.
Step 2: Generate initial population and update Pareto solution set.
Step 3: Set the number of iterations τ � 0 and evolution begins.
Step 4: Mutation and crossover operations are performed to produce progeny individuals.
Step 5: Execute a nondominated sorting on the collection C τ which merges the current population and the offspring population to obtain nondominated layers F i (i � 1, 2, 3, . . .).
Step 6: e individuals in each nondominated layer F i (i � 1, 2, 3, . . .) are put into the parent population in turn to replace the previous ones. For each F i , calculate whether it will exceed the population size before putting it in population: if it exceeds, calculate the number of individuals n' that can be put into the target population, calculate the crowding distance for each individual in F i , select the best n' individuals, and go to Step 7.
e flowchart of NSDE algorithm is expressed in Figure 5.

Computational Results
In order to verify the performance of the proposed NSDE algorithm, this section compares it with Gurobi and multiobjective tabu search (MOTS) algorithm [17], respectively. All algorithms are implemented in C++ programming language and run on a PC with an Intel Corei7 CPU (2.6 GHz × 12) and 8 GB RAM. e relevant parameters of NSDE algorithm are measured by orthogonal experimental design method: mutation coefficient: 0.1, crossover coefficient: 0.3, population size: 50, and maximum number of iterations: 400. e relevant parameters of MOTS algorithm are set by orthogonal experimental design method: tabu list length: 20, maximum number of iterations: 500. e process of the orthogonal experimental design for the NSDE algorithm is presented in Appendix . e performance indicators for multi-objective evolutionary algorithms are generally categorized into four groups based on performance criteria, i.e., counting indicators, convergence indicators, diversity indicators, and comprehensive indicators [45,46], in which each indicator is selected in this paper. Four evaluation indexes used to evaluate the Pareto solution set obtained by multi-objective algorithms are chosen as follows: (1) Overall nondominated vector generation (ONVG): number of vectors in Pareto front [47], and a larger ONVG value indicates more superior performance (this metric is actually applied to count the quantity of approximate Pareto solutions obtained by the correlation experimental algorithms); (2) Convergence metric (CM): an indicator to evaluate convergence. It measures the dominant relationship between the two Pareto fronts E 1 and E 2 . C(E 1 , E 2 ) reflects the percentage of solutions in E 2 dominated by E 1 or equal to E 1 . A larger C(E 1 , E 2 ) value indicates better convergence of Pareto fronts E 1 . e calculation formula is as follows [48]: (1) Set Algorithm parameters: population size Λ, mutation coefficient Z, crossover coefficient Y and maximum iteration τ max , and set the Pareto solution set to be null.
(2) Set τ= 0, generate initial population and update the Pareto solution set.
Perform mutation and crossover operations and produce progeny individuals (3) Spacing metric (SM): an indicator to evaluate the uniformity of the solution set distribution. e smaller the SM value, the more uniform the distribution. e calculation formula is as follows: where d a is the Euclidean distance between solution a and the solution nearest to a, E is the Pareto front, and d � a∈E d a /|E| is the average distance [49]. (4) Hypervolume (HV): the sum of the hypervolume of the hypercube composed of all nondominated solutions and nadir point. e convergence, distribution uniformity and universality of the algorithm results can be evaluated simultaneously [50].

Comparison of NSDE and Gurobi.
In the third-party optimizer evaluation conducted by Decision Tree for Optimization Software, the world's most famous professional optimizer evaluation website, Gurobi shows faster optimization speed and accuracy, which has been proved to have obvious advantages over other large-scale optimizers. e Gurobi optimization solver can be used to solve multi-objective optimization problems, which is one of its greatest strengths. is section compares the performance of NSDE and Gurobi. In order to adapt to the solution of the proposed problem, Gurobi's solving process is designed as follows: since the objective of AGV number cannot be solved directly by Gurobi's function, the number of AGVs is fixed and attempted from small to large, and the proposed problem is transformed into a single-objective optimization problem. en check the optimal solution and the feasibility of the corresponding path. Finally, form the Pareto solution set to output. e layout of the automated warehouse is set as follows: the number of columns of goods (excluding passageways) is fixed at 20, the number of rows of goods r is taken in the set {20, 30}, and the number of picking points n � {10, 20, 30, 40}. ese eight combination experiments are carried out, and five groups of random data are tested for each combination. Because the solution time of Gurobi is too long, after recording the running time of each NSDE experiment, Gurobi runs the same time to stop. e relevant experimental results are shown in Table 8 which records the number of Pareto solutions and running time obtained by the two algorithms under different scales.
When n � 10 to 40, the number of Pareto solutions obtained by NSDE algorithm is significantly more than that of Gurobi, and when n � 30 and 40, Gurobi cannot even obtain a feasible solution to solve the presented problem. It takes some time in Gurobi to preprocess and on judge the conflict of the searched solutions. erefore, Gurobi's solving capability is limited. For the scale of n � 30 and above, Gurobi cannot get the same number of solutions as NSDE or even a feasible solution in the same time, but NSDE still has good performance.

Comparison of NSDE and MOTS.
is section modifies the novel tabu search algorithm proposed by Xing et al. [17], so that it can solve the multi-objective problem, which can be used as a comparison algorithm of NSDE to evaluate the performance. e modified NTS algorithm, i.e., a MOTS algorithm, introducing the judgment of dominant relationship. In the comparative experiment between NSDE and MOTS, the layout of the automated warehouse is set as follows: the number of columns of goods (excluding passageways) is fixed at 20 columns, and the row number of goods r is taken in the set {20, 30}. Eight combination experiments with the number of pickup points n � {60, 120, 140, 180} are carried out, 10 groups of random data are tested for each combination, and the relevant experimental results are recorded. Four evaluation indexes are used for performance comparison, including ONVG, CM, SM, and HV. erefore, in the selection of HV reference point, the point when both objectives take the minimum value is selected and smaller HV value indicate the superiority of the algorithm. Because MOTS is time-consuming, after recording the running time of each NSDE experiment, MOTS runs the same time to stop. e relevant experimental results are shown in Tables 9-12 and Figures 6-9. ese tables, respectively, record ONVG, CM, SM, HV, and running time obtained by the two algorithms under different scales.          better Pareto solutions, avoid too many inferior individuals in the population, and then greatly reduce the calculation time. According to the average level of the four metrics in Figures 6-9, NSDE is superior to MOTS in most cases. In general, NSDE performed well in all aspects, proving the advantages of NSDE's population-based evolution mechanism and elite selection strategy.

Conclusions
is paper studies the biobjective multi-AGV routing problem in an automatic warehouse. In order to make use of the greater advantages of a population-based intelligent optimization algorithm in solving multi-objective optimization problems, a NSDE algorithm is designed. Firstly, a novel encoding and decoding method is designed to eliminate the influence of repeated path, prevent too many bad solutions, and save time.
en the fast nondominated sorting method and elite selection strategy of NSGA-II are introduced to improve the selection operation of differential evolution algorithm, so that better individuals in the offspring population can be retained to continuously update the Pareto solution set. Finally, the effectiveness of the algorithm is proved by numerical experiment.
Some limitations of the study include as follows: When the proposed NSDE algorithm runs at a late stage, the population may be updated slowly sometimes. In this case, future studies will consider the design of a population restart strategy to further improve computing performance. Due to the variability of practical problems, it will be continued to study the following aspects in the future: for the multi-AGV path optimization problem, the limited capacity of AGV when the quality or volume of goods is large, and solution method of the online state will be considered. In these cases, the types of conflict may be more complex.

18
Mathematical Problems in Engineering