Effective guide for behaviour of farmers in the withdrawal of rural homesteads: An evolutionary game-based study

1.   Introduction

The multi-point dynamic aggregation problem (MPDAP) ^[1] is widely used in resource allocation and robot routing scheduling problems, such as post-disaster rescue ^[2], e-commerce logistics, forest fire fighting ^[3], medical equipment scheduling ^[4], and other task allocation optimization problems ^[5,6]. As a representative of the problems, robot routing optimization with dynamic task scheduling is an interesting research issue in intelligent optimization community. For the purpose of improving the efficiency of completing tasks, the moving routes of the robots may often be adjusted according to the observed tasks. Another characteristic of the problem is that there are multiple robots involved, that is to say, this is a multi-robot task scheduling problem. Some interesting real-world applications, such as e-commerce robots ^[5] and agricultural robots ^[6], have been presented in which the multi-robot task allocation models are developed and solved. In these problems, one needs to consider that each may change over time, the path of each robot, and the cooperation with other robots.

To better illustrate the problem, a two-robot and five-task case is present in Figure 1 and Figure 2. Figure 1 is a schematic diagram of two robots cooperating to execute five tasks according to their own driving path from the depot, in which the rectangle is the depot, the black circle represents the robot, and the white circle stands for the task. The solid line is the driving path of robot 1, and the dashed line is the driving path of robot 2. Figure 2 shows the time-demand change diagram of task 3 which needs to be executed by the robots. In Figure 2, point $A$ gives the initial state value of task 3, point $B$ provides the total task demand facing robot 1 and robot 2 when they reach task point 3, and the task demand at point $C$ approaches 0, indicating that task 3 has been cooperatively completed by both robot 1 and robot 2. It can also be seen from Figure 2 that the cumulative time from the depot to completing task 3 is 5.3. One methodology to simplify the problem is to minimize the time cost at which all five tasks are completed by the two robots.

The above-mentioned multi-point dynamic aggregation problem is a combination of path optimization problems and the task assignment of multi-robot systems. It not only needs to plan the path of each robot to execute the task, but also reasonably adjust the cooperation between robots to complete the tasks at which the time-demanding amount may be dynamically changing. Therefore, the multi-point dynamic aggregation problem has the following challenges. First, in order to complete a task requires the cooperation of robots, each robot needs to consider its own path as well as plans of other robots during the task execution. Second, the time each robot requires to complete a certain task depends on the growth rate of task demand at each task, which is a dynamically decision-making procedure. Third, the optimization of both the traveling routes of robots and other objectives has to be made in a multi-factor interactive environment involving dynamic task allocation and cooperation among multiple individuals, etc.

Despite the computational complexity of the problem, the wide range of practical applications has stimulated the study of modeling and solution methods of the problem. In ^[7], a model of forest fire suppression was presented as a multiple permutation combinatorial optimization problem, named MPDAP, and a distribution estimation algorithm using an effective encoding method was proposed for this problem. The given experimental show that the proposed algorithm can find a better solution for MPDAP than other compared approaches. In ^[8], a nonlinear agent routing problem in multi-point dynamic task (ARP-MPDT) was developed, which was based on a node histogram model (NHM), an edge histogram model (EHM) and probabilistic modeling distribution estimation algorithm (EDA), in which the ratio of NHM and EHM probability models can be adjusted adaptively. Gao et al. ^[9] represented each robot task by multiple row access sequence coding which designed a kind of heuristic rule of implicit decoding strategy to simplify the problem from the mixed-discrete variable optimization, and was based on a one-step and multi-step local search and presented a memetic algorithm for solving the problem. Hao et al. ^[10] combined the advantages of differential evolution algorithms and distribution estimation algorithms, and then proposed a hybrid algorithm. The proposed algorithm was executed on MPDAP examples of different sizes, and the simulation results show that compared with single differential evolution or distribution estimation algorithms, the hybrid method can solve this problem more efficiently. In an ant colony framework, in order to improve the efficiency of the solution construction, a new pheromone matrix and pheromone updating mechanism were proposed ^[11]. In order to reduce the search space and delete some bad solutions, Gao et al. ^[11] designed a heuristic value evaluation function in solution construction, and a pheromone repair based mechanism was investigated to enhance the ability to build viable solutions, and a local search was adopted to improve the balance between exploration and exploitation.

It should be noted that most of existing models for MPDAPs work to minimize the completion time of all tasks, in which the number of robots is predetermined. A recurring question is "how may robots are best?", that is to say, from the point of view of controlling total costs involving robots, time costs and material damage, how does a decision-maker optimize the number of dispatching robots and the traveling lines of them? Different from the above literatures, in the manuscript, a new bilevel programming problem is provided and solved by designing a new encoding scheme and efficient evolutionary operators. Some main innovations of this manuscript are presented as follows:

1) In the proposed model, the total cost of performing all tasks, instead of only time cost, is taken into account, in which the number of robots need to be optimized. Hence, the existing model with a predetermined number of robots is simply a special case of the proposed model. In the proposed model, both the robot cost and completion time are located in different levels and hierarchically optimized.

2) A new matrix encoding scheme is developed for all follower's solutions, which is simpler in form and more efficient than most of the existing encoding methods.

3) Both novel two-strategy mutation and selection schemes are designed and adopted in different evolution stages, which helps to find better potential solutions.

The rest of the manuscript is arranged as follows. Section 2 introduces bilevel models and differential evolution methods. Section 3 introduces the MPDAP model. The detailed algorithm design is presented in Section 4. The experimental studies are implemented in Section 5. Section 6 provides a discussion. Finally, Section 7 concludes this manuscript.

2.   Preliminaries

2.1.   Bilevel optimization models

Bilevel programming is a hierarchical optimization problem, in which variables are divided and located at two levels, the leader's and follower's levels, respectively. A general bilevel optimization model can be formulated as

where

and

are called the leader's and follower's problems, respectively. $F(x, y)$ and $f(x, y)$ are the leader's and follower's objective functions, respectively, $G(x, y)$ and $g(x, y)$ are the constraints of the leader's and follower's problems, respectively, and, $x$ and $y$ are called the leader's and follower's variables, respectively. The bilevel programming problem is widely used in many practical fields, such as economics, engineering science, operation research, and optimal control ^[12], and has a decision-making system with a hierarchical (nested) optimization process, in which multiple decision-makers with different priorities locates at different decision-making levels. Generally speaking, the priority of decision-makers at the higher level is higher than that of decision-makers at the lower level. The decision-makers at the higher level regulate the decision-making process at the lower level through optimization procedures, and the decision of the follower can also affect the optimization result of the leader's problem.

The computational difficulty of solving the bilevel optimization problem lies in the large amount of computation caused by frequently solving the follower's problems. At present, the methods to solve bilevel optimization problems mainly include the traditional optimization method using gradient descent, swarm intelligent approaches, and hybrid technologies. The traditional optimization methods include the single level reduction method, the branch and bound method, the penalty function method, the trust region method, and so on. Sinha et al. ^[13] used Karush-Kuhn-Tucher (K-K-T) conditions to transform the bilevel optimization model into a single-level optimization problem, though the transformed model always causes a complex constraint restriction. For integer bilevel optimization problems, Liu et al. ^[14] proposed a branch and bound method based on a relax and check procedure. Li et al. ^[15] proposed a numerical method for solving bilevel optimization problems with non-convex follower's problems, in which, the follower's problem is replaced by an equivalent substitution through K-K-T conditions. Joseph et al. ^[16] proposed a bilevel optimization model for classification feature selection problems and designed a genetic algorithm based on derivatively-free optimization. Li et al. ^[17] studied a class of linear fraction bilevel optimization problems where the coefficients and the right vector of the objective function are interval numbers and proposed an evolutionary algorithm based on optimality conditions. Aboelnage et al. ^[18] proposed an improved algorithm based on a new selection method and a chaotic search method for bilevel optimization problems. Goshu et al. ^[19] proposed a meta-heuristic algorithm based on a random space system sampling technique to solve the stochastic bilevel optimization problem. The hybrid optimization method combines the advantages of various algorithms to deal with the bilevel optimization problem, and shows a good performance. Islam et al. ^[20] proposed a memetic algorithm combining global search and local search. Yousria et al. ^[21] proposed a hybrid algorithm combining an improved genetic algorithm and chaotic search techniques. Sinha et al. ^[22] proposed a bilevel optimization algorithm with a multi-valued iterative approximation, which uses the follower's reaction set mapping to reduce the computation. Molai ^[23] used an extension principle for fuzzy multi-objective linear bi-level optimization problems.

Most of the existing methods are efficient in solving bilevel programming problems with small scale convex and differentiable functions. However, when the scale of the problem is sharply increased and non-convex, either non-differential functions or discrete variables are involved and the performance of these approaches always degenerates rapidly and can't be used to deal with the problem of this kind.

2.2.   Differential evolution

Differential evolution (DE) is a very popular evolutionary algorithm paradigm proposed by Storn and Price ^[24]. It contains four processes: initialization, mutation, crossover, and selection.

Initialization: There are $N$ individuals generated randomly in a feasible region.

where $X_m$ is the $m$th individual.

Mutation: The mutation operation is the main way for populations to generate new individuals, which can increase population diversity. For each individual $X_m$, a mutation vector $V_m$ is generated. Some common mutation operators are as follows

DE/rand/1/bin:

DE/best/1/bin:

DE/rand/2/bin:

DE/best/2/bin:

DE/rand-to-best/1/bin:

where $r_1$, $r_2$, $r_3$, $r_4$, and $r_5$ are five different integers chosen randomly from $[1, N]$, $X_{best}$ is the best individual in the population, $rand(0, 1)$ is a uniformly distributed random number from $[0, 1]$, and $F_{sf}$, $F_{sf_1}$, and $F_{sf_2}$ are scaling factors.

Crossover: Using binomial crossover to generate crossover offspring $U_m$.

Selection: The better one between $X_m$ and $U_m$ is selected into the next generation.

DE uses the differential information of multiple individuals in a population to achieve perturbation of evolved individuals. As s one of the effective swarm intelligent algorithms, DE has been successfully applied to many practical problems ^[25]. Wang et al. ^[26] used a new encoding method based differential evolution to solve the distribution problem of wind farms. Aiming at power optimization problem of power system, Chi et al. ^[27] adopted a differential evolution algorithm based on different evolutionary operations. Additionally, other interesting applications can be found, such as economic or emission dispatch problem ^[28], image segmentation ^[29], and task scheduling in cloud computing environment ^[30]. As a bilevel programming model with discrete variables, the proposed MPDAP in the manuscript is very hard to solve. In view of the good performance of differential evolution in dealing with complex optimization problems, this paper developed an improved version of this algorithm to solve the proposed MPDAP model.

3.   Multi-point dynamic aggregation problem

3.1.   Demand state model for task in MPDAP

From the introduction, it can be seen that in a MPDAP, the task demand of each task (e.g. fire point) may dynamically increase over time. The relationship between task demand and time is suggested by the following equations in ^[11]:

where, $q_{i}^t$ represents the demand at task $i$ at time t, the initial demand of task $i$ is $q_{i}^0$, and $\alpha_i$ is the inherent increment rate of task $i$.

3.2.   The model of the MPDAP

The workload of the robot $r$ to solve the task $i$ can be expressed by the following formula:

where, $S_{{task}_i}$ represents the total amount of work that robot does, $v_{task}^r$ is the work efficiency, and $t_{{task}_i}$ stands for the time.

$S_{path_{ji}}$ is the path length from task $j$ to task $i$, the speed of the robot $r$ is $v_{path}^r$, and the time taken from task $j$ to task $i$ is $t_{path_{ji}}$. The formula is as follows:

In ^[11], the objective function in the MPDAP model is

where, $t_i$ is the completion time of task $i$. In fact, in a real-world planing problem, it is very important to control overall cost, not just time cost. In this manuscript, a more general case is taken into account, which involves the cost of operation time, damage caused by fire, and robot cost. In the model, once the number of robots are determined, one can program the route and assign a task for each robot. In this procedure, a hierarchical optimization process arises. As a result, a bilevel programming model is developed as follows:

where, the leader's objective function $F$ stands for the cost of all the robots to complete the tasks and the amount of damage in the tasks; the follower's objection function $f$ is the maximum completion time of all the tasks.

Equation (3.6) represents the total amount of damages at all $M_1$ tasks when all tasks are completed, in which the total amount of damage at task $i$ is $\lambda_1$ times the task demand at task $i$. Equation (3.7) is the state-determined cost of renting (or buying) $M_2$ robots used in the tasks, here, $c_r^0$ is the inherent cost of robot $r$. Equation (3.8) is the running cost of all the robots and $c_r$ is the unit running cost of robot $r$.

As is known, model (3.4) is simply the follower's problem of model (3.5). The number of robots and inherent costs in model (3.5) are fixed, which is degenerated to model (3.4). Therefore, model (3.5) is more general and extensive than model (3.4).

In order to make the model of multiple robots and tasks more explicit, relevant assumptions in ^[11] are also adopted in this manuscript:

1) There are $M_2$ robots, and their initial location is the depot. The cost of robots of the same size is the same, and they have the same traveling speed $v_{path}$ and the work efficiency $v_{task}$.

2) There are $M_1$ tasks in different positions, the tasks at different locations may have the same or different task demand.

3) In the working environment, there are no obstacles, and each robot does not have to consider the collision with others. Once the driving path is determined, the robot can drive.

4) The starting time that the robots leaves the depot is denoted as 0.

5) For each task, the number of ways into the task point equals to that of ways out.

6) Each task is performed by at least one robot.

7) Each task is executed by a robot at most once.

8) Once the current task is completed, the robot can move on to the next.

9) The moving path of each robot can be planned in advance before performing tasks.

4.   Algorithmic design

In the proposed bilevel model (3.5), the leader's problem is to optimize the number of robots assigned to tasks. Once the number of the robots is determined, the follower needs to optimize the moving routes of robots. Since the optimization procedure is hierarchical and most of the computational amounts mainly come from the solutions of the follower's problem, the manuscript is first focused on the follower's algorithm design. Finally, by embedding the follower's algorithm, a novel bilevel differential evolution (NDE) is presented for the bilevel model (3.5).

4.1.   Algorithmic design of the follower's problem

This section mainly introduces the main schemes, as well as the follower's algorithm. First of all, a new encoding scheme is developed, which can more efficiently express the robot's traveling path and completion tasks than others. In addition, a two-strategy mutation is provided by using the mean of locations of the top individuals, which can efficiently avoid the directional mislead from a single good point.

4.1.1.   Encoding / Decoding strategies

In an MPDAP tasks planning problem with $M_1$ tasks and $M_2$ robots, a feasible solution is encoded as an $M_2\times (M_1+1)$ matrix $X$:

where the elements in the first column of the matrix (4.1) are all 1, indicating that every robot starts from the depot. From the second column to the last column, the elements are 0 or 1, as shown in (4.2). (4.3) indicates that, for each task, the total ability of the robots executing it must be greater than its inherent increment rate. Otherwise, the task can never be completed.

The matrix encoding strategy is presented in a simpler form, where the elements of the matrix are either 0 or 1, that is to say, there are only two states, go or not. In the encoding scheme, the driving path of each robot is represented by a row, and a feasible solution is equivalent to a scheduling plan that can efficiently complete all tasks. Furthermore, a binary encoding structure can be conveniently used in evolutionary operations and can efficiently search for potential path combinations.

The decoding procedure is executed as follows: in each row of (4.1), value 1 means the robot visits the task, whereas value 0 represents the robot never visits the task. The order of visiting tasks is determined by

where $\alpha_i$ is the inherent increment rate of task $i$ and $s_{path_{i}}$ is the path length from depot to task $i$. A task with a larger $\alpha_i$ and a smaller $s_{path_{i}}$ should be visited preferentially, Hence, $VO_{task_i}$ can be taken as an index by which the tasks can be chosen by a probability choice scheme, one by one.

For example, the example in Figure 1 would be encoded as follows:

where the first row $[1 \; 1 \; 0 \; 1 \; 1 \; 0]$ means that tasks $[depot, task_1, task_3, task_4]$ needs to be visited sequentially. The values of the relevant parameters of the tasks in Table 1. According to the order $VO_{task_3} > VO_{task_1} > VO_{task_4}$, the visiting path of $robot_1$ is determined as $depot\rightarrow task_3\rightarrow task_1\rightarrow task_4$.

4.1.2.   Evolutionary operators

Ⅰ. Mutation operator

In the proposed follower's algorithm, the total two mutation strategies are adopted in early and later evolution stages, respectively. The first mutation operator is designed as follows:

Mutation operator 1, $DE/rand-to-mean_{better}/1/bin$:

In the mutation operator, a group of top individuals are chosen to provide an evolutionary direction for individuals. Different from the mutation operator 2 (as below), the mutation operator 1 can keep the diversity of offspring since the heuristic information comes from a group of better individuals, instead of a single point. The advantage of the scheme is that more than one better point can provide a more stable evolutionary direction than that a single point does.

In the later stage of evolution, the population trends to convergence and then only a single best individual can work well. As a result, the following mutation operation 2 is directly adopted.

Mutation operation 2, $DE/rand-to-best/1/bin$:

where $X_{r_1}$, $X_{r_2}$, and $X_{r_3}$ are randomly selected individuals in the population, $rand(0, 1)$ represents a random number that follows an uniform distribution in $[0, 1]$, $F_{sf}\in [0, 2]$ is the scaling factor, $X_{mean_{best}}$ is the average position of some top individuals, and $X_{best}$ is the best individual found so far.

The offspring of the mutation operator can be described as

Ⅱ. Crossover operator

Through the binomial crossover, a trial vector $U_m$ is generated based on $X_m$ and $V_m$.

where $m = 1, 2, \cdots, NP$, $CR \in [0, 1]$ is the crossover control parameter, and $rand_j$ is a random number $[0, 1]$.

Ⅲ. Individual repair operator

During the evolution of the individual, there may be cases where none of the robots can complete their allocated tasks, leading to an infeasible solution. For example, given an MPDAP instance with two robots and five tasks, the inherent increment rate of a task is larger than the individual ability of both robots. Thus, the robots have to go to the same task and complete it. If only a single robot is asked to execute a task with a large inherent increment rate separately, then such a task will not be completed. To address this issue, a individual repair operator is developed to reallocate the robots to make the solution feasible. The mutation and crossover operators only affect the second to last column of the individual matrix $X_m$ and $V_m$, so the individual repair operator will only target these columns as well.

In the offspring individual matrix $U_m$, the elements in the first column are all 1, but the elements in the other positions may not take values of 0 or 1, and perhaps do not satisfy the constraint (4.3). The repair process is as follows:

Step 1. The elements in the individual matrix are made feasible, i.e., the elements greater than 1 become 1, the elements less than 0 become 0, and others are rounded off.

Step 2. Verify whether the elements of each column satisfy constraints (4.3), if so, turn to step 4, otherwise, turn to step 3.

Step 3. Randomly select the element in each column that is 0 and change it to 1 until the element in that column satisfies constraint (4.3).

Step 4. End the repair process.

Ⅳ. Selection operator

In the process of evolution, two selection operators are used. One is the same as in the original differential evolution, which is denoted by Selection operator 1 and expressed as follows:

where, $f_{fit}$, as the fitness, is the follower's objective function.

The other is newly provided as a selection operator 2, which can be finished as follows: in all individuals, parents and offspring, some top individuals are first selected into the next generation of population, and then the rest are chosen by the procedure of selection operator 1. It follows that the selection operator 2 focuses on better individuals instead of diversity kept by parents. Hence the selection operator 2 is helpful to speed up the convergence of the algorithm.

In early evolution, the selection operator 1 is executed for the purpose of maintaining diversity. In order to improve the exploitation capability of the algorithm, the selection operator 2 is utilized in the later stage.

4.1.3.   Pseudo-code and flowchart of the follower's algorithm

Based on the mentioned-above design, the pseudo-code of the follower's algorithm is presented in Algorithm 1.

4.2.   A novel differential evolution for the bilevel model (NDE)

In the bilevel optimization algorithm for solving the MPDAP, the leader's objective function is taken as the fitness function, and a classical differential evolution is adopted to optimize the variables of the leader's problem. The pseudo-code of the bilevel optimization algorithm is presented in Algorithm 2.

5.   Simulation

5.1.   Comparison with other methods

All 50 benchmark instances are taken from ^[11]. In these 50 benchmark instances, the number of robots is fixed, and the task demand function is a linear function. The differences between groups of benchmark instances are shown in Table 2. $M_1$ is the number of tasks and $M_2$ is the number of robots. In order to make the experimental results well presented, the benchmark instances were divided into 3 groups, namely Group 1, Group 2, and Group 3. All the instances are called by their group name, number of robots, number of tasks, and the ratio between the sum of all the task inherent rates and the sum of all the robot abilities. For example, the benchmark instance 1 is named as $G_1$$\_$5$\_$4$\_$0.39, where $G_1$ denotes the group of the instance, 5 is the number of robots, 4 is the number of tasks, and 0.39 is the ratio. In order to make the comparison fair and verify the effectiveness of the NDE, the parameter settings are taken similar to the benchmark instances, and the stopping criterion is set is the same as the competitors. The proposed algorithm is run independently for 30 times, and the computational results are recorded in Tables 3-8. The data of MA ^[9] based on two variants are shown in columns MA-MLS and MA-OLS in Tables 3, 5, and 7. The data of EDA ^[10] are shown in column EDA in Tables 3, 5, and 7. The data of ILS ^[32] are presented in column ILS in Tables 4, 6, and 8. The data of AC-ACO ^[11] are provided in column AC-ACO in Table 4, 6, and 8. Column NDE provides the computational results by NDE.

In Tables 3-8, the comparison data include the means and standard deviations of the completion time of tasks, and $*$ means the corresponding method cannot obtain a feasible solution in a limited period of time. The numbers in parentheses indicate the rank-order of all compared approaches on a benchmark instance. The smaller the rank value, the better the algorithm. In addition, some visual comparisons are also provided by Figure 3 (the Group 1), Figure 4 (the Group 2), and Figure 5 (the Group 3). For the case that some algorithms fail to find the optimal solution on some instances, a maximum value is used on the figure to ensure the visuality of figures.

The results of the Group 1 are shown in Table 3 and Table 4. NDE is slightly worse than MA-MLS in $G_1$$\_$15$\_$20$\_$0.67, $G_1$$\_$17$\_$23$\_$1.71, $G_1$$\_$20$\_$20$\_$0.58, and $G_1$$\_$20$\_$20$\_$0.967. However, NDE is better than the comparison algorithm in the remaining 23 instances in Table 3. In Table 4, NDE provides worse results for 11 instances, but better results for the remaining 16 instances. Among the 27 benchmark instances in Group 1, the rank-order values of NDE are 31 and 41, respectively. In Table 4, the rank-order values of both AC-ACO and NDE are 41, indicating that they both have similar algorithm performances. Compared with other algorithms except AC-ACO, the rank values of NDE are the smallest.

The computational results on the Group 2 are shown in Table 5 and Table 6. In Table 5-6, NDE has slightly worse results for $G_2$$\_$40$\_$15$\_$0.67, but better results in the remaining 5 instances. For 6 benchmark instances in the Group 2, the rank-order of the NDE are 9 and 8, respectively. Compared with other algorithms, NDE has a better performance in Group 2.

The computational results on Group 3 are shown in Table 7 and Table 8. In Table 7, the results of NDE are all better than those of the comparison algorithm. Although Table 8 displays that NDE is slightly worse than AC-ACO in $G_3$$\_$60$\_$30$\_$1.14 and $G_3$$\_$80$\_$80$\_$0.54, NDE is better than the comparison algorithm in other 15 instances. For 17 benchmark instances in Group 3, the rank-orders of the NDE are 17 and 19, respectively. Compared with other algorithms, NDE has a better performance in Group 3.

NDE and AC-ACO have similar performance in solving benchmark instances in the Group 1; however NDE has better performance in solving benchmark instances in the Group 2 and the Group 3. Therefore, it is not difficult to see that the NDE is effective for solving these 50 benchmark instances.

5.2.   Experimental results under the bilevel optimization model

In this part, when the bilevel optimization model is taken into account, it is very necessary to analyse the effect on the total cost caused by the amount of robots. For the 50 benchmark instances in subsection 5.1., besides the same amount of robots as in literature, some additional robots are added to perform these tasks. The purpose of the procedure is to understand the optimal amount of robots, such that the total cost can be minimized. When robots need to be added, the robot with the lowest efficiency is first chosen since it can lead to the worse case, increase the cost, and lower the efficiency. The experimental results are shown in Tables 9 to 11. In Tables 9 to 11, $F$ represents the leader's objective function in the bilevel programming model, $f$ represents the follower's objective function, and $N_i$ represents the number of robots added to execute tasks.

Since there are multiple robots with different capabilities in the benchmark instances, in line with the working capability of robots, the state-determined cost of renting (or buying) robots is set as the value that equals to $\lambda_2$ times the working efficiency of the robots, and the running cost of the robot is the product of working efficiency and working time. Without any loss of generality, the values of both $\lambda_1$ and $\lambda_2$ are taken as 1 in the experiment.

Among the 27 benchmark instances in Group 1 in Table 9, based on the bilevel optimization model, shows that the optimal solution of 9 instances is better than that of the original instances. For example, for $G_1$$\_$5$\_$5$\_$1.66, the optimal solution needs to add 3 robots to the original instances. When the number of robots is increased, the state-determined cost of robots increases, but the time to complete the task can evidently decrease. As a result, the leader's objective function is decreased. So the decision maker can choose the better option.

For the 6 instances in Group 2 in Table 10, we find that on two instances ($G_2$$\_$20$\_$40$\_$3.61 and $G_2$$\_$15$\_$60$\_$3.64), the better solutions are obtained when an appropriate amount of robots are added. For the benchmark instances in Group 3, from Table 11, one can find that one instance $G_2$$\_$20$\_$80$\_$4.33 finds the better solution than original case.

6.   Discussion

In the field of task scheduling and route optimization, the multi-point dynamic aggregation problem is interesting but computationally complex. In the present research, in order to reduce decision-making cost, some key parameters are predetermined, which should be further optimized, such as the number of robots as well as tasks. Based on these assumptions, some simplified models have be given, which simply focus on optimizing the routes of robots to minimize the completion time of all tasks. In fact, from the decision maker's point of view, at least two additional costs, robots and property damage, need to be taken into account. It follows that optimization procedure should be divided into two levels. One is to determine the number of robots(leader's level), whereas the other is to optimize the routes of these robots(follower's level). As a result, a bilvel programming problem is proposed in the manuscript. When the leader's variables are given, the problem can degrade into a single-level problem just as in literature. Hence, the proposed bilevel problem is more universal than the existing models.

Bilevel programming problems are computationally hard. In order to efficiently deal with the hierarchical optimization problem, a two-level DE is developed. In the algorithm for the follower's program, a matrix encoding scheme is provided and followed by a newly designed decoding rule. The individual representation is in nature binary, which is simpler and more efficient than ever. In addition, the mean location of multiple high-qualify individuals can provide a stabler evolutionary direction than a single elite individual.

In spite of the fact that the bilevel model is more suitable for practice, because of the discrete variables and hierarchical structure involved, this problem is still very difficult to solve. More specifically, for large-scale cases in which too many robots and tasks are involved, the proposed algorithm still needs to pay expensive computation cost.

7.   Conclusion

In the real-world field, MPDAP is a challenging and difficult optimization problem. In this manuscript, the problem is re-modeled as a bilevel programming model by considering the total cost of completion tasks and the damages. In order to efficiently deal with the bilevel programming problem, a novel differential evolution with a nested structure is developed, in which a new matrix encoding method is designed to represent the routes of robots. Meanwhile, in order to solve the follower's problem efficiently, we also present two-strategy mutation and selection operators. From the experimental results, one can find that the proposed algorithm has better performance when compared with other similar algorithms. Also, based on the bilevel optimization model, the solutions with smaller cost are obtained by dispatching additional robots.

When some interesting topics, such as the dynamic task number and the priority of tasks, are taken into the model, the optimization procedure may become more complex than ever. In the future study, more efficient heuristical operators as well as optimization techniques need to be further developed for this kind of problems.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work is supported by the National Science Foundation of China (No.61966030).

Conflict of interest

The authors declare there is no conflict of interest.