# Intelligent Warehouse Robot Scheduling System Using a Modified Nondominated Sorting Algorithm

Authors: Jia Ma; Shujun Yang; Hao Jing
Journal: Discrete Dynamics in Nature and Society
 (2022)

---

## Abstract

In intelligent warehouse, the problem of transporting goods in intelligent warehouse is becoming increasingly complex, and the traditional way of automatically guiding vehicles (AGVs) is inefficient, so automated robot systems are introduced into intelligent warehouses. In this paper, a task assignment model for robots is presented with the transportation problem of robots in intelligent warehouse as the research background. To solve the robot task assignment problem in intelligent warehouse, a novel Pareto-based multiobjective optimization algorithm (MOEA) is proposed, and the aggregation function is invoked to replace the crowding distance; the brain storm operator is used for crossover and mutation. Finally, the ability of the algorithm to solve the benchmark test problem suite and real-world problems is experimentally confirmed.


---

## Body

## 1. Introduction

Recently, the logistics industry is facing more fierce competition, forcing the logistics industry to adopt cheap and efficient automatic robots to replace the traditional automatic guidance vehicles (AGVs) and adopt new intelligent robot scheduling system to reduce operating costs and improve storage efficiency [1]. From the perspective of warehouse management, it is important to quickly transport goods in the warehouse to reduce the production and transaction cycles. In other words, the basic task of intelligent warehouse multirobot scheduling system is to store, transport, and extract goods [2]. Therefore, the task allocation program needs to allocate tasks rationally and effectively. In the process of multirobot task scheduling, a group of tasks need to be scheduled in an optimal way (for example, allocation and execution).At present, there are some metaheuristic algorithms to handle the robot allocation problems. Whale optimization algorithm (WOA) is used to handle mobile robot assignment in the intelligent manufacturing system [3]. Hyperheuristic algorithms based on stigmergy and flocking are used to solve the robot search problem in the unmanned space [4]. A generalized graph-based heuristic algorithm is used to solve a dynamic route planning problem during robot movement [5]. Sam uses metaheuristic algorithm to provide carton manufacturers with an optimal robot transportation strategy for cases with more autonomous robots [6]. Faiza proposes a (Grey wolf optimization and particle swarm optimization) PSO-GWO algorithm to solve the problem of mobile robot avoidance of obstacles [7]. Liu proposes two dynamic order planning algorithms to reduce the order delay problem in smart warehouses [8]. Tan solves the vertical picking problem in the warehouse using the PSO algorithm [9].However, most recent work on task assignment in intelligent warehouse has focused only on single-objective minimization, i.e., it has considered only the minimization of overall time and not the minimization of individual autonomous robot time [10]. For this reason, we build optimization problems based on the characteristics of robot task assignment in intelligent warehouse, which includes two objective functions; the first objective function is minimization of overall time for the robots to perform tasks and the other objective function is to maximize the time a single robot spends performing tasks. Therefore, this is a multiobjective optimization problem (MOP) [11], where multiple objectives need to be optimized simultaneously. Meanwhile, with the increase of orders in warehouse, the scheduling of multiple robots in the intelligent warehouse will become more difficult [12].For real-world MOPs with multiple attributes, the multiobjective optimization algorithm (MOEA) appears [13–15]. Among these MOEAs, Pareto-based method has always been an effective method to solve MOPs [16], and the most representative algorithms are NSGA-II [17] and SPEA2 [18]. In these algorithms, nondominated sorting has great advantages in solving 2,3-objective optimization problems, it allows nondominated stratification of the population and then selects individuals which satisfy the conditions to the next generation.In summary, a novel algorithm which uses nondominated sorting and maximin aggregation function was proposed for the task allocation problem in intelligent warehouse. In this algorithm, nondominated sorting is used to quickly select solutions, and maximin function and one-by-one strategy are used to maintain the uniform distribution of the population. Finally, brain storm operator is used to generate new individuals [19]. Some effective optimization properties of maximin aggregation function will be introduced in the following sections. For the proposed algorithm, the maximin aggregation function can better replace the crowdedness distance to evaluate the contribution of individuals in the same nondominated layer to the new population, the one-by-one comparison strategy can more accurately select a more suitable new solution from the candidate solutions, and the randomness of the brain storm optimization operator helps the algorithm to jump out of the local optimum and enhance the global search ability of the algorithm.The contributions of this paper are as follows:(i)
The task assignment problem for the intelligent warehouse incorporating minimization of overall multirobot time and maximization of individual robot time is proposed. For this task scheduling model, a different crossover mutation operator is adopted.(ii)
For this multirobot scheduling problem, we propose a new Pareto-based algorithm, which also uses maximin function and one-by-one comparison strategy as a supplement, and uses Brain Storm operator to select the parents for crossover and mutation.(iii)
The results of simulation experiments on ZDT and DTLZ benchmark test suites and task assignment problem in intelligent warehouse verify the excellent performance of the algorithm in solving MOPs.The other sections are as follows. Section2 presents a mathematical model for robot task assignment in the intelligent warehouse. In Section 3, related technologies and the proposed algorithm MB-NSGA-II are shown. The experimental results are summarized and analyzed in Section 4. Finally, Section 5 gives the conclusion.


## 2. Problem Formulation

### 2.1. Background Description

For the intelligent warehouse, a batch of goods after delivery to the warehouse, the manager will provide a list of tasks to the intelligent warehouse system based on demand [20]. What the task assignment in the intelligent warehouse needs to do is to reasonably allocate this group of task sequences to the autonomous robots. Among them, the tasks are divided into the following categories: (1) Inbound task, transporting goods to the warehouse; (2) Transportation task is to move goods from one shelf to another according to the needs of the manager; (3) Outbound task, transporting goods out of the warehouse.The two-dimensional plane diagram of the intelligent warehouse is shown in Figure1, in which the neatly arranged look grid represents the shelves. The white cells in the shelves indicate that the shelves are empty because there are no goods stored on them. The automatic robot walks in the aisle to reach the position of the shelf to transport goods. In this way, the information such as the shelf of the intelligent warehouse and the walking range of the robot can be clearly displayed. In the example Figure 1, all goods enter the warehouse from the bottom right and leave the warehouse from the top left.Figure 1
Floor plan of intelligent warehouse.For the robot task assignment issues, the two objective functions of the task assignment model are minimizing the total time for robots to perform tasks and minimizing the time for a single robot to perform tasks, respectively. Assuming that there arem intact robots that can move freely and n tasks to be assigned in the intelligent warehouse, the robots can assign inbound tasks, transportation tasks, and outbound tasks, respectively. When a robot completes a series of tasks, it will encounter the following situations. For the inbound task and transportation task, the coordinates of task ti are assumed to be (xi, yi). When the time consumption between task and task is not considered, the time consumption of robot executing task ti is TCinbound(ti). Therefore, time required for outbound task TCoutbound(ti) is(1)TCoutboundti=xi−xin+yi.The time consumption of the inbound taskTCinbound(ti) is(2)TCinboundti=xi+yi−yout.Unlike the inbound task and the outbound task, the transportation task is different. In the transportation task, assuming that taskti needs to transport goods from (xi, yi) to (mi, ni), its formula is as follows:(3)TCtransti=xi−mi+yi−ni.In the robot task scheduling model, in addition to the time consumption of autonomous robot executing tasks, there is also the time consumption between tasks. Assuming that the destination (xi, yi) of a robot executing task ti and the starting point of the next task zi to be executed by the robot is (mi, ni), the time consumption between two tasks is(4)Cbetweenti,zi=mi−xi+ni−yi.Therefore, the multirobot task assignment problem is transformed into two objective functions: maximizing the time consumption of a single robot and minimizing the total time for robots to perform tasks. Maximizing the time consumption of a single robot refers to the maximum time consumption MRC of a single robot among the robots completing the task, as shown in the following formula:(5)MRC=maxTCtotalri,TCtotalri,Si=∑j=1kTCti,j+∑n=2kCbetweenti,n−1,ti,n,where TCtotal denotes the total time consumed by the robot ri to execute a task sequence Si: ti1⟶ti2⟶ti3⋯⟶tik, tik denotes the ith robot performing the kth task. Cbetween(ti, n-1, ti,n) represents the time consumed by the ith robot when executing the task, k denotes the total number of tasks performed by the ith robot, and n denotes the nth task in execution.Total robot time consumption minimization (MTC) is the second objective function; we need to minimize the total time for robots to perform tasks, and the formula is shown below:(6)MTC=∑i=1NTCtotalri,Si,where N denotes the total number of robots in the intelligent warehouse.All the above time consumption functions can be calculated using the Manhattan distance, since the speed of the robot in the warehouse is set to 1 m/s. The total distance traveled by the robots is also equal to the total time consumed by the robots.


### 2.2. Generation of Offspring

In the iterative optimization of discrete decision variables, some special crossover and mutation methods are used to form new offspring.To solve the robot task assignment problem, we redesigned chromosome part and divided the chromosome into two parts. The first part is a set of task sequences, and the other part contains the number of tasks undertaken by each robot, and the sum of numbers is the total number of tasks in the task list. Figure2 illustrates the relationship between the first part of the chromosome and the second part of the chromosome. As shown in Figure 2, robot 1 performs three tasks, robot 2 performs two tasks, 4 and 5, and robot 3 performs the remaining tasks.Figure 2
Example diagram of task assignment.For discrete decision variables, the conventional crossover and mutation operators [21] are not appropriate. Therefore, for this particular chromosome, two crossover operators are used. For these two parts of the chromosomes, we use order crossover [22] and simulated binary crossover.Mutation operation can usually increase the random exploration ability of the algorithm, preventing algorithms from falling into local optimal. For discrete decision variables, we use slight mutation to replace polynomial mutation. The mutation process is shown in Figure3. Firstly, we select a point from these points as the mutation point, then select a substring from the chromosome, and finally insert the substring behind the mutation point to form a new string of chromosomes.Figure 3
Example diagram of slight mutation.


## 2.1. Background Description

For the intelligent warehouse, a batch of goods after delivery to the warehouse, the manager will provide a list of tasks to the intelligent warehouse system based on demand [20]. What the task assignment in the intelligent warehouse needs to do is to reasonably allocate this group of task sequences to the autonomous robots. Among them, the tasks are divided into the following categories: (1) Inbound task, transporting goods to the warehouse; (2) Transportation task is to move goods from one shelf to another according to the needs of the manager; (3) Outbound task, transporting goods out of the warehouse.The two-dimensional plane diagram of the intelligent warehouse is shown in Figure1, in which the neatly arranged look grid represents the shelves. The white cells in the shelves indicate that the shelves are empty because there are no goods stored on them. The automatic robot walks in the aisle to reach the position of the shelf to transport goods. In this way, the information such as the shelf of the intelligent warehouse and the walking range of the robot can be clearly displayed. In the example Figure 1, all goods enter the warehouse from the bottom right and leave the warehouse from the top left.Figure 1
Floor plan of intelligent warehouse.For the robot task assignment issues, the two objective functions of the task assignment model are minimizing the total time for robots to perform tasks and minimizing the time for a single robot to perform tasks, respectively. Assuming that there arem intact robots that can move freely and n tasks to be assigned in the intelligent warehouse, the robots can assign inbound tasks, transportation tasks, and outbound tasks, respectively. When a robot completes a series of tasks, it will encounter the following situations. For the inbound task and transportation task, the coordinates of task ti are assumed to be (xi, yi). When the time consumption between task and task is not considered, the time consumption of robot executing task ti is TCinbound(ti). Therefore, time required for outbound task TCoutbound(ti) is(1)TCoutboundti=xi−xin+yi.The time consumption of the inbound taskTCinbound(ti) is(2)TCinboundti=xi+yi−yout.Unlike the inbound task and the outbound task, the transportation task is different. In the transportation task, assuming that taskti needs to transport goods from (xi, yi) to (mi, ni), its formula is as follows:(3)TCtransti=xi−mi+yi−ni.In the robot task scheduling model, in addition to the time consumption of autonomous robot executing tasks, there is also the time consumption between tasks. Assuming that the destination (xi, yi) of a robot executing task ti and the starting point of the next task zi to be executed by the robot is (mi, ni), the time consumption between two tasks is(4)Cbetweenti,zi=mi−xi+ni−yi.Therefore, the multirobot task assignment problem is transformed into two objective functions: maximizing the time consumption of a single robot and minimizing the total time for robots to perform tasks. Maximizing the time consumption of a single robot refers to the maximum time consumption MRC of a single robot among the robots completing the task, as shown in the following formula:(5)MRC=maxTCtotalri,TCtotalri,Si=∑j=1kTCti,j+∑n=2kCbetweenti,n−1,ti,n,where TCtotal denotes the total time consumed by the robot ri to execute a task sequence Si: ti1⟶ti2⟶ti3⋯⟶tik, tik denotes the ith robot performing the kth task. Cbetween(ti, n-1, ti,n) represents the time consumed by the ith robot when executing the task, k denotes the total number of tasks performed by the ith robot, and n denotes the nth task in execution.Total robot time consumption minimization (MTC) is the second objective function; we need to minimize the total time for robots to perform tasks, and the formula is shown below:(6)MTC=∑i=1NTCtotalri,Si,where N denotes the total number of robots in the intelligent warehouse.All the above time consumption functions can be calculated using the Manhattan distance, since the speed of the robot in the warehouse is set to 1 m/s. The total distance traveled by the robots is also equal to the total time consumed by the robots.


## 2.2. Generation of Offspring

In the iterative optimization of discrete decision variables, some special crossover and mutation methods are used to form new offspring.To solve the robot task assignment problem, we redesigned chromosome part and divided the chromosome into two parts. The first part is a set of task sequences, and the other part contains the number of tasks undertaken by each robot, and the sum of numbers is the total number of tasks in the task list. Figure2 illustrates the relationship between the first part of the chromosome and the second part of the chromosome. As shown in Figure 2, robot 1 performs three tasks, robot 2 performs two tasks, 4 and 5, and robot 3 performs the remaining tasks.Figure 2
Example diagram of task assignment.For discrete decision variables, the conventional crossover and mutation operators [21] are not appropriate. Therefore, for this particular chromosome, two crossover operators are used. For these two parts of the chromosomes, we use order crossover [22] and simulated binary crossover.Mutation operation can usually increase the random exploration ability of the algorithm, preventing algorithms from falling into local optimal. For discrete decision variables, we use slight mutation to replace polynomial mutation. The mutation process is shown in Figure3. Firstly, we select a point from these points as the mutation point, then select a substring from the chromosome, and finally insert the substring behind the mutation point to form a new string of chromosomes.Figure 3
Example diagram of slight mutation.


## 3. The Proposed Method

In recent research, most of the literature studies tend to solve MOPs based on Pareto-optimal, allowing the algorithm to keep approaching the true Pareto optimal front through nondominated sorting and many improved Pareto-based methods [23]. In addition, some population diversity conservation mechanisms ensure that populations remain well distributed in the objective space after nondominated sorting.


### 3.1. Construction of Nondominated Solution Set

To reduce the high time complexity of constructing nondominated solution sets in NSGA, NSGAII proposed a new nondominated sorting method to select solutions for the new population. Given the good performance of NSGAII in dealing with MOPs, here we use this hierarchical nondominated solution set construction scheme to construct new populations.In the nondominated sorting process, the nondominated individuals are first selected into the first stratum, then, the second stratum is the set of nondominant individuals obtained after removing the first stratum individuals from the population, and so on. In the final selection, individuals in the first stratum are considered first, and then individuals in the second stratum is considered until the new population size is satisfied. Specific details can be found in reference [17]. Algorithm 1 is the process of constructing the nondominated solution set.


### 3.2. Maximin Fitness Function

In the continuous development of multiobjective optimization, some fitness functions are used as indicators to evaluate individuals in the population, among which the maximin function is applied to multiobjective optimization with its own characteristics. The formula is as follows:(7)fitnessi=maxj≠iminkfkxi−fkxj,where k denotes objective from 1 to m, i represents the ith individual, and j denotes any individual except ith individual. And the properties of the maximin function are as follows [24, 25]:(1)
The maximin fitness can reflect the dominance between individuals. The maximin value greater than, equal to, or less than zero means that the individual is a dominated individual, the individual is weakly dominated, and individual is nondominated in the population, respectively.(2)
The maximin fitness reflects the presence or absence of clustering around an individual, as shown in Figure4.(3)
The maximin function can be used to reflect the relationship between the dominated individual and the nondominated front. (see Figure4).Figure 4
Properties of the maximin function.Algorithm 1: Nondominated sorting.
Input: P (population)Output: P1,P2,…,Pn (stratification results)(1)
∀p∈P,Sp=0,Dp=∅,i = 1;// Sp denotes the size of set of solutions that dominate p, Dpi S denotes the set of solutions dominated by p(2)
for∀p∈P(3)
for∀q∈P(4)
ifp≻q then Dp=Dp∪q(5)
elseifq≻p then Sp=Sp+1(6)
end for q(7)
ifSp=0 then P1=P1∪p(8)
end for p(9)
whilePi≠∅ // Pi indicates the number of non-dominated layers(10)
{POP=∅ ;(11)
for∀p∈Pi(12)
for∀q∈sp, nq=nq−1 ;(13)
ifnq=0POP=POP∪q(14)
end for p(15)
i = i + 1;(16)
Pi=POP ;(17)
end for while(18)
end(19)
returnP1,P2,…,Pn


### 3.3. One-by-One Comparison Strategy

The contribution of individual to the population can be effectively assessed by the maximin aggregation function. However, if only the maximin aggregation function is used to select suitable individuals, then the selection process may encounter the equivalence selection dilemma, that is, some individuals have equal maximin fitness, and there is no way to select a better individual from them, and then it is necessary to use a one-by-one comparison strategy to select a better individual from these individuals. The following figure shows an example of using the one-by-one comparison strategy to further select individuals [26].Through one-by-one comparison strategy, the problem of equivalence selection can be solved efficiently. Suppose that three individuals are selected from the four individuals in Figure5. First, the maximin aggregation function is used to evaluate these individuals. A with the minimum value is selected from these individuals to become the first candidate individual, then continue to select individuals from the remaining individuals through the maximin aggregation function, and the individual D is selected into the new population. Next, individuals B and C are compared with the new population by maximin function, and we select individual B to join the new population. Finally, the new population is formed by individuals A, B, D. As can be seen in the figure, the selected individuals remain well distributed.Figure 5
Diagram of the process of solving the equivalence selection problem by one-by-one comparison strategy.
(a)(b)(c)(d)


### 3.4. Brain Storm Operator

Unlike most multiobjective algorithms, we use brain storm operator as a strategy for selecting parents, which increases the ability of the algorithm in terms of exploration and makes it less prone to fall into local optima. We demonstrate this part of the process with Algorithm 2.Algorithm 2: Brain storm operator (P, S).
(1)
Objective normalization (P)(2)
fori = 1: S − 1(3)
Clustering: Dividing the population into m clusters using the k-means method.(4)
pgen=rand0,1(5)
ifpgen < pfixed(6)
pone=rand0,1 //Select a random cluster from m clusters(7)
ifpone < pone-fixed//pone-fixed is a predetermined value(8)
Select the clustering center and the individualPopi as parent generation of population(9)
else(10)
Select a random individual and the individualPopi as parent generation of population(11)
end(12)
end(13)
Compare the newly generated individuals with the present ones by maximin function, where the better one is preserved.(14)
end


### 3.5. Proposed Algorithm: MB-NSGA-II

The flow and framework of the algorithm are described below. In the first step, the populationP is randomly initialized. And the brain storm operator is used to cluster parent population and then select appropriate chromosomes as parent chromosomes based on suitable conditions and generate the new population Q by crossover and mutation. After this, parent P and offspring Q are combined to form P'. Finally, the new population is regenerated by comprehensive selection. Keep cycling through the above steps until the MAX_ Fitness Evaluations is reached. The general framework of MB-NSGA-II is shown in Algorithm 4.


#### 3.5.1. Normalization

In multiobjective optimization algorithms, normalization can effectively solve the problem that different objectives have widely varying ranges of values [27]. The formula is as follows:(8)fm’x=fmx−zmlowerzmupper−zmlower,where zmlower denotes the lower bound of the mth objective function, zmupper denotes the upper bound of the mth objective function, m means this is the mth objective function, and m∈1,2,…,M.Algorithm 3: Comprehensive Selection (Off, S).
Input: Off (Offspring population), S (population size)Output: Pop (new population)(1)
P1,P2,…,Pn  = Nondominated sorting (Off)(2)
Pop=∅ ;(3)
ifn = 1(4)
Pop  = One-by-one comparison strategyOff,z∗,S;(5)
else(6)
whilem(7)
Pop=Pop∪Pi(8)
i = i+1(9)
end for while(10)
Ps  = One-by-one comparison strategySi,Pop,S−Pop ;(11)
Pop=P″∪Ps(12)
end for if(13)
returnPopAlgorithm 4: General Framework of MB-NSGA-II.
(1)
Initialization: Initial populationP containing N randomly individuals(2)
whileFitness Evaluations < MAX_Fitness Evaluationsdo(3)
Q = Brain Storm OperatorP,s ;(4)
Off = P ∪ Q(5)
Pop = Comprehensive _Selection (Off, S);(6)
end for while(7)
returnPopAlgorithm 5: One-by-one comparison strategy S,P″,size.
Input: S: individuals in stratum I, P″: new population, size: number of individuals needed for new populationOutput: Ps: the set of individuals provided to the new population(1)
Pn=∅;(2)
Objective normalization (Q);(3)
if |S|>size(4)
fori = 1: size-1(5)
MF = argmini=1,2,…,Sfitnessi;//individual with minimal maximin fitness between individuals in Si and the new population P(6)
Pn  =  Pn, ≤MF;(7)
end for(8)
end if(9)
ReturnPs


#### 3.5.2. Comprehensive Selection

As shown in Algorithm 3, comprehensive selection consists of the following steps. First, nondominated sorting for the populationPop is executed to divide Pop into strata from 1 to n. If there is only one layer of nondominated layer, the maximin fitness value of the population and the ideal point is calculated, and individuals are selected into the new population one by one. Then, individuals were selected layer by layer starting from the first layer until the new population size is satisfied. When proceeding to stratum i, if the predetermined S is less than the number of individuals in the ith stratum, the set of individuals Ps is selected using a one-by-one comparison strategy (Line 8 Algorithm 5), and finally Ps is merged with Pop to obtain the new population Pop.


#### 3.5.3. One-by-One Comparison Strategy

The one-by-one comparison strategy is shown in Algorithm 5. One-by-one comparison strategy is used to compare individuals at the same stratum and select a subset of individuals from them. Firstly, letPn be empty. First, each objective of the population P″ is normalized by (7).Then, if the population size of S is larger than size, comparing individual fitness values in S using the maximin aggregation function, select the individual with the smallest maximin fitness from it to join the MF, and so on, until Pn reaches the new population size.


### 3.6. Computational Complexity Analysis of MB-NSGA-II

To study the computational cost consumption of the algorithm, we summarize and discuss the time complexity of one generation in Algorithm 3. In addition to crossover and mutation operators, computational complexity is affected by the nondominated sorting and one-by-one comparison strategy.First, normalization requires a computational complexity of O(mN). And the nondominated sorting’s computational complexity in Algorithm 1 is O(mN2). The complexity of calculating the maximin fitness value is O(mN). The step of one-by-one comparison strategy needs O(mN2). In a word, the overall complexity of the algorithm is O(mN2) in one generation.


## 3.1. Construction of Nondominated Solution Set

To reduce the high time complexity of constructing nondominated solution sets in NSGA, NSGAII proposed a new nondominated sorting method to select solutions for the new population. Given the good performance of NSGAII in dealing with MOPs, here we use this hierarchical nondominated solution set construction scheme to construct new populations.In the nondominated sorting process, the nondominated individuals are first selected into the first stratum, then, the second stratum is the set of nondominant individuals obtained after removing the first stratum individuals from the population, and so on. In the final selection, individuals in the first stratum are considered first, and then individuals in the second stratum is considered until the new population size is satisfied. Specific details can be found in reference [17]. Algorithm 1 is the process of constructing the nondominated solution set.


## 3.2. Maximin Fitness Function

In the continuous development of multiobjective optimization, some fitness functions are used as indicators to evaluate individuals in the population, among which the maximin function is applied to multiobjective optimization with its own characteristics. The formula is as follows:(7)fitnessi=maxj≠iminkfkxi−fkxj,where k denotes objective from 1 to m, i represents the ith individual, and j denotes any individual except ith individual. And the properties of the maximin function are as follows [24, 25]:(1)
The maximin fitness can reflect the dominance between individuals. The maximin value greater than, equal to, or less than zero means that the individual is a dominated individual, the individual is weakly dominated, and individual is nondominated in the population, respectively.(2)
The maximin fitness reflects the presence or absence of clustering around an individual, as shown in Figure4.(3)
The maximin function can be used to reflect the relationship between the dominated individual and the nondominated front. (see Figure4).Figure 4
Properties of the maximin function.Algorithm 1: Nondominated sorting.
Input: P (population)Output: P1,P2,…,Pn (stratification results)(1)
∀p∈P,Sp=0,Dp=∅,i = 1;// Sp denotes the size of set of solutions that dominate p, Dpi S denotes the set of solutions dominated by p(2)
for∀p∈P(3)
for∀q∈P(4)
ifp≻q then Dp=Dp∪q(5)
elseifq≻p then Sp=Sp+1(6)
end for q(7)
ifSp=0 then P1=P1∪p(8)
end for p(9)
whilePi≠∅ // Pi indicates the number of non-dominated layers(10)
{POP=∅ ;(11)
for∀p∈Pi(12)
for∀q∈sp, nq=nq−1 ;(13)
ifnq=0POP=POP∪q(14)
end for p(15)
i = i + 1;(16)
Pi=POP ;(17)
end for while(18)
end(19)
returnP1,P2,…,Pn


## 3.3. One-by-One Comparison Strategy

The contribution of individual to the population can be effectively assessed by the maximin aggregation function. However, if only the maximin aggregation function is used to select suitable individuals, then the selection process may encounter the equivalence selection dilemma, that is, some individuals have equal maximin fitness, and there is no way to select a better individual from them, and then it is necessary to use a one-by-one comparison strategy to select a better individual from these individuals. The following figure shows an example of using the one-by-one comparison strategy to further select individuals [26].Through one-by-one comparison strategy, the problem of equivalence selection can be solved efficiently. Suppose that three individuals are selected from the four individuals in Figure5. First, the maximin aggregation function is used to evaluate these individuals. A with the minimum value is selected from these individuals to become the first candidate individual, then continue to select individuals from the remaining individuals through the maximin aggregation function, and the individual D is selected into the new population. Next, individuals B and C are compared with the new population by maximin function, and we select individual B to join the new population. Finally, the new population is formed by individuals A, B, D. As can be seen in the figure, the selected individuals remain well distributed.Figure 5
Diagram of the process of solving the equivalence selection problem by one-by-one comparison strategy.
(a)(b)(c)(d)


## 3.4. Brain Storm Operator

Unlike most multiobjective algorithms, we use brain storm operator as a strategy for selecting parents, which increases the ability of the algorithm in terms of exploration and makes it less prone to fall into local optima. We demonstrate this part of the process with Algorithm 2.Algorithm 2: Brain storm operator (P, S).
(1)
Objective normalization (P)(2)
fori = 1: S − 1(3)
Clustering: Dividing the population into m clusters using the k-means method.(4)
pgen=rand0,1(5)
ifpgen < pfixed(6)
pone=rand0,1 //Select a random cluster from m clusters(7)
ifpone < pone-fixed//pone-fixed is a predetermined value(8)
Select the clustering center and the individualPopi as parent generation of population(9)
else(10)
Select a random individual and the individualPopi as parent generation of population(11)
end(12)
end(13)
Compare the newly generated individuals with the present ones by maximin function, where the better one is preserved.(14)
end


## 3.5. Proposed Algorithm: MB-NSGA-II

The flow and framework of the algorithm are described below. In the first step, the populationP is randomly initialized. And the brain storm operator is used to cluster parent population and then select appropriate chromosomes as parent chromosomes based on suitable conditions and generate the new population Q by crossover and mutation. After this, parent P and offspring Q are combined to form P'. Finally, the new population is regenerated by comprehensive selection. Keep cycling through the above steps until the MAX_ Fitness Evaluations is reached. The general framework of MB-NSGA-II is shown in Algorithm 4.


### 3.5.1. Normalization

In multiobjective optimization algorithms, normalization can effectively solve the problem that different objectives have widely varying ranges of values [27]. The formula is as follows:(8)fm’x=fmx−zmlowerzmupper−zmlower,where zmlower denotes the lower bound of the mth objective function, zmupper denotes the upper bound of the mth objective function, m means this is the mth objective function, and m∈1,2,…,M.Algorithm 3: Comprehensive Selection (Off, S).
Input: Off (Offspring population), S (population size)Output: Pop (new population)(1)
P1,P2,…,Pn  = Nondominated sorting (Off)(2)
Pop=∅ ;(3)
ifn = 1(4)
Pop  = One-by-one comparison strategyOff,z∗,S;(5)
else(6)
whilem(7)
Pop=Pop∪Pi(8)
i = i+1(9)
end for while(10)
Ps  = One-by-one comparison strategySi,Pop,S−Pop ;(11)
Pop=P″∪Ps(12)
end for if(13)
returnPopAlgorithm 4: General Framework of MB-NSGA-II.
(1)
Initialization: Initial populationP containing N randomly individuals(2)
whileFitness Evaluations < MAX_Fitness Evaluationsdo(3)
Q = Brain Storm OperatorP,s ;(4)
Off = P ∪ Q(5)
Pop = Comprehensive _Selection (Off, S);(6)
end for while(7)
returnPopAlgorithm 5: One-by-one comparison strategy S,P″,size.
Input: S: individuals in stratum I, P″: new population, size: number of individuals needed for new populationOutput: Ps: the set of individuals provided to the new population(1)
Pn=∅;(2)
Objective normalization (Q);(3)
if |S|>size(4)
fori = 1: size-1(5)
MF = argmini=1,2,…,Sfitnessi;//individual with minimal maximin fitness between individuals in Si and the new population P(6)
Pn  =  Pn, ≤MF;(7)
end for(8)
end if(9)
ReturnPs


### 3.5.2. Comprehensive Selection

As shown in Algorithm 3, comprehensive selection consists of the following steps. First, nondominated sorting for the populationPop is executed to divide Pop into strata from 1 to n. If there is only one layer of nondominated layer, the maximin fitness value of the population and the ideal point is calculated, and individuals are selected into the new population one by one. Then, individuals were selected layer by layer starting from the first layer until the new population size is satisfied. When proceeding to stratum i, if the predetermined S is less than the number of individuals in the ith stratum, the set of individuals Ps is selected using a one-by-one comparison strategy (Line 8 Algorithm 5), and finally Ps is merged with Pop to obtain the new population Pop.


### 3.5.3. One-by-One Comparison Strategy

The one-by-one comparison strategy is shown in Algorithm 5. One-by-one comparison strategy is used to compare individuals at the same stratum and select a subset of individuals from them. Firstly, letPn be empty. First, each objective of the population P″ is normalized by (7).Then, if the population size of S is larger than size, comparing individual fitness values in S using the maximin aggregation function, select the individual with the smallest maximin fitness from it to join the MF, and so on, until Pn reaches the new population size.


## 3.5.1. Normalization

In multiobjective optimization algorithms, normalization can effectively solve the problem that different objectives have widely varying ranges of values [27]. The formula is as follows:(8)fm’x=fmx−zmlowerzmupper−zmlower,where zmlower denotes the lower bound of the mth objective function, zmupper denotes the upper bound of the mth objective function, m means this is the mth objective function, and m∈1,2,…,M.Algorithm 3: Comprehensive Selection (Off, S).
Input: Off (Offspring population), S (population size)Output: Pop (new population)(1)
P1,P2,…,Pn  = Nondominated sorting (Off)(2)
Pop=∅ ;(3)
ifn = 1(4)
Pop  = One-by-one comparison strategyOff,z∗,S;(5)
else(6)
whilem(7)
Pop=Pop∪Pi(8)
i = i+1(9)
end for while(10)
Ps  = One-by-one comparison strategySi,Pop,S−Pop ;(11)
Pop=P″∪Ps(12)
end for if(13)
returnPopAlgorithm 4: General Framework of MB-NSGA-II.
(1)
Initialization: Initial populationP containing N randomly individuals(2)
whileFitness Evaluations < MAX_Fitness Evaluationsdo(3)
Q = Brain Storm OperatorP,s ;(4)
Off = P ∪ Q(5)
Pop = Comprehensive _Selection (Off, S);(6)
end for while(7)
returnPopAlgorithm 5: One-by-one comparison strategy S,P″,size.
Input: S: individuals in stratum I, P″: new population, size: number of individuals needed for new populationOutput: Ps: the set of individuals provided to the new population(1)
Pn=∅;(2)
Objective normalization (Q);(3)
if |S|>size(4)
fori = 1: size-1(5)
MF = argmini=1,2,…,Sfitnessi;//individual with minimal maximin fitness between individuals in Si and the new population P(6)
Pn  =  Pn, ≤MF;(7)
end for(8)
end if(9)
ReturnPs


## 3.5.2. Comprehensive Selection

As shown in Algorithm 3, comprehensive selection consists of the following steps. First, nondominated sorting for the populationPop is executed to divide Pop into strata from 1 to n. If there is only one layer of nondominated layer, the maximin fitness value of the population and the ideal point is calculated, and individuals are selected into the new population one by one. Then, individuals were selected layer by layer starting from the first layer until the new population size is satisfied. When proceeding to stratum i, if the predetermined S is less than the number of individuals in the ith stratum, the set of individuals Ps is selected using a one-by-one comparison strategy (Line 8 Algorithm 5), and finally Ps is merged with Pop to obtain the new population Pop.


## 3.5.3. One-by-One Comparison Strategy

The one-by-one comparison strategy is shown in Algorithm 5. One-by-one comparison strategy is used to compare individuals at the same stratum and select a subset of individuals from them. Firstly, letPn be empty. First, each objective of the population P″ is normalized by (7).Then, if the population size of S is larger than size, comparing individual fitness values in S using the maximin aggregation function, select the individual with the smallest maximin fitness from it to join the MF, and so on, until Pn reaches the new population size.


## 3.6. Computational Complexity Analysis of MB-NSGA-II

To study the computational cost consumption of the algorithm, we summarize and discuss the time complexity of one generation in Algorithm 3. In addition to crossover and mutation operators, computational complexity is affected by the nondominated sorting and one-by-one comparison strategy.First, normalization requires a computational complexity of O(mN). And the nondominated sorting’s computational complexity in Algorithm 1 is O(mN2). The complexity of calculating the maximin fitness value is O(mN). The step of one-by-one comparison strategy needs O(mN2). In a word, the overall complexity of the algorithm is O(mN2) in one generation.


## 4. Experimental Study

The experimental part contains two sections; the first section discusses and analyzes the effectiveness of MB-NSGA-II on the benchmark test suites; the second part describes how we use our proposed algorithm to solve practical problems in intelligent warehouse.


### 4.1. Experimental Results on Benchmarks

This section provides the experimental results of MB-NSGAII with two other start-of-art multiobjective algorithms, NSGAII and IBEA, on the benchmark test suites DTLZ [28] and ZDT [29]. On DTLZ, each test instance was experimented on the 2,3,5- objectives, and on ZDT, the experiments were run on the 2-objective. Each algorithm needs to be run 20 times on each test issue. We can discuss the results of this experiment in terms of comprehensive performance indicator IGD.


#### 4.1.1. Results and Discussion on ZDTs and DTLZs

The ZDT test suite contains six test problems with different characteristics, each test instance involving a feature that would cause the algorithm to have difficulty converging during optimization (e.g., multimodality). Through these test instances, it is possible to analyze which problems the algorithm excels in. DTLZ designs the problems through a systematic approach, and its decision variables and objectives are scalable, facilitating the algorithm to experiment in solving MaOPs. DTLZ increases the difficulty of the test problem by introducing manageable difficulties, facilitating the algorithm to experiment on test problems with any number of objectives. To more fully validate the effectiveness of several compared algorithms in solving MOPs, we chose to conduct experiments on 2,3, and 5 objectives.M in Table 1 indicates the objective number. The numbers of decision variable are set, as shown in Table 1. Details can be found in reference [28,29].Table 1
Setting of decision variables.
ProblemNumber of decision variablesZDT1-330ZDT4,610ZDT580DTLZ1M + 4DTLZ2-4M + 9To more visually reflect the ability of the algorithm on the benchmark test instances, we choose the comprehensive evaluation indicator HV and IGD as the basis for judging the performance of the algorithms. These two metrics are composed of different formulas, so the final population obtained by the algorithm can be evaluated from different perspectives. The formula and detailed description can be found in [30,31].


#### 4.1.2. Experimental Settings

(i)
Population size and termination conditions: we use population of size 100 for DTLZ with 2,3 objectives and 200 for DTLZ with 5 objectives. For ZDTs, we select population of size 100 for the three algorithms. For DTLZ1-4 and ZDT1-6, the algorithms will end after 500 generations.(ii)
Crossover and mutation: in this experiment, binary crossover and polynomial mutation were used as operators for crossover and mutation. Set 1.0 as the crossover probability and set 1/D as the probability of mutation, where D is the number of decision variables.(iii)
Selection of indicator: for IBEA, we selectIε+ as a performance indicator to evaluate the strengths and weaknesses of individuals. The details of the indicator can be seen in [32].


#### 4.1.3. Analysis of Experimental Results

Comparison data of MB-NSGA-II, NSGA-II, and IBEA on DTLZ and ZDT are presented in Tables2–5. These tables show the IGD and HV means and standard deviations of these compared algorithms on the benchmark test problems, and highlights the algorithm that worked best on these test problems. The Wilcoxon rank sum test with a significance level of 0.05 was used for the standard deviation analysis when conducting the experiment. The symbols “+,” “−,” and “=” are used to indicate that the algorithms are significantly better, significantly worse, or not significantly different between the compared algorithms and MB-NSGA-II.Table 2
Experimental results of IGD on ZDTs (mean and standard deviation).
ProblemMNSGAIIIBEAMBNSGAIIZDT124.7879e-3 (1.31e-4) -4.0901e-3 (1.01e-4) -3.9627e-3 (6.92e-5)ZDT224.9504e-3 (2.24e-4) +8.4212e-3 (7.74e-4) -5.3033e-3 (1.32e-4)ZDT325.4102e-3 (1.31e-4) -1.6020e-2 (8.05e-4) -5.0404e-3 (1.24e-4)ZDT425.0667e-3 (4.10e-4) +2.2761e-2 (1.59e-2) -6.3474e-3 (2.30e-4)ZDT525.1086e-1 (8.31e-2) +1.8513e+0 (4.39e-1)=1.6593e+0 (7.13e-1)ZDT623.7200e-3 (1.20e-4)=4.4425e-3 (1.14e-4) -3.7006e-3 (5.32e-5)Table 3
Experimental results of IGD on several DTLZs (mean and standard deviation).
ProblemMNSGAIIIBEAMBNSGAIIDTLZ122.2068e-3 (7.35e-5) -7.9965e-2 (7.99e-3) -2.0230e-3 (5.45e-5)32.8447e-2 (1.33e-3) -1.6616e-1 (2.51e-2) -2.2376e-2 (3.68e-4)51.7146e-1 (1.25e-1) -1.8372e-1 (2.16e-2) -5.1833e-2 (3.37e-4)DTLZ225.0483e-3 (1.75e-4) +1.6671e-2 (1.60e-3) -8.5238e-3 (2.92e-4)37.3058e-2 (2.57e-3) =8.4162e-2 (2.51e-3) -7.3973e-2 (1.57e-3)52.0563e-1 (4.67e-3) -1.9251e-1 (1.54e-3) -1.7906e-1 (2.44e-3)DTLZ327.1268e-3 (1.39e-3) +3.3913e-1 (8.73e-3) -9.8870e-3 (1.10e-3)31.4200e-1 (2.06e-1) =4.7891e-1 (4.51e-3) -7.8834e-2 (2.64e-3)57.5535e-1 (7.31e-1) -5.9365e-1 (9.76e-3) -1.8301e-1 (2.58e-3)DTLZ421.5241e-1 (3.11e-1) -4.5174e-1 (3.75e-1) -8.1849e-2 (2.32e-1)31.5811e-1 (2.77e-1) +8.2918e-2 (2.68e-3) =2.6157e-1 (2.41e-1)52.0630e-1 (3.41e-3) -2.1267e-1 (7.27e-2) -1.8232e-1 (1.93e-3)Table 4
Experimental results of HV on several DTLZs (mean and standard deviation).
ProblemMNSGAIIIBEAMBNSGAIIDTLZ125.8121e-1 (3.60e-4) -3.9969e-1 (1.87e-2) -5.8162e-1 (3.86e-4)38.2082e-1 (5.02e-3) -4.7837e-1 (5.50e-2) -8.3698e-1 (1.86e-3)56.9302e-1 (3.39e-1) -7.5172e-1 (4.51e-2) -9.7599e-1 (1.40e-3)DTLZ223.4654e-1 (1.97e-4) -3.4613e-1 (2.17e-4) -3.4731e-1 (1.09e-4)35.2551e-1 (3.99e-3) -5.5477e-1 (1.48e-3) -5.5689e-1 (9.61e-4)56.7753e-1 (7.73e-3) -8.0904e-1 (1.36e-3) +8.0024e-1 (2.49e-3)DTLZ323.4146e-1 (2.73e-3) =1.6879e-1 (3.86e-3) -3.4251e-1 (2.27e-3)34.5747e-1 (1.56e-1) -2.3645e-1 (9.09e-3) -5.4725e-1 (3.20e-3)54.0235e-1 (3.46e-1) -3.8077e-1 (6.39e-3) -8.0179e-1 (4.10e-3)DTLZ422.9554e-1 (1.08e-1) -1.9305e-1 (1.32e-1) -3.2164e-1 (8.11e-2)34.8619e-1 (1.39e-1) =5.5574e-1 (9.75e-4) =4.7330e-1 (1.09e-1)56.8256e-1 (8.53e-3) -8.0150e-1 (2.60e-2) =8.0780e-1 (1.73e-3)Table 5
Experimental results of HV on ZDTs (mean and standard deviation).
ProblemMNSGAIIIBEAMBNSGAIIZDT127.1925e-1 (1.76e-4) -7.2017e-1 (1.36e-4) -7.2031e-1 (8.51e-5)ZDT224.4399e-1 (2.03e-4) -4.4410e-1 (1.58e-4) -4.4485e-1 (6.70e-5)ZDT325.9938e-1 (8.34e-5) -5.9816e-1 (1.17e-4) -5.9973e-1 (4.53e-5)ZDT427.1800e-1 (1.09e-3) =7.0312e-1 (9.85e-3) -7.1735e-1 (5.54e-4)ZDT528.1676e-1 (1.20e-2) =8.1119e-1 (9.16e-3)8.1451e-1 (3.26e-4)ZDT623.8826e-1 (1.18e-4) =3.8766e-1 (1.17e-4) -3.8832e-1 (5.60e-5)In Table2 and Table 3, MB-NSGA-II has an outstanding performance on both DTLZ test problems and ZDT test problems, and it ranks first on 11 out of 18 test instances, and it is also competitive on other test instances. In particular, it achieved the best results on DTLZ1 for both 2,3,5 objectives, and also ranked first on DTLZ3 for two objective problems. Of course, NSGA-II also had good results, obtaining the best results on the six problems. In Table 3, MB-NSGAII ranked first in performance on 8 test problems, while NSGAII values achieved good results on three problems. It can be seen that MB-NSGA-II achieved the best results on all test problems on DTLZ1 with hyperplane PF, which is a significant improvement over NSGAII. On DTLZ2, MB-NSGAII achieves the best results on 5 objectives, which also proves the effectiveness of maximin aggregation function in high-dimensional objective space. For DTLZ3, which is easily trapped in the local optimum and difficult to converge to the global optimum, MB-NSGAII achieves the best results on all three test problems except for the 2-objective. On DTLZ4, obtained after modifying the DTLZ2 function, MB-NSGAII achieves the best results on the 2-objective and 5-objective. Next two tables show the results of the algorithm under HV performance indicators. MB-NSGA-II has outstanding performance on most test instances. It can be seen from this that the validity of using the aggregation function to evaluate the populations. Of course, IBEA and NSGA-II also achieved the first in a few problems. The results in the table also prove that the strategy used in MB-NSGA-II is effective compared to NSGA-II.


### 4.2. Results and Discussion on Task Assignment Model

To confirm the ability of the algorithm to solve real-world problems, we built the robot assignment problem in the warehouse and solved it using the algorithm MB-NSGA-II. We complement the intelligent warehouse’s situation mentioned above by first setting a two-dimensional coordinate system with the coordinates of the shipping gate of the warehouse as (0,100) and the coordinates of the incoming gate of the warehouse as (100,0). The intelligent warehouse has a certain number of autonomous robots in perfect condition and equal shipping speed of 1 m/s. All tasks to be performed in the warehouse are generated randomly. For a uniform comparison, we set the population of 500 for the compared algorithms and the termination condition to 200 generations.To fully validate the performance of the algorithm in solving the task scheduling model, we assume that 5 or 10 robots are used to solve 100 randomly generated tasks; and 10 or 20 robots are used to solve 500 randomly generated tasks.For real multiobjective optimization problems, we usually use the comprehensive performance indicator HV to measure the merit of the algorithm. To use HV as an indicator, the maximum values on each objective are selected as a reference points. In this real-world MOP, we select the maximum value that can be reached on each objective separately, and the specific details of HV can be found in [30]. Tables 6 and 7 show the HV values that the algorithm can obtain on 100 and 500 tasks, respectively. From the tables, the results show that MB-NSGAII obtains the maximum HV on all these tasks. That is, the populations obtained by MB-NSGAII have better convergence and distribution. Through the above experiments, it can be proved that MB-NSGAII can give managers better solutions.Table 6
Experimental results of HV when robots perform 100 tasks.
Number of robotsIBEANSGA-IIMB-NSGAII50.13790.14260.1445100.25140.23410.2549Table 7
Experimental results of HV when robots perform 500 tasks.
Number of robotsIBEANSGA-IIMB-NSGAII100.13880.13740.1452200.18050.18040.1853Tables8 and 9 show the average time consumption of robots when solving 100 tasks with 5 robots and 500 tasks with 10 robots, respectively. Since each robot runs at a speed of 1 m/s, their average running distance can also be represented by Tables 8 and 9. As can be seen from the tables, the time consumption of each robot is close to each other and load balancing is achieved. Figure 6 shows the average time consumption of the robot when performing 100 tasks and 500 tasks.Table 8
Average time consumption of robots when using 5 robots for 100 tasks.
IBEANSGA-IIMB-NSGAIIRobot 13429.33358.33424.3Robot 23589.937053719.2Robot 33524.43520.33418.6Robot 43424.33103.93031.7Robot 53316.43301.33128Table 9
Average time consumption of robots when using 10 robots for 500 tasks.
IBEANSGA-IIMB-NSGAIIRobot 18894.59296.38227.2Robot 28677.28882.18693.7Robot 38516.47974.58681.3Robot 48094.87616.17594.8Robot 59015.481698170.3Robot 678787454.17684.7Robot 78511.77680.77364.1Robot 88244.38189.89170Robot 98665.49522.69233.9Robot 108739.782108620.1Figure 6
Schematic diagram of the average time consumed by robots to perform tasks.


## 4.1. Experimental Results on Benchmarks

This section provides the experimental results of MB-NSGAII with two other start-of-art multiobjective algorithms, NSGAII and IBEA, on the benchmark test suites DTLZ [28] and ZDT [29]. On DTLZ, each test instance was experimented on the 2,3,5- objectives, and on ZDT, the experiments were run on the 2-objective. Each algorithm needs to be run 20 times on each test issue. We can discuss the results of this experiment in terms of comprehensive performance indicator IGD.


### 4.1.1. Results and Discussion on ZDTs and DTLZs

The ZDT test suite contains six test problems with different characteristics, each test instance involving a feature that would cause the algorithm to have difficulty converging during optimization (e.g., multimodality). Through these test instances, it is possible to analyze which problems the algorithm excels in. DTLZ designs the problems through a systematic approach, and its decision variables and objectives are scalable, facilitating the algorithm to experiment in solving MaOPs. DTLZ increases the difficulty of the test problem by introducing manageable difficulties, facilitating the algorithm to experiment on test problems with any number of objectives. To more fully validate the effectiveness of several compared algorithms in solving MOPs, we chose to conduct experiments on 2,3, and 5 objectives.M in Table 1 indicates the objective number. The numbers of decision variable are set, as shown in Table 1. Details can be found in reference [28,29].Table 1
Setting of decision variables.
ProblemNumber of decision variablesZDT1-330ZDT4,610ZDT580DTLZ1M + 4DTLZ2-4M + 9To more visually reflect the ability of the algorithm on the benchmark test instances, we choose the comprehensive evaluation indicator HV and IGD as the basis for judging the performance of the algorithms. These two metrics are composed of different formulas, so the final population obtained by the algorithm can be evaluated from different perspectives. The formula and detailed description can be found in [30,31].


### 4.1.2. Experimental Settings

(i)
Population size and termination conditions: we use population of size 100 for DTLZ with 2,3 objectives and 200 for DTLZ with 5 objectives. For ZDTs, we select population of size 100 for the three algorithms. For DTLZ1-4 and ZDT1-6, the algorithms will end after 500 generations.(ii)
Crossover and mutation: in this experiment, binary crossover and polynomial mutation were used as operators for crossover and mutation. Set 1.0 as the crossover probability and set 1/D as the probability of mutation, where D is the number of decision variables.(iii)
Selection of indicator: for IBEA, we selectIε+ as a performance indicator to evaluate the strengths and weaknesses of individuals. The details of the indicator can be seen in [32].


### 4.1.3. Analysis of Experimental Results

Comparison data of MB-NSGA-II, NSGA-II, and IBEA on DTLZ and ZDT are presented in Tables2–5. These tables show the IGD and HV means and standard deviations of these compared algorithms on the benchmark test problems, and highlights the algorithm that worked best on these test problems. The Wilcoxon rank sum test with a significance level of 0.05 was used for the standard deviation analysis when conducting the experiment. The symbols “+,” “−,” and “=” are used to indicate that the algorithms are significantly better, significantly worse, or not significantly different between the compared algorithms and MB-NSGA-II.Table 2
Experimental results of IGD on ZDTs (mean and standard deviation).
ProblemMNSGAIIIBEAMBNSGAIIZDT124.7879e-3 (1.31e-4) -4.0901e-3 (1.01e-4) -3.9627e-3 (6.92e-5)ZDT224.9504e-3 (2.24e-4) +8.4212e-3 (7.74e-4) -5.3033e-3 (1.32e-4)ZDT325.4102e-3 (1.31e-4) -1.6020e-2 (8.05e-4) -5.0404e-3 (1.24e-4)ZDT425.0667e-3 (4.10e-4) +2.2761e-2 (1.59e-2) -6.3474e-3 (2.30e-4)ZDT525.1086e-1 (8.31e-2) +1.8513e+0 (4.39e-1)=1.6593e+0 (7.13e-1)ZDT623.7200e-3 (1.20e-4)=4.4425e-3 (1.14e-4) -3.7006e-3 (5.32e-5)Table 3
Experimental results of IGD on several DTLZs (mean and standard deviation).
ProblemMNSGAIIIBEAMBNSGAIIDTLZ122.2068e-3 (7.35e-5) -7.9965e-2 (7.99e-3) -2.0230e-3 (5.45e-5)32.8447e-2 (1.33e-3) -1.6616e-1 (2.51e-2) -2.2376e-2 (3.68e-4)51.7146e-1 (1.25e-1) -1.8372e-1 (2.16e-2) -5.1833e-2 (3.37e-4)DTLZ225.0483e-3 (1.75e-4) +1.6671e-2 (1.60e-3) -8.5238e-3 (2.92e-4)37.3058e-2 (2.57e-3) =8.4162e-2 (2.51e-3) -7.3973e-2 (1.57e-3)52.0563e-1 (4.67e-3) -1.9251e-1 (1.54e-3) -1.7906e-1 (2.44e-3)DTLZ327.1268e-3 (1.39e-3) +3.3913e-1 (8.73e-3) -9.8870e-3 (1.10e-3)31.4200e-1 (2.06e-1) =4.7891e-1 (4.51e-3) -7.8834e-2 (2.64e-3)57.5535e-1 (7.31e-1) -5.9365e-1 (9.76e-3) -1.8301e-1 (2.58e-3)DTLZ421.5241e-1 (3.11e-1) -4.5174e-1 (3.75e-1) -8.1849e-2 (2.32e-1)31.5811e-1 (2.77e-1) +8.2918e-2 (2.68e-3) =2.6157e-1 (2.41e-1)52.0630e-1 (3.41e-3) -2.1267e-1 (7.27e-2) -1.8232e-1 (1.93e-3)Table 4
Experimental results of HV on several DTLZs (mean and standard deviation).
ProblemMNSGAIIIBEAMBNSGAIIDTLZ125.8121e-1 (3.60e-4) -3.9969e-1 (1.87e-2) -5.8162e-1 (3.86e-4)38.2082e-1 (5.02e-3) -4.7837e-1 (5.50e-2) -8.3698e-1 (1.86e-3)56.9302e-1 (3.39e-1) -7.5172e-1 (4.51e-2) -9.7599e-1 (1.40e-3)DTLZ223.4654e-1 (1.97e-4) -3.4613e-1 (2.17e-4) -3.4731e-1 (1.09e-4)35.2551e-1 (3.99e-3) -5.5477e-1 (1.48e-3) -5.5689e-1 (9.61e-4)56.7753e-1 (7.73e-3) -8.0904e-1 (1.36e-3) +8.0024e-1 (2.49e-3)DTLZ323.4146e-1 (2.73e-3) =1.6879e-1 (3.86e-3) -3.4251e-1 (2.27e-3)34.5747e-1 (1.56e-1) -2.3645e-1 (9.09e-3) -5.4725e-1 (3.20e-3)54.0235e-1 (3.46e-1) -3.8077e-1 (6.39e-3) -8.0179e-1 (4.10e-3)DTLZ422.9554e-1 (1.08e-1) -1.9305e-1 (1.32e-1) -3.2164e-1 (8.11e-2)34.8619e-1 (1.39e-1) =5.5574e-1 (9.75e-4) =4.7330e-1 (1.09e-1)56.8256e-1 (8.53e-3) -8.0150e-1 (2.60e-2) =8.0780e-1 (1.73e-3)Table 5
Experimental results of HV on ZDTs (mean and standard deviation).
ProblemMNSGAIIIBEAMBNSGAIIZDT127.1925e-1 (1.76e-4) -7.2017e-1 (1.36e-4) -7.2031e-1 (8.51e-5)ZDT224.4399e-1 (2.03e-4) -4.4410e-1 (1.58e-4) -4.4485e-1 (6.70e-5)ZDT325.9938e-1 (8.34e-5) -5.9816e-1 (1.17e-4) -5.9973e-1 (4.53e-5)ZDT427.1800e-1 (1.09e-3) =7.0312e-1 (9.85e-3) -7.1735e-1 (5.54e-4)ZDT528.1676e-1 (1.20e-2) =8.1119e-1 (9.16e-3)8.1451e-1 (3.26e-4)ZDT623.8826e-1 (1.18e-4) =3.8766e-1 (1.17e-4) -3.8832e-1 (5.60e-5)In Table2 and Table 3, MB-NSGA-II has an outstanding performance on both DTLZ test problems and ZDT test problems, and it ranks first on 11 out of 18 test instances, and it is also competitive on other test instances. In particular, it achieved the best results on DTLZ1 for both 2,3,5 objectives, and also ranked first on DTLZ3 for two objective problems. Of course, NSGA-II also had good results, obtaining the best results on the six problems. In Table 3, MB-NSGAII ranked first in performance on 8 test problems, while NSGAII values achieved good results on three problems. It can be seen that MB-NSGA-II achieved the best results on all test problems on DTLZ1 with hyperplane PF, which is a significant improvement over NSGAII. On DTLZ2, MB-NSGAII achieves the best results on 5 objectives, which also proves the effectiveness of maximin aggregation function in high-dimensional objective space. For DTLZ3, which is easily trapped in the local optimum and difficult to converge to the global optimum, MB-NSGAII achieves the best results on all three test problems except for the 2-objective. On DTLZ4, obtained after modifying the DTLZ2 function, MB-NSGAII achieves the best results on the 2-objective and 5-objective. Next two tables show the results of the algorithm under HV performance indicators. MB-NSGA-II has outstanding performance on most test instances. It can be seen from this that the validity of using the aggregation function to evaluate the populations. Of course, IBEA and NSGA-II also achieved the first in a few problems. The results in the table also prove that the strategy used in MB-NSGA-II is effective compared to NSGA-II.


## 4.1.1. Results and Discussion on ZDTs and DTLZs

The ZDT test suite contains six test problems with different characteristics, each test instance involving a feature that would cause the algorithm to have difficulty converging during optimization (e.g., multimodality). Through these test instances, it is possible to analyze which problems the algorithm excels in. DTLZ designs the problems through a systematic approach, and its decision variables and objectives are scalable, facilitating the algorithm to experiment in solving MaOPs. DTLZ increases the difficulty of the test problem by introducing manageable difficulties, facilitating the algorithm to experiment on test problems with any number of objectives. To more fully validate the effectiveness of several compared algorithms in solving MOPs, we chose to conduct experiments on 2,3, and 5 objectives.M in Table 1 indicates the objective number. The numbers of decision variable are set, as shown in Table 1. Details can be found in reference [28,29].Table 1
Setting of decision variables.
ProblemNumber of decision variablesZDT1-330ZDT4,610ZDT580DTLZ1M + 4DTLZ2-4M + 9To more visually reflect the ability of the algorithm on the benchmark test instances, we choose the comprehensive evaluation indicator HV and IGD as the basis for judging the performance of the algorithms. These two metrics are composed of different formulas, so the final population obtained by the algorithm can be evaluated from different perspectives. The formula and detailed description can be found in [30,31].


## 4.1.2. Experimental Settings

(i)
Population size and termination conditions: we use population of size 100 for DTLZ with 2,3 objectives and 200 for DTLZ with 5 objectives. For ZDTs, we select population of size 100 for the three algorithms. For DTLZ1-4 and ZDT1-6, the algorithms will end after 500 generations.(ii)
Crossover and mutation: in this experiment, binary crossover and polynomial mutation were used as operators for crossover and mutation. Set 1.0 as the crossover probability and set 1/D as the probability of mutation, where D is the number of decision variables.(iii)
Selection of indicator: for IBEA, we selectIε+ as a performance indicator to evaluate the strengths and weaknesses of individuals. The details of the indicator can be seen in [32].


## 4.1.3. Analysis of Experimental Results

Comparison data of MB-NSGA-II, NSGA-II, and IBEA on DTLZ and ZDT are presented in Tables2–5. These tables show the IGD and HV means and standard deviations of these compared algorithms on the benchmark test problems, and highlights the algorithm that worked best on these test problems. The Wilcoxon rank sum test with a significance level of 0.05 was used for the standard deviation analysis when conducting the experiment. The symbols “+,” “−,” and “=” are used to indicate that the algorithms are significantly better, significantly worse, or not significantly different between the compared algorithms and MB-NSGA-II.Table 2
Experimental results of IGD on ZDTs (mean and standard deviation).
ProblemMNSGAIIIBEAMBNSGAIIZDT124.7879e-3 (1.31e-4) -4.0901e-3 (1.01e-4) -3.9627e-3 (6.92e-5)ZDT224.9504e-3 (2.24e-4) +8.4212e-3 (7.74e-4) -5.3033e-3 (1.32e-4)ZDT325.4102e-3 (1.31e-4) -1.6020e-2 (8.05e-4) -5.0404e-3 (1.24e-4)ZDT425.0667e-3 (4.10e-4) +2.2761e-2 (1.59e-2) -6.3474e-3 (2.30e-4)ZDT525.1086e-1 (8.31e-2) +1.8513e+0 (4.39e-1)=1.6593e+0 (7.13e-1)ZDT623.7200e-3 (1.20e-4)=4.4425e-3 (1.14e-4) -3.7006e-3 (5.32e-5)Table 3
Experimental results of IGD on several DTLZs (mean and standard deviation).
ProblemMNSGAIIIBEAMBNSGAIIDTLZ122.2068e-3 (7.35e-5) -7.9965e-2 (7.99e-3) -2.0230e-3 (5.45e-5)32.8447e-2 (1.33e-3) -1.6616e-1 (2.51e-2) -2.2376e-2 (3.68e-4)51.7146e-1 (1.25e-1) -1.8372e-1 (2.16e-2) -5.1833e-2 (3.37e-4)DTLZ225.0483e-3 (1.75e-4) +1.6671e-2 (1.60e-3) -8.5238e-3 (2.92e-4)37.3058e-2 (2.57e-3) =8.4162e-2 (2.51e-3) -7.3973e-2 (1.57e-3)52.0563e-1 (4.67e-3) -1.9251e-1 (1.54e-3) -1.7906e-1 (2.44e-3)DTLZ327.1268e-3 (1.39e-3) +3.3913e-1 (8.73e-3) -9.8870e-3 (1.10e-3)31.4200e-1 (2.06e-1) =4.7891e-1 (4.51e-3) -7.8834e-2 (2.64e-3)57.5535e-1 (7.31e-1) -5.9365e-1 (9.76e-3) -1.8301e-1 (2.58e-3)DTLZ421.5241e-1 (3.11e-1) -4.5174e-1 (3.75e-1) -8.1849e-2 (2.32e-1)31.5811e-1 (2.77e-1) +8.2918e-2 (2.68e-3) =2.6157e-1 (2.41e-1)52.0630e-1 (3.41e-3) -2.1267e-1 (7.27e-2) -1.8232e-1 (1.93e-3)Table 4
Experimental results of HV on several DTLZs (mean and standard deviation).
ProblemMNSGAIIIBEAMBNSGAIIDTLZ125.8121e-1 (3.60e-4) -3.9969e-1 (1.87e-2) -5.8162e-1 (3.86e-4)38.2082e-1 (5.02e-3) -4.7837e-1 (5.50e-2) -8.3698e-1 (1.86e-3)56.9302e-1 (3.39e-1) -7.5172e-1 (4.51e-2) -9.7599e-1 (1.40e-3)DTLZ223.4654e-1 (1.97e-4) -3.4613e-1 (2.17e-4) -3.4731e-1 (1.09e-4)35.2551e-1 (3.99e-3) -5.5477e-1 (1.48e-3) -5.5689e-1 (9.61e-4)56.7753e-1 (7.73e-3) -8.0904e-1 (1.36e-3) +8.0024e-1 (2.49e-3)DTLZ323.4146e-1 (2.73e-3) =1.6879e-1 (3.86e-3) -3.4251e-1 (2.27e-3)34.5747e-1 (1.56e-1) -2.3645e-1 (9.09e-3) -5.4725e-1 (3.20e-3)54.0235e-1 (3.46e-1) -3.8077e-1 (6.39e-3) -8.0179e-1 (4.10e-3)DTLZ422.9554e-1 (1.08e-1) -1.9305e-1 (1.32e-1) -3.2164e-1 (8.11e-2)34.8619e-1 (1.39e-1) =5.5574e-1 (9.75e-4) =4.7330e-1 (1.09e-1)56.8256e-1 (8.53e-3) -8.0150e-1 (2.60e-2) =8.0780e-1 (1.73e-3)Table 5
Experimental results of HV on ZDTs (mean and standard deviation).
ProblemMNSGAIIIBEAMBNSGAIIZDT127.1925e-1 (1.76e-4) -7.2017e-1 (1.36e-4) -7.2031e-1 (8.51e-5)ZDT224.4399e-1 (2.03e-4) -4.4410e-1 (1.58e-4) -4.4485e-1 (6.70e-5)ZDT325.9938e-1 (8.34e-5) -5.9816e-1 (1.17e-4) -5.9973e-1 (4.53e-5)ZDT427.1800e-1 (1.09e-3) =7.0312e-1 (9.85e-3) -7.1735e-1 (5.54e-4)ZDT528.1676e-1 (1.20e-2) =8.1119e-1 (9.16e-3)8.1451e-1 (3.26e-4)ZDT623.8826e-1 (1.18e-4) =3.8766e-1 (1.17e-4) -3.8832e-1 (5.60e-5)In Table2 and Table 3, MB-NSGA-II has an outstanding performance on both DTLZ test problems and ZDT test problems, and it ranks first on 11 out of 18 test instances, and it is also competitive on other test instances. In particular, it achieved the best results on DTLZ1 for both 2,3,5 objectives, and also ranked first on DTLZ3 for two objective problems. Of course, NSGA-II also had good results, obtaining the best results on the six problems. In Table 3, MB-NSGAII ranked first in performance on 8 test problems, while NSGAII values achieved good results on three problems. It can be seen that MB-NSGA-II achieved the best results on all test problems on DTLZ1 with hyperplane PF, which is a significant improvement over NSGAII. On DTLZ2, MB-NSGAII achieves the best results on 5 objectives, which also proves the effectiveness of maximin aggregation function in high-dimensional objective space. For DTLZ3, which is easily trapped in the local optimum and difficult to converge to the global optimum, MB-NSGAII achieves the best results on all three test problems except for the 2-objective. On DTLZ4, obtained after modifying the DTLZ2 function, MB-NSGAII achieves the best results on the 2-objective and 5-objective. Next two tables show the results of the algorithm under HV performance indicators. MB-NSGA-II has outstanding performance on most test instances. It can be seen from this that the validity of using the aggregation function to evaluate the populations. Of course, IBEA and NSGA-II also achieved the first in a few problems. The results in the table also prove that the strategy used in MB-NSGA-II is effective compared to NSGA-II.


## 4.2. Results and Discussion on Task Assignment Model

To confirm the ability of the algorithm to solve real-world problems, we built the robot assignment problem in the warehouse and solved it using the algorithm MB-NSGA-II. We complement the intelligent warehouse’s situation mentioned above by first setting a two-dimensional coordinate system with the coordinates of the shipping gate of the warehouse as (0,100) and the coordinates of the incoming gate of the warehouse as (100,0). The intelligent warehouse has a certain number of autonomous robots in perfect condition and equal shipping speed of 1 m/s. All tasks to be performed in the warehouse are generated randomly. For a uniform comparison, we set the population of 500 for the compared algorithms and the termination condition to 200 generations.To fully validate the performance of the algorithm in solving the task scheduling model, we assume that 5 or 10 robots are used to solve 100 randomly generated tasks; and 10 or 20 robots are used to solve 500 randomly generated tasks.For real multiobjective optimization problems, we usually use the comprehensive performance indicator HV to measure the merit of the algorithm. To use HV as an indicator, the maximum values on each objective are selected as a reference points. In this real-world MOP, we select the maximum value that can be reached on each objective separately, and the specific details of HV can be found in [30]. Tables 6 and 7 show the HV values that the algorithm can obtain on 100 and 500 tasks, respectively. From the tables, the results show that MB-NSGAII obtains the maximum HV on all these tasks. That is, the populations obtained by MB-NSGAII have better convergence and distribution. Through the above experiments, it can be proved that MB-NSGAII can give managers better solutions.Table 6
Experimental results of HV when robots perform 100 tasks.
Number of robotsIBEANSGA-IIMB-NSGAII50.13790.14260.1445100.25140.23410.2549Table 7
Experimental results of HV when robots perform 500 tasks.
Number of robotsIBEANSGA-IIMB-NSGAII100.13880.13740.1452200.18050.18040.1853Tables8 and 9 show the average time consumption of robots when solving 100 tasks with 5 robots and 500 tasks with 10 robots, respectively. Since each robot runs at a speed of 1 m/s, their average running distance can also be represented by Tables 8 and 9. As can be seen from the tables, the time consumption of each robot is close to each other and load balancing is achieved. Figure 6 shows the average time consumption of the robot when performing 100 tasks and 500 tasks.Table 8
Average time consumption of robots when using 5 robots for 100 tasks.
IBEANSGA-IIMB-NSGAIIRobot 13429.33358.33424.3Robot 23589.937053719.2Robot 33524.43520.33418.6Robot 43424.33103.93031.7Robot 53316.43301.33128Table 9
Average time consumption of robots when using 10 robots for 500 tasks.
IBEANSGA-IIMB-NSGAIIRobot 18894.59296.38227.2Robot 28677.28882.18693.7Robot 38516.47974.58681.3Robot 48094.87616.17594.8Robot 59015.481698170.3Robot 678787454.17684.7Robot 78511.77680.77364.1Robot 88244.38189.89170Robot 98665.49522.69233.9Robot 108739.782108620.1Figure 6
Schematic diagram of the average time consumed by robots to perform tasks.


## 5. Conclusion

In this paper, we have proposed a corresponding mathematical model for robot task assignment problem in an intelligent warehouse and solve it by MB-NSGAII, which uses nondominated sorting, maximin aggregation function, and brain storm operator. MB-NSGA-II is validated on the DTLZ and ZDT test suites, and it achieves satisfactory results in solving the real-world MOP with robot task assignment.Next, we will continue to investigate the use of maximin aggregation function and brain storm operator to solve MOPs. We will also use this approach to solve more real-world MOPs.



---

*Source: 2021535-2022-06-15.xml*