Multirobot task allocation based on an improved particle swarm optimization approach

Abstract

Due to its complexity and non-deterministic polynomial-time hard characteristic, multirobot task allocation problem remains a challenging issue in the field of cooperative robotics. Thanks to its easy implementation and promising convergence speed, the particle swarm optimization method has recently aroused increasing research interest in the area of multirobot task allocation problem. However, the efficiency of the standard particle swarm optimization is hindered by several deficiencies such as the inefficient capabilities in balancing exploration and exploitation, as well as the high likelihood of plunging into stagnation. Aiming at enhancing the performance of particle swarm optimization via remedying these two drawbacks, this paper proposes an improved particle swarm optimization method, which integrates standard particle swarm optimization 2011 with evolutionary game theory. To prevent particles being locked into stagnation, particles in the proposed particle swarm optimization first adopt the updating rules of standard particle swarm optimization 2011 to undertake their movements. Subsequently, attempting to well trade off the exploration and exploitation capabilities of particles, a novel self-adaptive strategy, which is determined by the evolutionary stable strategies of evolutionary game theory and the iteration number of particle swarm optimization, is presented to adaptively adjust the main control parameters of particles in the proposed particle swarm optimization. Since the convergence of particle swarm optimization remains paramount and dramatically affects the performance of particle swarm optimization, this paper also analytically investigates the convergence of the proposed particle swarm optimization and provides a convergence-guaranteed parameter selection principle for the proposed method. Finally, leveraging the development of the proposed particle swarm optimization, this paper completes the design of a new particle swarm optimization–based multirobot task allocation method. The performance of the new particle swarm optimization–based multirobot task allocation method is tested through three different allocation cases against four well-known evolutionary methods. Experimental results confirm that the proposed method generally outperforms its contenders in terms of the solution quality. Moreover, the proposed method performs slightly better than the majority of its peers as far as the computation time is concerned.

Introduction

Since the cooperation of multiple robots can enhance the overall performance of the multirobot system, the application of teams of multiple robots to accomplish some important unmanned tasks has become increasingly common over the past few decades. ^1,2 A crucial and fundamental issue of the cooperative robotics is the task allocation problem (TAP), which aims to determine robot-to-task assignments, while optimizing a set of predefined performance criteria. ³ In virtue of its importance in the multirobot system, TAP has persistently aroused great research interest in recent years. ^3,4 However, the TAP has naturally been proven to be a non-deterministic polynomial-time hard (NP-hard) problem, and the NP-hard nature of TAP results in difficulties in solving this problem. ^3,4 Therefore, far more effective optimization methods are always demanded in the field of TAP.

Based on the best knowledge of the authors, the market-based allocation approaches can be one of main streams in solving TAP. The most representative market-based allocation approach could be the contract network protocol (CNP) method which was proposed by Smith in 1980. ⁵ Since the initial development of CNP, many researchers have committed themselves to developing some other market-based allocation methods in order to further efficiently solve TAP. A market-based allocation method for the energy-efficient execution problem has been presented by Haghighi in the study by Mo. ⁶ The authors in the study by Gautham et al. ⁷ have developed a distributed task allocation method based on auction and consensus principles for the multirobot system in healthcare facilities. A distributed market-based method where each robot bids for each task has been established in the study by Sahar et al.. ⁸ Thanks to their formidable reliability and scalability, the market-based allocation methods are particularly suitable to find optimal solutions for the distributed multirobot system. ⁴ But, in order to determine optimal robot-to-task assignments, robots in the market-based allocation approach often need to cooperate with each other through explicit communication. Once the commutation network of the multirobot system is interrupted or encounters some other potential issues, the performance of the market-based allocation approach would significantly degrade, which may constraint the implementation of the market-based method on large-scale allocation problems. ⁴

Probably because of their population-based nature and promising parallel search power to find satisfied solutions to some complex NP-hard problems within tractable time, many different evolutionary algorithms (EAs) such as modified differential evolution (mDE), ³ genetic algorithm (GA), ⁹ harmony search algorithm (HSA), ¹⁰ and hybrid ant colony optimization (HACO) ¹¹ have been proposed for solving the TAP within the last decade. Among the currently existing EAs, due to its simplicity and fast convergence speed, particle swarm optimization (PSO) algorithm has been widely used to solve TAP in recent years. A modified binary PSO (mbPSO) task allocation method has been proposed to handle the real-time task allocation in the study by Prescilla et al.. ¹² Since the cognitive and social acceleration parameters are fixed to 2 in mbPSO, this method may have difficulties in adaptively adjust the global and local search abilities of particles. The authors in the literature ¹³ have developed a motivated guaranteed convergence PSO (M-GCPSO) to handle TAP under communication constraints. Because different models of motivations of robots are incorporated into M-GCPSO, this method may be problem-dependent, which may restrict its implementation on some other allocation problems. An intrinsically motivated PSO (IMPSO) method, in which an intrinsic motivation for the fitness function is integrated with PSO, has been presented to solve TAP in the study by Klyne and Merrick. ¹⁴ As values of the inertia weight of particles are either fixed to 0.729 or 0 by the noisy of the fitness function in IMPSO, this method may be unable to well strengthen the exploration capabilities of particles in the early stages of the evolution.

Proverbially, the optimization efficiency of the standard PSO is restricted by its some typical drawbacks such as the inefficient ability in adjusting exploration and exploration as well as the high possibility of trapping into stagnation. ^15,16 In order to enhance of performance of PSO, these two typical flaws of the standard PSO need to be remedied or overcome. For this purpose, from the initial development of the standard PSO, ¹⁷ a considerable amount of research work has been devoted to developing different PSO methods. Aiming at avoiding particles being locked into stagnation, Clerc ¹⁸ has developed a modified PSO, called standard PSO 2011 (SPSO 2011). In order to guarantee the particle can search the solution space with a non-zero velocity, an additional random point is added to the velocity of the particle in SPSO 2011. Yet, since the three control parameters of particles are fixed, SPSO 2011 may still suffer from the difficulties in trade offing exploration and exploitation. ¹⁹

Because the exploration and exploitation capabilities of PSO heavily rely on its three control parameters that are the inertia weight and the cognitive and social acceleration parameters, developing PSO through different parameter settings has been extensively studied in recent years. In an attempt to well balance exploration and exploitation, Akbari and Ziarati ²⁰ have proposed a rank-based PSO (RBPSO), in which the γ best particles are considered to contribute to the updating of the three control parameters of particles. In order to strike a good balance between exploration and exploitation, some other terrific research work focusing on improving PSO through different parameter setting principles can be also found in literatures. ^{21 –23} Reviewing the excellent work done in literatures, ^{20 –23} despite improving the exploration and exploitation abilities of PSO through some advanced parameter setting principles, the convergence of their proposed PSO methods remains unknown or uncertain. Since the three main control parameters of particles not only influence the exploration and exploitation capabilities of PSO but also determine the convergence of PSO, and the convergence property is of great importance in the field of PSO development, the issue of guaranteeing the convergence of PSO must be addressed when improving PSO through different parameter setting principles. ²⁴ Unfortunately, as a stochastic method, like the most of the stochastic methods, the stochastic nature of PSO imposes on difficulties in the analytical investigation of the convergence of PSO. ²⁴

Aiming at remedying the two aforementioned shortcomings of the standard PSO to enhance the performance of PSO, this paper proposes a new enhanced PSO method, called self-adaptive evolutionary game particle swarm optimization (SAEGPSO), which integrates SPSO 2011 with evolutionary game theory (EGT). In order to drag particles in SAEGPSO out of stagnation, SAEGPSO first borrows the moving rules defined in SPSO 2011 to update the velocities and positions of particles. Then, in an attempt to adaptively adjust the exploration and exploitation powers of SAEGPSO, a novel self-adaptive strategy determined by the evolutionary stable strategies of EGT and the iteration number of SAEGPSO is presented to fine-tune the three main control parameters of particles. Moreover, since the three control parameters have profound impacts on the convergence of PSO, this paper also analytically investigates the convergence of SAEGPSO and provides a convergence-guaranteed parameter selection principle for the proposed method.

Leveraging the development of SAEGPSO, this paper completes the design of a SAEGPSO-based task allocation method for the multirobot task allocation problem (MRTAP). In order to reduce the optimization burden and add more diversification to solutions, the feasibility-based rule presented by Kalyanmoy et al. ²⁵ is adopted to handle constraints of the MRTAP in the SAEGPSO-based allocation method. Finally, the developed SAEGPSO-based task allocation method is verified through three different task allocation cases against mDE, ³ GA, ⁹ SPSO 2011, ¹⁸ and improved linearly decreasing weight particle swarm optimization (ILDWPSO). ²⁶ The simulation results indicate that the proposed method performs superior to its competitors as far as the solution quality is considered. Moreover, the computation time of the proposed method is comparable with that of SPSO 2011 and slightly better than those of the other compared methods.

The rest of this paper is organized as follows. The second section states and formulates the studied MRTAP. After simply reviewing SPSO 2011 and EGT, SAEGPSO is presented in the third section. The fourth section mainly analytically investigates the convergence of SAEGPSO. The fifth section describes the SAEGPSO-based allocation method for the MRTAP. Simulations and comparisons are conducted in the sixth section. The seventh section ends this study by drawing conclusions and showing some potential future work.

Formulation of MRTAP

This study concentrates on the multitask single-robot (MT-SR) allocation instance, ²⁷ in which each robot can execute multiple tasks and each task only needs a single robot to execute. Given a set of $\mathit{T}=\{T_1,T_2\mathrm{,...},T_M\}$ tasks and $\mathit{R}=\{R_1,R_2\mathrm{,...},R_N\}$ robots, the goal of the MRTAP is to find a conflict-free task execution sequence for each robot, while optimizing the defined global objective function of the multirobot system. Here, an assignment is said to be free of conflict if each task is assigned no more than once. ¹

In this paper, the global objective function of the multirobot system is the summation of objective functions of each robot. Similar to the study by Shin and Zhang, ²⁸ three objective functions that are the traveling distance, denoted as J_di, the paid cost, denoted as J_ci, and the gained benefit, denoted as J_bi, are considered for each robot i in this article. The traveling distance J_di denotes the total traveling distance that robot i needs to move after the robot accomplishing all assigned tasks. ²⁸ The paid cost J_ci stands for the cost that robot i needs to pay when it fails to execute the assigned tasks. ²⁸ The gained benefit J_bi represents the benefit that robot i gains after the robot executing the assigned tasks. ²⁸

In real-world applications, when each robot can successfully execute some tasks, the multirobot system can obtain a reward. Therefore, in the studied MRTAP, the gained benefit J_bi is designed to simulate the reward that the multirobot system can obtain when the robot successfully completes the assigned tasks. On the other hand, in real-world applications, since the robot may fail to accomplish some tasks, especially in some complex or hazard environments, the multirobot system needs to pay some costs in the case where some tasks are not assigned or executed. As a result, in the studied MRTAP, the paid cost J_ci is designed to simulate the cost that the multirobot system needs to pay when some tasks are not assigned or executed. Obviously, the more tasks the team of robots execute, the more gained benefit the multirobot system may gain; however, the more risks the multiple-robot system has to take to accomplish the assigned tasks, that is, the more paid cost the multirobot system may need to pay. Therefore, there exist conflicts between the paid cost J_ci and the gained benefit J_bi in the MRTAP. ²⁸

Similar to the study by Shin and Zhang, ²⁸ this study also focuses on the numerical value of each objective function. Using a binary decision variable x_ij to indicate whether task j is assigned to robot i or not, the three objective functions of each robot are calculated as follows ²⁸

$\begin{array}{l}{J}_{di}={\displaystyle \sum _{j=1}^M}{x}_{ij}\cdot d_{ij}\hfill \\ J_{ci}=[1-{\displaystyle \prod _{j=1}^M}(1-x_{ij}\cdot V{P}_j)]\cdot A{V}_i\hfill \\ J_{bi}={\displaystyle \sum _{j=1}^M}{x}_{ij}\cdot E{P}_i\cdot A{T}_j\hfill \end{array}$

1

where AV_i denotes the tactical value of robot i. VP_j denotes the violated execute probability of task j to all robots. AT_j is the tactical value of the j th task. EP_i is the execution probability of robot i to all tasks. d_ij represents the Euclidean distance between the current position of robot i and the position of task j. Here, it is important to note that since the mission of the MRTAP is to determine a conflict-free task execution sequence for each robot i in the multirobot system, the current position of robot i is determined by the positions of tasks in the task execute sequence assigned to the robot. For example, the execution sequence of robot i is given as [Task _A , Task _B ], namely, robot i is first assigned to execute Task _A , and then execute Task _B . When robot i is first allocated to complete Task _A , the current position of robot i is the initial position of this robot in the working environment. Next, after robot i completing Task _A , the robot is then assigned to execute Task _B from the position of Task _A , that is, the current position of robot i is now changed into the position of Task _A . Therefore, in the given example, the total traveling distance J_di of robot i equals to the summation of the distance between the initial position of the robot to Task _A and the distance between Task _A and Task _B .

Clearly, J_di and J_ci need to be minimized, while J_bi must be maximized in order to guarantee the multirobot system to gain maximum benefit with minimum traveling distance and paid cost. Therefore, the MRTAP can be formulated into a combinatorial optimization problem using the model described in equation (1) as follows

$\text{min} J={\displaystyle \sum _{i=1}^N}{a}_1\cdot J_{di}+a_2\cdot J_{ci}-a_3\cdot J_{bi}$

2

$\text{Subject}\text{to}:\{\begin{array}{l}{\displaystyle \sum _{i=1}^N{x}_{ij}\le 1,\forall j\in M(3)}\hfill \\ {\displaystyle \sum _{j=1}^N{x}_{ij}\le 1,\forall j\in N(4)}\hfill \\ x_{ij}\in \{0,1\},\forall (i,j)\in N\times M(5)\hfill \end{array}$

where a₁, a₂, and a₃ are three weighting parameters indicating the relative importance of J_di, J_ci, and J_bi ( $a_1=a_2=a_3=1$ in this paper). M and N denote the total number of tasks and robots, respectively. x_ij is a binary decision variable, where $x_{ij}=1$ indicates assigning robot i to task j, otherwise, $x_{ij}=0$ . Due to the limitation of the payload of each robot, the maximum number of tasks that each robot can execute is constrained, that is, the payload constraint of each robot is limited. Here, L_i is used to denote the payload constraint of robot i, which stands for the maximum number of tasks that robot i can complete.

The proposed SAEGPSO method

Since SPSO 2011 and EGT are incorporated into the proposed SAEGPSO method, both these methods are first reviewed in this section. Afterwards, the proposed SAEGPSO is presented at the end of this section.

Review of SPSO 2011

Assume that the position and velocity of the i th particle at iteration t are denoted as $X_i(t)=[x_{i1}(t)\mathrm{,...},x_{in}(t)]$ and $V_i(t)=[v_{i1}(t)\mathrm{,...},v_{in}(t)]$ , respectively. Let ${\text{Pbest}}_i(t)=[{\text{pbest}}_{i1}(t)\mathrm{,...},{\text{pbest}}_{in}(t)]$ represent the personal best solution of particle i at iteration t. Let $\text{Gbest}(t)=[{\text{gbest}}_1(t),\dots ,{\text{gbest}}_n(t)]$ stand for the global best solution of the swarm at iteration t, where n is the dimension of the search space. Suppose $B{C}_i(t)$ denotes the iso-barycentre of particle i at iteration t. The authors in the literature ¹⁸ suggested to calculate the barycenter of ${\text{BC}}_i(t)$ as follows

${\text{BC}}_i(t)=X_i(t)+\frac{c_1[{\text{Pbest}}_i(t)-X_i(t)]+c_2[\text{Gbest}(t)-X_i(t)]}{3}$

6

where c₁ and c₂ are two positive parameters denoting the cognitive and social acceleration parameters, respectively.

In SPSO 2011, a random point, denoted as ${X^{\prime }}_t(t)$ , is first generated in a hypersphere, represented as $\text{HP}({\text{BC}}_i(t),||{\text{BC}}_i(t)-X_i(t)||)$ . Here, the center of the hypersphere $\text{HP}({\text{BC}}_i(t),||{\text{BC}}_i(t)-X_i(t)||)$ is located at ${\text{BC}}_i(t)$ and the radius of $\text{HP}({\text{BC}}_i(t),||{\text{BC}}_i(t)-X_i(t)||)$ is determined by the Euclidean distance between ${\text{BC}}_i(t)$ and the particle’s current position $X_i(t)$ . Then, ${X^{\prime }}_t(t)$ is regarded as a disturbance added to the velocity of the particle in order to guarantee the particle can keep exploring the search space, so that the stagnation issue can be avoided in SPSO 2011. ⁸ After ${X^{\prime }}_t(t)$ is produced in the hypersphere, particles in SPSO 2011 update their velocities and positions based on the moving rules as follows ¹⁸

$\{\begin{array}{c}{V}_i(t+1)=\omega V_i(t)+{X^{\prime }}_i-X_i(t)(7)\\ X_i(t+1)=V_i(t+1)+X_i(t)(8)\end{array}$

where ω is a real coefficient denoting the inertia weight of the particle.

Also, it was suggested in the literature ¹⁸ to fix the three control parameters of particles in SPSO 2011 as follows

$\{\begin{array}{l}\omega =0.7213\hfill \\ c_1=c_2=1.1931\hfill \end{array}$

9

As given in equation (7), because ${X^{\prime }}_i(t)$ is introduced to the velocity of the particle, the particle can keep exploring the search space with a non-null velocity. Therefore, the non-stagnation property can be achieved in SPSO 2011. ¹⁸ Due to this main advantage of SPSO 2011, this method has been widely investigated in literatures. ^{29 –31}

Review of EGT

Suppose a population has a certain amount of players who are playing a game in competition with one another. The purpose of each player in the competition game is to maximum its fitness value by selecting different strategies both based on the experience of the player and those of its competitors. The basic idea of playing this kind of game is that the more fit a strategy is at the current moment, the more likely it will be employed by more players in the future. ³² In order to answer how strategies adopted by players change over time in this kind of game, modeling over different types of animal conflicts, EGT was originally proposed in the context of biology. ³³ The requirement on modeling the evolution phenomena leads to the application of mathematical theory of games to interpret aspects of the evolutionary process. ³⁴

Although EGT was originally developed in the context of the biologic science, it has increasingly aroused research interest of the other researchers from different fields such as the economists, sociologists, and the philosophers over the last few decades. Such an interpretation to the wide application of EGT on different fields may be explained by the following facts. Firstly, the notion of the evolution needs to be understood as the change of beliefs and norms over time. ³⁴ Secondly, the modeling of the variations of strategies offers a social aspect which exactly matches the social system interactions. ³⁴ Lastly, EGT can dynamically model the social interactions of different players within a population, which is one of missing elements of the classical game theory. ³⁴ In EGT, the final output of a game is determined by the equilibrium point, called the evolutionary stable strategy (ESS). The principle of EGT is both based on the performances gained by a single strategy and the presence of other strategies. In addition to the ESS, the replicator dynamic equation (RDE) is another key concept in EGT, which is used to describe the adaptation of the population over time. Considering their importance in EGT, both ESS and RDE are described in the herein subsections.

Description of ESS

In EGT, a strategy can be considered as an ESS such that, when all players in a population play this strategy, no mutant strategy can invade the population under the influence of the natural selection. ³² Suppose a population has some players following the strategy p₁ and a few players adopting strategy p₂. Assume $f(p_1,\overline{S_d})$ represents the payoff, that is, the fitness value of a specific player using strategy p₁. Then, for a minimization optimization problem (i.e., players in EGT tend to minimize their fitness values), p₁ is said to be an ESS if $f(p_2,\overline{S_d})>f(p_1,\overline{S_d})$ for all $p_2\ne p_1$ and any sufficiently small μ ³² . Here, $\overline{S_d}$ denotes the strategy distribution of the population and can be calculated as follows ³²

$\overline{S_d}=(1-\mu )p_1+\mu p_2$

10

where μ denotes the frequency of p₂ in the population.

Description of RDE

In order to discuss the stability of EGT, it is essential to define a dynamic for the EGT ³² to describe the impacts of the variations of strategies on the population. ³² For this purpose, Taylor and Jonker ³² have developed the replicator dynamic equation to model the adaptation of strategies in a population. From the work done by Taylor and Jonker, ³² it can be found that the matrix game is a most commonly used way to depict the strategy interactions among players, which can be described using the following notations: (1) e_i represents the i th pure strategy, where $i=1,2,\mathrm{...},m$ and m is the number of pure strategies; (2) $K=f(e_i,e_j)$ denotes the payoff matrix with the size of $m\times m$ ; (3) ${\nabla }^m\equiv \{p=(p_1\mathrm{,...},p_m)|{\displaystyle {\sum }_{i=1}^m}{p}_i=1,0\le p_i\le 1\}$ is the set of mixed strategies (probability distribution over the pure strategy e_i); (4) $f(p,q)=pK{q}^T$ stands for the payoff, that is, the fitness value of players using strategy p facing other players adopting strategy q. In EGT, RDE is used to describe the difference between the payoff of the i th strategy and the mean payoff of the population, which is defined as follows ³²

${\dot{p}}_i=-p_i(e_i\cdot K{p}^T-p\cdot K{p}^T)$

11

The RED defined by equation (11) is one of the most widely accepted evolutionary dynamics used to describe the evolution of the strategy frequency p_i. ³⁴ Thanks to its promising search stability and population-based nature, EGT has been widely studied in literatures. ^{24,34 –37} For more details about EGT, the reader is also suggested to refer to these excellent works done in literatures. ^{24,34 –37}

Modeling of SAEGPSO

Universally, striking a good balance between the exploration and exploitation capabilities of PSO is of great importance in the design of a PSO-based optimizer in order to enhance the performance of PSO. Ideally, on one hand, the exploration ability of PSO must be strengthened in the early stages of the evolution in order to encourage particles to search through the whole search space, rather than converging toward the current population-best solution. ^15,20,22,38 On the other hand, the exploitation capability of PSO needs to be promoted in the latter phases of the evolution, so that particles are more likely to be encouraged to search carefully in a local region to increase the likelihood of finding the global optimal solution. ^15,20,22,38 As is known to all, such two abilities of PSO are dramatically affected by its three control parameters. The basic principles regarding how the three control parameters of PSO influence such two capabilities can be summarized as follows: (1) the exploration power of PSO benefits from a large inertia weight, while the exploitation power benefits from a small inertia weight ^15,20,22,38 ; (2) comparing with the social component, a large cognitive acceleration parameter leads particles to global searches and, consequently, promotes the exploration ability of PSO ^15,38,39 ; (3) comparing with the cognitive component, a large social acceleration parameter turns particles into local searches and thus facilitates the exploitation ability of PSO. ^15,38,39

Based on the basic principles stated above, despite achieving the non-stagnation property in SPSO 2011, this method may still suffer from the potential difficulties in well trade offing the exploration and exploitation abilities, since the three control parameters of particles in SPSO 2011 are fixed and there does not exit any distinction between the cognitive and social acceleration parameters. This paper concentrates on improving PSO through remedying the two aforementioned shortcomings of the standard PSO and proposes an improved method, called SAEGPSO, which combines SPSO 2011 with EGT. The motivation of integrating SPSO 2011 with EGT is that since both these two approaches are population-based methods, combining them together can naturally leverage their advantages to improve the performance of the enhanced method. Motivated by this consideration, in order to avoid particles falling into stagnation, the moving rules defined by equations (6) to (8) in SPSO 2011 are applied to guide the movements of particles in SAEGPSO. In order to well trade off the exploration and exploitation powers of SAEGPSO, a newly developed self-adaptive strategy determined by the evolutionary stable strategies of EGT and the iteration number of SAEGSPO is proposed to fine-tune the three control parameters of particles in SAEGPSO.

Prior to presenting how to combine the EGT principle with the iteration number of SAEGPSO in the developed self-adaptive strategy, there is necessity to recall some important notations such as “player”, “strategy,” and “payoff matrix” of EGT, since these notations will be analogized in the newly developed self-adaptive strategy presented in SAEGPSO. Players in EGT denote candidate solutions to an optimization problem, which push the evolution process of EGT. Strategies in EGT denote the candidate principles that each player can select to maximize its playoff, that is, the fitness value. The payoff matrix is a matrix which is constructed by the payoffs, that is, the fitness values of all players in the population. After recalling these notations in EGT, before the introduction of the proposed self-adaptive strategy, we first present the analogy between SAEGPSO and EGT as follows: (1) each player in EGT analogizes each particle in SAEGPSO; (2) each particle in SAEGPSO has three candidate strategies, which are moving only based on its inertia weight, moving only based on its personal best solution, and moving only based on the global best solution of the swarm, respectively; (3) the average of performance gained by the particle in SAEGPSO playing a specific strategy among the three candidate strategies consists of the payoff matrix of EGT.

Let e₁, e₂, and e₃, respectively, denote that the particle follows the strategy only based on its inertia weight, the strategy only based on its personal best solution, and the strategy only based on global best solution of the swarm. Then, the payoff matrix used in this article is defined as follows

$K=\left\{\begin{array}{ccc}f(e_1)& \frac{f(e_1)-f(e_2)}{2}& \frac{f(e_1)-f(e_3)}{2}\\ \frac{f(e_2)-f(e_1)}{2}& f(e_2)& \frac{f(e_2)-f(e_3)}{2}\\ \frac{f(e_3)-f(e_1)}{2}& \frac{f(e_3)-f(e_2)}{2}& f(e_3)\end{array}\right\}$

12

where $f(e_i)$ ( $i=1,2,3$ ) is the payoff, that is, the fitness value that the particle obtains by only using the strategy e_i ( $i=1,2,3$ ).

At each iteration t, the ESS of EGT is used to weight the mean of the fitness obtained by each particle, which denotes the ratio of each strategy when the swarm of SAEGPSO reaches a stable equilibrium point. At each iteration t, the ESS applied in this study is given as follows

$\text{ESS}(t)=[Z_1(t),Z_2(t),Z_3(t)]$

13

subject to (s.t.)

$(\forall t\in \mathbb{N})({\displaystyle \sum _{i=1}^3}{Z}_i(t)=1)$

14

The value of $f(e_i)$ ( $i=1,2,3$ ) is gained from the previous experience of the particle and calculated as follows

$f(e_i)=\frac{1}{t}{\displaystyle \sum _{t_1=1}^t}{Z}_i(t_1-1)F(X(t_1-1))$

15

where t ∈ ℕ denotes the iteration number of SAEGPSO. $F(X(t_1-1))$ is the previously obtained fitness value of the particle X at iteration $(t_1-1)$ , where $t_1\in \mathbb{N}$ and $1\le t_1\le t$ . $Z_i(t_1-1)$ is the ratio of the i th strategy at iteration ( $t_1-1$ ). Note that in order to guarantee each particle in SAEGPSO to equally follow the three strategies at the initial step of the evolution of SAEGPSO, we set that $\text{ESS}(0)=[1/3,1/3,1/3]$ .

Once the payoff $f(e_i)$ ( $i=1,2,3$ ) of each particle is obtained using equation (15), it is then used to fill the payoff matrix defined by equation (12). Once the payoff matrix given by equation (12) is filled, the RDE given by equation (11) is implemented to obtain the associated ESS, that is, $\text{ESS}(t)=[Z_1(t),Z_2(t),Z_3(t)]$ at each iteration t. After $\text{ESS}(t)=[Z_1(t),Z_2(t),Z_3(t)]$ is obtained, the three ratios $Z_1(t)$ , $Z_2(t)$ , and $Z_3(t)$ are then applied to adjust the three control parameters of each particle in the developed self-adaptive strategy in SAEGPSO as follows

$\omega (t+1)=({\omega }_s-{\omega }_f)\text{exp}(-\frac{{\delta }_{\omega }t}{\beta })+{\omega }_f$

16

$c_1(t+1)=(c_{1s}-c_{1f})\text{exp}(-\frac{{\delta }_{c_1}t}{\beta })+c_{1f}$

17

$c_2(t+1)=(c_{2s}-c_{2f})\text{exp}\left(\frac{{\delta }_{c_2}t}{\beta }\right)+c_{2f}$

18

where

${\delta }_{\omega }=\frac{{\omega }_s-{\omega }_f}{t_{\text{max}}}$

19

${\delta }_{c_1}=\frac{c_{1s}-c_{1f}}{t_{\text{max}}}$

20

${\delta }_{c_2}=\frac{c_{2s}-c_{2f}}{t_{\text{max}}}$

21

$\beta =\frac{Z_1(t)+Z_2(t)}{Z_3(t)+\Delta }$

22

where the subscripts “s” and “f” in ω, c₁, and c₂ denote their initial and final values, respectively. t_max represents the maximum iteration number. Δ is a sufficiently small positive parameter and empirically set to be $\Delta =1e-10$ in this paper in order to avoid the denominator of β given in equation (22) equals to zero. It is important to note that ${\omega }_s>{\omega }_f$ , $c_{1s}>c_{1f}$ , and $c_{2f}>c_{2s}$ in the above self-adaptive strategy in SAEGPSO.

The implementation of the three ratios $Z_1(t)$ , $Z_2(t)$ , and $Z_3(t)$ are used to determine the control parameters of particles in SAEGPSO can be interpreted by the fact that these three ratios denote a stable search direction that particles need to follow. Since the three ratios are obtained based on the previous experience of particles, when they are implemented to determine the control parameters of particles, particles may adapt the “shape” of the search space to optimize the search direction of the swarm. In addition, the stability characteristic of ESS in EGT may be applied to face the potential irregularity of the search space to help particles to avoid clustering into some local optimal, without damaging the search capabilities of particles. As authors in the literature ²² have discovered that updating the control parameters of PSO based on the exponential or double exponential manner can accelerate the convergence of PSO, the three control parameters of particles in SAEGPSO are adaptively adjusted using the exponential function, as defined by equations (16) to (18) in the developed self-adaptive strategy.

Parametric analysis for SAEGPSO

In order to intuitively interpret the rationalities of the developed self-adaptive strategy in SAEGPSO, this subsection conducts a parametric analysis for SAEGPSO. It is clear from equations (16) to (18) that the inertia weight ω and the cognitive acceleration parameter c₁ decrease, while the social acceleration parameter c₂ increases as the iteration number t increases. Thus, based on the basic principles noted above, the exploration ability of SAEGPSO is more likely to be maintained in the early phases of the evolutionary process, which would be reduced over time, so that the exploitation capability may take over the exploration ability quickly in the latter states of the evolution. It is noticeable that, adopting the updating rule of a constant β, trade-offs between these two abilities only vary with respect to the iteration number t.

In addition to the iteration number t, the balance of the search behavior in SAEGPSO is also adapted based on the parameter β. From equations (16) to (18), one can easily observe that ω and c₁ decrease, while the variation in c₂ becomes larger as β increases. According to the basic principles stated above, this implies that the exploration tendency of SAEGPSO will be more preserved for a large β. From equation (22), a large β indicates a relatively large value of ( $Z_1(t)+Z_2(t)$ ). For a relatively large ( $Z_1(t)+Z_2(t)$ ), the search directions of particles could be more stable when particles follow the strategies of their inertia weight and personal best experience, which implies that the population may benefit more through particles following these two strategies. Therefore, according to equations (16) and (17), it is natural to increase ω and c₁ to promote the exploration power of SAEGPSO for a large β. Contrarily, a small β indicates that the value of $Z_3(t)$ is relatively large based on equation (22) (ignoring the sufficiently small constant Δ), which implies that the search directions of particles may be more stable and the swarm could benefit more when particles play the strategy of following the global best experience of the swarm. Hence, as shown in equation (18), for a small β, the variation of c₂ increases in order to enhance the exploitation ability of SAEGPSO.

To summarize, through implementing the self-adaptive strategy defined by equations (16) to (22), the three main control parameters of particles in SAEGPSO are adaptively adjusted, complying with the basic principles of PSO development. Thus, particles may be promoted to search for high-quality solutions and the performance of SAEGPSO can be enhanced. The variation trends of these three parameters of particles in SAEGPSO under different values of β are illustrated in Figure 1.

Figure 1.

Parameter changes with respect to different β in SAEGPSO. (a) Changes of inertia weight ω under different β. (b) Changes of cognitive acceleration parameter c₁ under different β. (c) Changes of social acceleration parameter c₂ under different β. SAEGPSO: self-adaptive evolutionary game particle swarm optimization.

Theoretical investigations of SAEGPSO

As stated in the third section, in order to enhance the performance of PSO via remedying the two typical flaws of the standard PSO, this paper proposes an enhanced PSO, called SAEGPSO. Aiming at dragging particles in SAEGPSO out of stagnation, the updating rules of SPSO 2011 defined by equations (7) and (8) are borrowed by SAEGPSO to update the velocities and positions of particles. On the other hand, focusing on well balancing the exploration and exploitation capabilities of SAEGPSO, a newly developed self-adaptive strategy determined by EGT and the iteration number of SAEGPSO is presented to tune the three control parameters of particles in SAEGPSO. Since one main modification of SAEGPSO lies on proposing a parameter setting principle to update the three control parameters of particles and the three control parameters influence the convergence of PSO, there is necessity to analytically investigate the convergence of SAEGPSO with respect to different values of the three control parameters. Therefore, we carry out a theoretical investigation, including the convergence analysis, the equilibrium point, and the convergence-guaranteed parameter selection principle, pertaining to SAEGPSO in following contents of this section.

Convergence analysis for SAEGPSO

The convergence analysis of PSO is used to determine the boundaries of the three control parameters which guarantee the convergence of PSO. In this subsection, we focus on investigating the convergence of SAEGPSO based on the deterministic model convergence analysis, ³⁰ that is, in the case where ${X^{\prime }}_i={\text{BC}}_i(t)$ . For simplicity and without loss of generality, omitting the subscript i of each variable in equations (7) and (8), the moving rules of particles in SAEGPSO can be rewritten into the matrix form as follows

$\left[\begin{array}{c}X(t+1)\\ V(t+1)\end{array}\right]=\left[\begin{array}{cc}1-c& \omega \\ -c& \omega \end{array}\right]\left[\begin{array}{c}X(t)\\ V(t)\end{array}\right]+\left[\begin{array}{c}c\\ c\end{array}\right]P$

23

where

$c=\frac{c_1+c_2}{3}$

24

$P=\frac{c_1\cdot \text{Pbest}+c_2\cdot \text{Gbest}}{c_1+c_2}$

25

The characteristic equation of system (23) is easily obtained as follows

${\eta }^2-(1+\omega -c)\eta +\omega =0$

26

Assume two roots to equation (26) are denoted as η₁ and η₂. Then, it is easily obtained that

${\eta }_{1,2}=\frac{1+\omega -c\pm \sqrt{{(1+\omega -c)}^2-4\omega }}{2}$

27

Based on the dynamic system theory, the system (23) converges iff the magnitudes of the two characteristic roots are less than 1 (“iff” means “if and only if” in this paper). ⁴⁰ In another word, system (23) converges, iff

$\text{Max}\{\left|{\eta }_1\right|,\left|{\eta }_2\right|\}<1$

28

From equation (27), since η₁ and η₂ may be two real or complex roots, both these two cases are discussed separately so as to easily study the convergence of system (23).

Case 1. η₁ and η₂ are two complex roots, presented as ${\eta }_{1,2}\in ℂ$ , where ℂ denotes the imaginary domain.

Lemma 1

For system (23), ${\eta }_{1,2}\in ℂ$ , iff

$\{\begin{array}{l}1+\omega -2\sqrt{\omega }<c<1+\omega +2\sqrt{\omega }\hfill \\ \omega \ge 0\hfill \end{array}$

29

Proof

It is trivial from equation (26) that ${\eta }_{1,2}\in ℂ$ , iff

${(1+\omega -c)}^2-4\omega <0$

30

The proof of lemma 1 is completed by solving equation (30) with classical mathematical method.

Then, we need to find the conditions on ω and c which guarantee the convergence of system (23) in the case where ${\eta }_{1,2}\in ℂ$ . Here, please recall that the convergence of system (23) is guaranteed iff $\text{Max}\{\left|{\eta }_1\right|,\left|{\eta }_2\right|\}<1$ .

Lemma 2

For ${\eta }_{1,2}\in ℂ$ , system (23) converges, iff

$\{\begin{array}{l}1+\omega -2\sqrt{\omega }<c<1+\omega +2\sqrt{\omega }\hfill \\ 0\le \omega <1\hfill \end{array}$

31

Proof

Notice that the magnitude of any imaginary number H can be calculated as $\left|H\right|=\sqrt{H_r^2+H_c^2}$ , where H_r and H_c are the real and imaginary parts of H, respectively. Therefore, in the case where ${\eta }_{1,2}\in ℂ$ , it can be easily obtained that

$\begin{array}{l}\text{Max}\{\left|{\eta }_1\right|,\left|{\eta }_2\right|\}=\left|{\eta }_1\right|=\left|{\eta }_2\right|=\sqrt{\omega }\hfill \end{array}$

32

Therefore

$\text{Max}\{\left|{\eta }_1\right|,\left|{\eta }_2\right|\}<1\iff \sqrt{\omega }<1$

33

In the case where ${\eta }_{1,2}\in ℂ$ , based on lemma 1, equation (29) must hold. Hence, both taking two conditions that ${\eta }_{1,2}\in ℂ$ and $\text{Max}\{\left|{\eta }_1\right|,\left|{\eta }_2\right|\}<1$ into account, system (23), that is, SAEGPSO, converges in the case where ${\eta }_{1,2}\in ℂ$ , iff

$\{\begin{array}{l}1+\omega -2\sqrt{\omega }\le c\le 1+\omega +2\sqrt{\omega }\hfill \\ 0\le \omega <1\hfill \end{array}$

34

This completes the proof of lemma 2.

Figure 2 displays the three-dimensional convergence domain of SAEGPSO and the corresponding projections of the convergence domain on different parameter planes in the case where ${\eta }_{1,2}\in ℂ$ . For any ${\eta }_{1,2}\in ℂ$ , SAEGPSO converges only if the selections of ω and c locate in the convergence region demonstrated in Figure 2(b).

Figure 2.

The convergence domain of SAEGPSO for ${\eta }_{1,2}\in ℂ$ . (a) 3-D presentation of the convergence domain. (b) Projection of convergence domain on plane (ω, c). (c) Projection of convergence domain on plane (ω, Max $\left\{\left|{\eta }_1\right|,\left|{\eta }_2\right|\right\}$ ). (d) Projection of convergence domain on plane (c, Max $\left\{\left|{\eta }_1\right|,\left|{\eta }_2\right|\right\}$ ). SAEGPSO: self-adaptive evolutionary game particle swarm optimization.

1. Case 2. η₁ and η₂ are two real roots, denoted as ${\eta }_{1,2}\in ℝ$ , where ℝ stands for real-valued domain.

Lemma 3

For equation (23), ${\eta }_{1,2}\in ℝ$ , iff

$\{\begin{array}{l}c\in ℝ,\omega <0\hfill \\ c\le 1+\omega -2\sqrt{\omega }\text{or}c\ge 1+\omega +2\sqrt{\omega },\omega \ge 0\hfill \end{array}$

35

Proof

It is clear from equation (26) that ${\eta }_{1,2}\in ℝ$ , iff

${(1+\omega -c)}^2-4\omega \ge 0$

36

Solving equation (36) using classical mathematical method, the proof of lemma 3 is completed.

Then, let us discover conditions on ω and c which can guarantee the convergence of system (23) under the situation where ${\eta }_{1,2}\in ℝ$ .

Lemma 4

System (23) converges in the case where ${\eta }_{1,2}\in ℝ$ , iff

$\{\begin{array}{l}0<c<2\omega +2,-1<\omega <0\hfill \\ 0<c\le 1+\omega -2\sqrt{\omega }\text{or}1+\omega +2\sqrt{\omega }\le c<2\omega +2,0\le \omega <1\hfill \end{array}$

37

Proof

It is trivial from equations (27) and (28) that $\text{Max}\{\left|{\eta }_1\right|,\left|{\eta }_2\right|\}<1$ holds for ${\eta }_{1,2}\in ℝ$ , iff

$-1<\frac{1+\omega -c\pm \sqrt{{(1+\omega -c)}^2-4\omega }}{2}<1$

38

Namely, iff

$c-\omega -3<\pm \sqrt{{(1+\omega -c)}^2-4\omega }<c-\omega +1$

39

As ${\eta }_{1,2}\in ℝ$ , it is trivial from equation (39) that

$(38)\iff \{\begin{array}{l}c-\omega -3<-\sqrt{{(1+\omega -c)}^2-4\omega }\hfill \\ \sqrt{{(1+\omega -c)}^2-4\omega }<c-\omega +1\hfill \end{array}$

40

Simplifying the right-hand inequalities in equation (40) using classical mathematical method, generates

$(38)\iff \{\begin{array}{l}2\omega +2-c>0\hfill \\ c>0\hfill \end{array}$

41

Based on lemma 3, equation (35) must hold when ${\eta }_{1,2}\in ℝ$ . Thereafter, both considering two conditions that ${\eta }_{1,2}\in ℝ$ and $\text{Max}\{\left|{\eta }_1\right|,\left|{\eta }_2\right|\}<1$ , system (23), that is, SAEGPSO, converges in the case where ${\eta }_{1,2}\in ℝ$ , iff

$\{\begin{array}{l}0<c<2\omega +2,-1<\omega <0\hfill \\ 0<c\le 1+\omega -2\sqrt{\omega }\text{or}1+\omega +2\sqrt{\omega }\le \phi <2\omega +2,0\le \omega <1\hfill \end{array}$

42

For ${\eta }_{1,2}\in ℝ$ , the three-dimensional convergence domain of SAEGPSO and the corresponding projections of the convergence domain on different parameter planes are visualized in Figure 3. Only if any parameter choice of ω and c locates within the region illustrated in Figure 3(b), the convergence of SAEGPSO can be guaranteed in the case where ${\eta }_{1,2}\in ℝ$ .

Figure 3.

The convergence domain of SAEGPSO in the case where ${\eta }_{1,2}\in ℝ$ . (a) 3-D presentation of the convergence domain. (b) Projection of convergence domain on plane (ω, c). (c) Projection of convergence domain on plane (ω, Max $\left\{\left|{\eta }_1\right|,\left|{\eta }_2\right|\right\}$ ). (d) Projection of convergence domain on plane (c, Max $\left\{\left|{\eta }_1\right|,\left|{\eta }_2\right|\right\}$ ). SAEGPSO: self-adaptive evolutionary game particle swarm optimization.

Then, let us combine the conclusions drawn in lemma 3 and lemma 4 together, we can conclude that the system (23), that is, SAEGPSO, converges, iff

$\{\begin{array}{l}0<c<2\omega +2\hfill \\ -1<\omega <1\hfill \end{array}$

43

Figure 4 shows the three-dimensional convergence domain of SAEGPSO and the corresponding projections of the convergence domain on different parameter planes. The convergence of SAEGPSO can be guaranteed as long as any parameter choice of ω and c locates within the triangle region shown in Figure 4(b).

Figure 4.

The convergence domain of SAEGPSO. (a) 3-D presentation of the convergence domain. (b) Projection of convergence domain on plane (ω, c). (c) Projection of convergence domain on plane (ω, Max $\left\{\left|{\eta }_1\right|,\left|{\eta }_2\right|\right\}$ ). (d) Projection of convergence domain on plane (c, Max $\left\{\left|{\eta }_1\right|,\left|{\eta }_2\right|\right\}$ ). SAEGPSO: self-adaptive evolutionary game particle swarm optimization.

The equilibrium point of SAEGPSO

The equilibrium point of PSO is a stable point that particles converge toward at the end of the evolution. After studying the convergence of SAEGPSO in ‘Convergence analysis for SAEGPSO’ section, we still need to find the equilibrium point of this method. Calculating limits on both sides of equation (23) yields

$\{\begin{array}{l}\underset{t\to \infty }{{\displaystyle \text{lim}}}X(t+1)=\omega \underset{t\to \infty }{{\displaystyle \text{lim}}}V(t)+c\underset{t\to \infty }{{\displaystyle \text{lim}}}(P-X(t))\hfill \\ \underset{t\to \infty }{{\displaystyle \text{lim}}}V(t+1)=\underset{t\to \infty }{{\displaystyle \text{lim}}}X(t)+\underset{t\to \infty }{{\displaystyle \text{lim}}}V(t)\hfill \end{array}$

44

Obviously, $\underset{t\to \infty }{{\displaystyle \text{lim}}}X(t+1)=\underset{t\to \infty }{{\displaystyle \text{lim}}}X(t)$ and $\underset{t\to \infty }{{\displaystyle \text{lim}}}V(t+1)=\underset{t\to \infty }{{\displaystyle \text{lim}}}V(t)$ when the convergence of SAEGPSO is guaranteed. Hence, substituting these two equations into equation (44), produces

$\{\begin{array}{c}\underset{t\to \infty }{\text{lim}}X(t)=P=\frac{c_1\cdot \text{Pbest}+c_2\cdot \text{Gbest}}{c_1+c_2}(45)\\ \underset{t\to \infty }{\text{lim}}V(t)=0(46)\end{array}$

where c₁ and c₂, respectively, denote the cognitive and social acceleration parameters of the particle. $\text{Pbest}$ and $\text{Gbest}$ represent the personal best solution of the particle and global best solution of the swarm, respectively.

Convergence-guaranteed parameter selection principle for SAEGPSO

In ‘Convergence analysis for SAEGPSO’ section, it has been proven that SAEGPSO converges iff the condition given by equation (43) is satisfied. Here, the remaining task is to answer how to set the initial and final values of the three control parameters of particles in the newly developed self-adaptive strategy to satisfy the condition given by equation (43), so that the convergence of SAEGPSO can be sufficiently guaranteed. Thus, this subsection provides a convergence-guaranteed parameter selection principle for SAEGPSO.

Lemma 5

SAEGPSO converges, only if the initial and final values of the three control parameters meet the following conditions

$\{\begin{array}{l}0<c_{1s}+c_{1f}<6{\omega }_f+6\hfill \\ -1<{\omega }_f<{\omega }_s<1\hfill \\ c_{1s}=c_{2f}>c_{1f}=c_{2s}>0\hfill \end{array}$

47

Proof

Substituting $c=(c_1+c_2)/3$ into the condition given by equation (43), we can obtain that SAEGPSO converges, iff

$\{\begin{array}{l}0<c_1+c_2<6\omega +6\hfill \\ -1<\omega <1\hfill \end{array}$

48

If $c_{1s}=c_{2f}>c_{1f}=c_{2s}>0$ , it appears from equations (17), (18), (20), and (21) that $c_1+c_2=c_{1s}+c_{1f}>0$ can be trivially satisfied for any particle at any iteration in SAEGPSO. Besides, we can easily make an observation from equations (17) to (19) that ${\omega }_f\le \omega \le {\omega }_s$ , $c_{1f}\le c_1\le c_{1s}$ and $c_{2s}\le c_2\le c_{2f}$ for any particle at any iteration in SAEGPSO. Thus, it is clear that

$\{\begin{array}{l}0<c_{1s}+c_{1f}<6{\omega }_f+6\hfill \\ -1<{\omega }_f<{\omega }_s<1\hfill \\ c_{1s}=c_{2f}>c_{1f}=c_{2s}>0\hfill \end{array}\Rightarrow \{\begin{array}{l}0<c_1+c_2<6\omega +6\hfill \\ -1<\omega <1\hfill \end{array}$

49

Because the right-hand side of equation (49) is the necessary and sufficient condition for the convergence of SAEGPSO, the proof of lemma 5 is completed based on the relationship given by equation (49).

It is important to note that values of ω_s, ω_f, $c_{1s}$ , $c_{1f}$ , $c_{2s}$ , and $c_{2f}$ are given beforehand. Therefore, the convergence condition given by equation (47) is easy to be met through setting properly values of these predefined parameters. The convergent position and velocity trajectories of the particle in SAEGPSO are demonstrated in Figure 5, where values of these parameters are empirically suggested to be as follows: ${\omega }_s=0.9$ , ${\omega }_f=0.1$ , $c_{1s}=c_{2f}=2.5$ , and $c_{1f}=c_{2s}=0.1$ .

Figure 5.

Convergent trajectory of the particle in SAEGPSO. (a) Convergent position trajectory. (b) Convergent velocity trajectory. SAEGPSO: self-adaptive evolutionary game particle swarm optimization.

The SAEGPSO-based allocation method for the MRTAP

When using the proposed SAEGPSO method to develop a SAEGPSO-based allocation method for the MRTAP, two issues need to be first addressed. The first one is to how to encode particles, that is, how to build a mapping between positions of particles and solutions to the MRTAP. The second issue is to answer how to handle constraints of the MRTAP in the SAEGPSO-based allocation method. The main mission of this section is to answer these two questions.

Encoding particles

As stated in the second section, the purpose of the studied MRTAP is to find a task execution sequence for each robot. To achieve this goal, as illustrated in Figure 6, the position of the particle in SAEGPSO is encoded as a $1\times M$ vector, where M is the total number of tasks. As shown in Figure 6, let $[x_j]$ and $\left\{x_j\right\}$ denote the integer and decimal parts of the element x_j ( $j=1,2,\mathrm{...},M$ ) in the position vector of the particle, respectively. Here, $[x_j]=i$ indicates assigning the j th task to the i th robot. Besides, the magnitude of the decimal part $\left\{x_j\right\}$ of x_j represents the sequence that the j th task is executed by the i th robot. The smaller the decimal part of the element x_j has, the more prior the corresponding task is performed by the corresponding assigned robot. For the position vector of any considered particle, if we sort all elements having the same integer part in ascending order according to the decimal parts of those elements, then we can easily obtain the task execution order for each robot.

Figure 6.

Representation of the particle in SAEGPSO. SAEGPSO: self-adaptive evolutionary game particle swarm optimization.

To concretely explain how to encode and decode particles based on the way described above, an illustrative MRTAP example in which three robots are assigned to execute seven tasks is given here. Assume the position vector of a particle is randomly given as [3.5211, 1.312, 2.23416, 1.45016, 1.01001, 2.36143, 3.25]. According to the way of decoding particles depicted above, the task execution sequences for robots 1, 2, and 3 are: $\text{Task}5\to \text{Task}2\to \text{Task}4$ , $\text{Task}3\to \text{Task}6$ , and $\text{Task}7\to \text{Task}1$ , respectively. The first advantage of representing particles using the way described above is that the storage of the optimizer can be decreased since each particle is encoded as a $1\times M$ vector. Moreover, $[x_j]=i$ formally indicates assigning the i th robot to conduct the j th task, which implies that the binary decision variable $x_{ij}=1$ , and $x_{hj}=0$ for $\forall h\ne i$ ( $h=1,2,\mathrm{...},N$ ) in the MRTAP, where N is the total number of robots. Thus, the second advantage of representing particles applying the way presented above is that the constraints given by equations (3) and (5) of the MRTAP can be inherently satisfied.

Handling constraints of the MRTAP

Because the MRTAP studied in this paper is a constrained optimization problem (COP), the task regarding how to handle constraints of the problem must be addressed in order to efficiently solve this problem. Among the currently existing constraint handling techniques, the penalty function approach ⁴¹ could be one of the most popular approaches. The penalty method can transform a COP into an unstrained optimization problem through introducing several penalty terms to the objective function of the COP. However, properly setting the penalty factors of the penalty function method is difficult and problem-dependent, which is the main drawback of this method. ⁴²

In our SAEGPSO-based task allocation method, we adopt the feasibility-based rule ²⁵ to tackle constraints of the MRTAP in order to release the burden of the optimizer, as well as increase diversification of solutions. As explained in ‘Encoding particles’ section, the constraints given by equations (3) and (5) of the MRTAP are trivially met using the way presented in this paper to encode particles. Therefore, only the constraint given by equation (4) of the MRTAP needs to be tackled. In order to handle the constraint given by equation (4), the constraint violation degree of each particle is first calculated as follows

$\text{viol}={\displaystyle \sum _{i=1}^N}(0,\text{max}(\left|T{S}_i\right|-L_i))$

50

where N is the total number of robots. $\left|{\text{TS}}_i\right|$ is the length of the task execution set of robot i. L_i is the payload constraint of robot i. Note that the information of the task execution set ${\text{TS}}_i$ of each robot can be obtained based on the way of encoding particles as explained in ‘Encoding particles’ section. Again, take the illustrative example as exampled in ‘Encoding particles’ section as an example to introduce how to calculate the constraint violation degrees of particles based on equation (50). Since the task execution sequences of robots 1, 2, and 3 are $\text{Task}5\to \text{Task}2\to \text{Task}4$ , $\text{Task}3\to \text{Task}6$ , and $\text{Task}7\to \text{Task}1$ , respectively, we have that $\left|T{S}_1\right|=3$ , $\left|T{S}_2\right|=\left|T{S}_3\right|=2$ . Assume the payload constraints of the three robots are 1, 2, and 4, respectively. Then, it can be easily obtained that $\text{viol}=2$ for the particle as exampled in ‘Encoding particles’ section.

After calculating constraint violation degrees of particles based on equation (50), the feasibility-based rule presented by Kalyanmoy et al. ²⁵ is then used to choose the dominated solution between any two candidate solutions in the SAEGPSO-based allocation method. For a minimization optimization problem, the feasibility-based rule is depicted as follows: (1) for any two candidate solutions having the same constraint violation degree, the solution with smaller fitness value dominates the solution with larger fitness value; (2) for any two candidate solutions having different constraint violation degrees, the solution with smaller constraint violation degree is preferred over the solution with larger constraint violation degree.

Because the objective function and the constraint violation degree of each particle are considered and compared separately in the feasibility-based rule, no additional penalty factor is required when implementing this rule to handle constraints. Therefore, this rule can reduce the burden of the optimizer. ²⁵ In addition, although the non-feasible solutions violate some constraints, they may also contain some useful information for finding high-quality solutions. ⁴² Consequently, when these non-feasible solutions are considered in the feasibility-based rule, the diversification of solutions and the likelihood of finding better solutions may be increased.

In the proposed SAEGPSO-based task allocation method, in order to assure each task can assign a robot, $[x_j]$ must satisfy the boundary constraint $1\le [x_j]\le N$ , where N is the total number of robots. For any integer part $[x_j]$ of the element x_j in the position vector of the particle violates the boundary constraint, the saturation strategy presented by Khalili-Damghani et al. ⁴¹ is used to modify x_j as follows

$x_j=\{\begin{array}{l}1+\left\{x_j\right\}\text{if}[x_j]<1\hfill \\ x_j\text{if}1\le [x_j]\le N\hfill \\ N+\left\{x_j\right\}\text{otherwise}\hfill \end{array}$

51

where $[x_j]$ and $\left\{x_j\right\}$ are the integer and decimal parts of the element x_j in the position vector of the particle, respectively. For any element x_j in the position vector of the particle, the aim of the saturation strategy given by equation (51) is to only modify the integer part of x_j, so that the integer part of x_j satisfy $1\le [x_j]\le N$ and the decimal part of x_j remains unchanged.

Application of SAEGPSO-based allocation method on the MRTAP

The pseudocode of SAEGPSO for solving the MRTAP is shown in Table 1, where $t_{\text{max}}$ denotes the maximum iteration number of SAEGPSO.

Table 1.

The SAEGPSO-based framework for MRTAP.

1. Randomly generate an initial swarm

2. Evaluate the initial swarm and obtain $\text{Pbest}$ and $\text{Gbest}$ , as well as $\text{ESS}$ for the initial swarm

3. while $t\le t_{\text{max}}$ do

4. Get the previous experience, including $Z_1(t-1)$ , $Z_2(t-1)$ , $Z_3(t-1)$ , and the fitness value $F(X(t-1))$ , of each particle

5. Calculate $\pi (e_1)$ , $\pi (e_2)$ , and $\pi (e_3)$ based on equation (15)

6. Compute the payoff matrix K based on equation (12)

7. Calculate $Z_1(t)$ , $Z_2(t)$ , and $Z_3(t)$ based on the RDE given by equation (11)

8. Calculate the three control parameters ω, c1, and c2 of each particle based on equations (16) to (22)

9. Calculate the barycenter $\text{BC}(t)$ of each particle based on equation (11)

10. Randomly obtain a point $X^{\prime }(t)$ within the hypersphere $\text{HP}(\text{BC}(t),||\text{BC}(t)-X(t)||)$ for each particle

11. Update the velocity of each particle based on equation (7)

12. Update the position of each particle based on equation (8)

13. Modify the position vector of each particle based on the saturation strategy given by equation (51)

14. Calculate the fitness value of each particle based on equations (1) and (2)

15. Calculate the constraint violation degree of each particle based on equation (50)

16. Update $\text{Bbest}$ of each particle based on the feasibility-based rule

17. Update $\text{Gbest}$ of the swarm based on the feasibility-based rule

18. Increase the iteration number t by 1

19. end while

20. Output $\text{Gbest}$

It is noticeable from Table 1 that, similar to ω_s, ω_f, $c_{1s}$ , $c_{1f}$ , $c_{2s}$ , and $c_{2f}$ , the value of $t_{\text{max}}$ is empirically predefined by the users. Also, please note that how to empirically set the values of these parameters may be problem-dependent, which requires the previous experience of the users. However, it is important to note that settings of these parameters must abide by several generally rules of PSO development. Firstly, settings of ω_s, ω_f, $c_{1s}$ , $c_{1f}$ , , $c_{2s}$ , and $c_{2f}$ must meet the convergence condition of SAEGPSO given by equation (47) in order to sufficiently guarantee the convergence of SAEGPSO. Secondly, in order to facilitate the exploration ability of SAEGPSO in the early stages of the evolution and promote the exploitation ability of SAEGPSO in the latter stages of the evolution, it has to satisfy that ${\omega }_s>{\omega }_f$ , $c_{1s}>c_{1f}$ , and $c_{2f}>c_{2s}$ in the newly proposed self-adaptive strategy presented in SAEGPSO. Thirdly, the setting of $t_{\text{max}}$ must consider the solution optimality and the computation time of SAEGPSO in a compromising way, since $t_{\text{max}}$ has influences on the solution optimality and computation time of SAEGPSO. The bigger value $t_{\text{max}}$ has, the higher-quality solution SAEGPSO may find; however, the more computation time SAEGPSO has to consume. Since the decision makers may prefer to obtain a satisfied a solution with low computation cost, $t_{\text{max}}$ is empirically set to be 2000 when using SAEGPSO to solve the MRTAP in this paper. Also, in order to guarantee the convergence of SAEGPSO and well balance the trade-offs between exploration and exploitation, when using SAEGPSO to solve the MRTAP, we empirically set that ${\omega }_s=0.9$ , ${\omega }_f=0.1$ , $c_{1s}=c_{2f}=2.5$ , and $c_{1f}=c_{2s}=0.1$ based on the analysis results in ‘Convergence-guaranteed parameter selection principle for SAEGPSO’ section.

Simulations, comparisons, and analysis

In this subsection, SAEGPSO is evaluated through three different-sized task allocation cases against mDE, ³ GA, ⁹ SPSO 2011, ¹⁸ and ILDWPSO. ²⁶ For each allocation case, after each method conducting a Monte Carlo experiment with 30 runs, this study reports and compares the statistical results with respect to the solution optimality and computation time of each method. In every single run of the Monte Carlo test, the best solution of each method is output after 40 particles evolve 2000 iterations. The simulation parameters for SAEGPSO are set as follows: ${\omega }_s=0.9$ , ${\omega }_f=0.1$ , $c_{1s}=c_{2f}=2.5$ , and $c_{1f}=c_{2s}=0.1$ based on the analysis results in ‘Convergence-guaranteed parameter selection principle for SAEGPSO’ section. The simulation parameters for the other compared methods are exacted from their corresponding literature and summarized in Table 2.

Table 2.

The simulation parameters for different compared methods.

Methods	Parameter setting
GA 9	$\text{CR}=0.8$
mDE 3	${\text{CR}}_{\text{min}}=0.1$ , ${\text{CR}}_{\text{max}}=0.9$
SPSO 2011 18	$\omega =0.7213$ , $c_1=1.1931$ , $c_2=1.1931$
ILDWPSO 26	$\alpha =4$ , ${\omega }_{\text{start}}=0.9$ , ${\omega }_{\text{end}}=0.1$ , and ${\mu }_{\text{m}}=0.02$

Case study 1

This subsection shows the simulation results of all methods in solving the first allocation instance which assigns three robots to 12 tasks. The properties of robots and tasks are given and shown in Tables 3 and 4, respectively. After conducting the Monte Carlo experiment, Tables 5 and 6 report the statistical results of the solution and the computation time of each method, respectively. In these two tables, the best mean results with respect to the solution and computation time are highlighted in boldface. For concretely showing the allocation results and better understanding the evolutionary processes of different methods, the allocation schemes and fitness curves of the best solutions obtained by different methods in the Monte Carlo test are exhibited in Table 7 and Figure 7, respectively. The robots’ trajectories of the best solutions of SAEGPSO and GA are illustrated in Figure 8, respectively. Note that because positions of robots and tasks are given in Tables 3 and 4, and the allocation schemes of the best solutions of all methods are shown in Table 7, the robots’ trajectories of the best solutions of the other methods (except SAEGPSO and GA) are not displayed in order not to enlarge the size of the paper.

Table 3.

Properties of robots for the first studied case.

Robot index	Robot location (m)	Tactical value	Execution probability	Payload constraint
1	(50,100)	76	0.55	3
2	(50,200)	100	0.78	4
3	(50,300)	156	0.85	5

Table 4.

Properties of tasks for the first studied case.

Task index	Task location (m)	Tactical value	Vio-execution probability	Task index	Task location	Tactical value	Vio-execution probability
1	(100,50)	120	0.45	7	(153,285)	123	0.40
2	(94,152)	110	0.40	8	(142,381)	101	0.35
3	(102,252)	105	0.43	9	(200,20)	98	0.35
4	(95,363)	123	0.51	10	(210,120)	109	0.45
5	(150,75)	104	0.45	11	(202,240)	123	0.32
6	(145,185)	112	0.38	12	(230,340)	111	0.30

vio-execution probability: violated execution probability of the task to robots.

Table 5.

Simulation results of the allocated solutions obtained by each method for the first studied case.

	SAEGPSO	ILDWPSO	SPSO 2011	mDE	GA
Best	1.8131E+02	1.8796E+02	2.7455E+02	1.8388E+02	2.3649E+02
Worst	3.8292E+02	4.4714E+02	9.9014E+02	5.3785E+02	8.3413E+02
Mean	2.5549E+02	3.2747E+02	6.3294E+02	4.0792E+02	4.9294E+02
Std.	5.3881E+01	7.1334E+01	1.7039E+02	9.4914E+01	1.5759E+02

Std.: standard deviation. The best mean results with respect to the solution and the computation time among the four methods are highlighted in boldface.

Table 6.

Simulation results of the computation time (s) of each method for the first studied case.

	SAEGPSO	ILDWPSO	SPSO 2011	mDE	GA
Best	32.82	37.89	21.98	36.33	33.39
Worst	37.47	44.38	33.43	43.45	37.13
Mean	35.65	41.47	29.55	40.17	36.13
Std.	3.401	3.001	4.394	4.031	4.781

Std.: standard deviation. The best mean results with respect to the solution and the computation time among the four methods are highlighted in boldface.

Table 7.

Allocation schemes of the best solutions searched by different methods for the first studied case.^a

Methods	Allocation scheme of individual robot
SAEGPSO	$R_1:T_1\to T_5\to T_9$
	$R_2:T_2\to T_6\to T_{11}\to T_{10}$
	$R_3:T_3\to T_7\to T_{12}\to T_8\to T_4$
ILDWPSO	$R_1:T_1\to T_5\to T_9$
	$R_2:T_2\to T_6\to T_{11}\to T_{10}$
	$R_3:T_4\to T_8\to T_{12}\to T_7\to T_3$
SPSO 2011	$R_1:T_1\to T_5\to T_9$
	$R_2:T_2\to T_4\to T_8\to T_{12}$
	$R_3:T_3\to T_7\to T_{11}\to T_6\to T_{10}$
mDE	$R_1:T_1\to T_5\to T_9$
	$R_2:T_2\to T_6\to T_{11}\to T_{10}$
	$R_3:T_3\to T_7\to T_4\to T_8\to T_{12}$
GA	$R_1:T_1\to T_5\to T_2$
	$R_2:T_3\to T_6\to T_{10}\to T_9$
	$R_3:T_4\to T_8\to T_{12}\to T_7\to T_{11}$

^aHere, R_i and T_j denote the i th robot and j th task, and “ $T_A\to T_B$ ”means that Task B is executed after Task A.

Figure 7.

The fitness curves of the best runs of different methods for the first studied case.

Figure 8.

Robots’ trajectories of the best solutions searched by SAEGPSO and GA for the first studied case. (a) Robots’ trajectories of the best solution searched by SAEGPSO. (b) Robots’ trajectories of the best solution searched by GA. SAEGPSO: self-adaptive evolutionary game particle swarm optimization; GA: genetic algorithm.

From Table 5, the proposed SAEGPSO method is followed by ILDWPSO, mDE, GA, and SPSO 2011 in terms of the average solution optimality, which implies that SAEGPSO generally outperforms its contenders in terms of the solution optimality in solving the first test case. From Table 5, comparing with ILDWPSO, mDE, GA, and SPSO 2011, SAEGPSO averagely improves the solution optimality by 21.98%, 37.37%, 48.17%, and 59.63%, respectively. From Table 6, one can easily observe that SPSO 2011 is the most efficient method, while SAEGPSO method is ranked second in terms of the average computation time. Such observation is reasonable because the three control parameters of SPSO 2011 remain constant, no additional computation resource is needed to update these control parameters of SPSO 2011. Hence, in comparison with SAEGPSO, the computation time of SPSO 2011 is naturally reduced. However, it is important to note from Tables 5 and 6 that although SPSO 2011 is the most efficient method in the computation time, it obtains the worst performance in terms of the solution quality, and even if SAEGPSO is the second most efficient method as far as the computation time is concerned, it has the best performance in terms of solution quality. Since the studied MRTAP is completed off-line, the computation time may weigh less than the solution optimality in evaluating a task allocation optimizer. Thus, considering all of the concerns raised, it allows us to conclude that, comparing with the other four methods, SAEGPSO is highly competitive in solving the first test case. Moreover, from Table 7, it can be easily observed that the best solutions obtained by all methods are feasible since all the best solutions satisfy the constraints of the first studied allocation case, which, to some extent, can reveal the feasibility of these methods in solving the small-scale MRTAP.

Case study 2

The second task allocation instance is a medium-scale TAP, in which five robots are assigned to complete 22 tasks. The properties of robots and tasks are given and shown in Tables 8 and 9, respectively. The statistical results with respect to the solution and computation time of different methods for this allocation case are summarized in Tables 10 and 11, respectively. With the same reasons explained in ‘Case study 1’ section, the allocation schemes and fitness curves of the best solutions searched by all considered methods are displayed in Table 12 and Figure 9, respectively.

Table 8.

Properties of robots for the second studied case.

Robot index	Robot location (m)	Tactical value	Execution probability	Payload constraint
1	(16,11)	20	0.40	3
2	(32,21)	46	0.55	5
3	(45,65)	55	0.60	6
4	(72,64)	41	0.50	4
5	(11,77)	31	0.45	4

Table 9.

Properties of tasks for the second studied case.

Task index	Task location (m)	Tactical value	Vio-execution probability	Task index	Task location(m)	Tactical value	Vio-execution probability
1	(10,220)	50	0.38	12	(132,54)	37	0.58
2	(32,157)	37	0.45	13	(263,157)	40	0.45
3	(60,173)	39	0.59	14	(280,38)	37	0.34
4	(72,181)	40	0.50	15	(102,241)	36	0.46
5	(98,189)	35	0.47	16	(216,144)	42	0.44
6	(176,184)	46	0.40	17	(220,156)	46	0.55
7	(61,162)	44	0.60	18	(224,179)	47	0.67
8	(368,178)	42	0.57	19	(228,192)	50	0.37
9	(300,70)	41	0.60	20	(232,214)	37	0.57
10	(296,97)	41	0.65	21	(234,242)	44	0.44
11	(288,192)	38	0.40	22	(246,224)	46	0.41

vio-execution probability: violated execution probability of the task to robots.

Table 10.

Simulation results of the allocated solutions obtained by each method for the second studied case.

	SAEGPSO	ILDWPSO	SPSO 2011	mDE	GA
Best	1.7721E+03	1.7848E+03	1.9678E+03	1.8843E+03	2.0112E+03
Worst	2.3078E+03	2.3854E+03	2.7121E+03	2.5108E+03	2.4547E+03
Mean	1.9949E+03	2.1239E+03	2.3103E+03	2.1625E+03	2.2625E+03
Std.	1.2752E+02	1.6013E+02	1.9789E+02	1.8154E+02	1.2339E+02

The best mean results with respect to the solution and the computation time among the four methods are highlighted in boldface.

Table 11.

Simulation results of the computation time (seconds) of each method for the second studied case.

	SAEGPSO	ILDWPSO	SPSO 2011	mDE	GA
Best	51.95	59.36	49.35	59.17	50.45
Worst	62.34	68.43	57.25	65.35	64.73
Mean	56.78	65.48	51.59	63.71	58.47
Std.	3.892	3.744	4.064	3.021	4.715

The best mean results with respect to the solution and the computation time among the four methods are highlighted in boldface.

Table 12.

Allocation schemes of the best solutions searched by different methods for the second studied case.

Methods	Allocation scheme of individual robot
SAEGPSO	$R_1:T_{18}\to T_{19}\to T_{11}$
	$R_2:T_{14}\to T_9\to T_{10}\to T_{17}\to T_8$
	$R_3:T_2\to T_6\to T_{20}\to T_{21}\to T_{22}\to T_{13}$
	$R_4:T_1\to T_4\to T_5\to T_{16}$
	$R_5:T_{12}\to T_7\to T_3\to T_{15}$
ILDWPSO	$R_1:T_{17}\to T_{13}\to T_{21}$
	$R_2:T_7\to T_3\to T_6\to T_{11}\to T_8$
	$R_3:T_2\to T_{12}\to T_{14}\to T_9\to T_{10}\to T_{18}$
	$R_4:T_{20}\to T_{22}\to T_{16}\to T_{19}$
	$R_5:T_4\to T_5\to T_{15}\to T_1$
SPSO 2011	$R_1:T_2\to T_3\to T_{20}$
	$R_2:T_{10}\to T_{17}\to T_{18}\to T_{22}\to T_{11}$
	$R_3:T_1\to T_{15}\to T_4\to T_7\to T_5\to T_6$
	$R_4:T_{12}\to T_{14}\to T_{13}\to T_8$
	$R_5:T_9\to T_{16}\to T_{19}\to T_{21}$
mDE	$R_1:T_{16}\to T_{13}\to T_{14}$
	$R_2:T_{21}\to T_{20}\to T_{19}\to T_{18}\to T_{15}$
	$R_3:T_7\to T_4\to T_6\to T_{22}\to T_{11}\to T_{10}$
	$R_4:T_{12}\to T_3\to T_2\to T_1$
	$R_5:T_5\to T_{17}\to T_9\to T_8$
GA	$R_1:T_{13}\to T_{22}\to T_{21}$
	$R_2:T_{20}\to T_{11}\to T_8\to T_{18}\to T_{14}$
	$R_3:T_2\to T_3\to T_{17}\to T_{19}\to T_9\to T_{10}$
	$R_4:T_{12}\to T_{16}\to T_6\to T_4$
	$R_5:T_5\to T_7\to T_{15}\to T_1$

Figure 9.

The fitness curves of the best runs of different methods for the second studied case.

From Table 10, it is clear that SAEGPSO and ILDWPSO are the most and second most efficient approaches, while SPSO 2011 obtains the worst performance in terms of the solution quality for the second simulation scenario. Based on the results shown in Table 10, comparing with ILDWPSO, mDE, GA, and SPSO 2011, SAEGPSO generally enhances the solution optimality by 6.07%, 9.28%, 11.83%, and 13.65%, respectively. From Table 11, SPSO 2011 and SAEGPSO are ranked the first and second in terms of the computation time, respectively. It is notable from Tables 10 and 11 that, despite having the best performance in the computation time, SPSO 2011 is the least efficient method in the solution optimality, and even if SAEGPSO is ranked second in terms of the computation time, it is the most efficient approach among the five methods in terms of the solution quality. Hence, we can still conclude that SAEGPSO is highly competitive in solving the second allocation instance. Moreover, it can easily seen from Table 12 that the best solutions found by the five methods are feasible because these best solutions can satisfy all the constraints of the second allocation case, which, to a certain degree, indicates the feasibility of these methods in the medium-scale MRTAP.

From Figure 9, it is interesting to observe that the best fitness curves of some methods such as GA first rise in the early of the evolution and then drop or remain stable in the later stages of the evolution. This observation is due to that we use the feasibility-based rule ²⁵ to choose the dominated solution between any two candidate solutions in the SAEGPSO-based task allocation method. Because the initial population of each method is randomly produced for any simulation scenario, it may happen that all the initial individuals of a method violate some constraints of this scenario, which means that all individuals of the method are initially non-feasible. In this case, the initial global best solution of the method is thus also non-feasible. As the evolutionary process of the method carrying on, if any particle of the method can find a solution having a smaller constraint violation degree but a larger fitness value than the initial non-feasible global best solution, the newly produced solution will then replace the initial non-feasible global best solution according to the feasibility-based rule. Therefore, in such a case, the fitness curves of the method first rise in the early stages of the evolutionary process. However, this phenomenon will vanish in the later stages of the evolutionary process, namely the fitness curves of the method either keep dropping or stable in the latter stages of the evolutionary process. This is because, with the evolution of the method continuing, some individuals of the method can find some feasible solutions in the latter stages of the evolution. In this case, after these newly found feasible solutions replace the historically non-feasible global best solution of the method, the current global best solution of the method becomes feasible, which can be only replaced by the feasible solutions that have less fitness values than the current feasible global best solution based on the feasibility-based rule. Therefore, as shown in Figure 9, the fitness curves of the five methods keep dropping or stable in the latter phases of the evolution.

Note that, since the first test case studied in ‘Case study 1’ section is a small-scale TAP, in which only 12 tasks are involved, it has a high chance that each method can obtain a feasible initial global best solution after randomly initializing the swarm of each method. Because all the initial global best solutions of all methods are feasible in the first test case studied in ‘Case study 1’ section, all the initial global best solutions of these methods can be only replaced by newly produced feasible solutions with smaller fitness values according to the feasibility-based rule. Therefore, as shown in Figure 7, the best fitness curves of the five methods first keep going down in the early stages of the evolution and then remain unchanged in the latter stages of the evolution.

Case study 3

The third test case is a relatively large-sized allocation instance, in which six robots are assigned to 40 tasks. The properties of robots and tasks are shown in Tables 13 and 14, respectively. Tables 15 and 16 show the statistical results of the allocation solution and the computation time of different methods for this studied case, respectively. The allocation results and fitness curves of the best runs of different methods in the Monte Carlo experiments conducted for this scenario are illustrated in Table 17 and Figure 10, respectively.

Table 13.

Properties of robots for the third studied case.

Robot index	Robot Location (m)	Tactical value	Execution probability	Payload constraint
1	(127.21,247.06)	24.75	0.41	7
2	(88.310,162.41)	22.62	0.53	9
3	(76.470,35.030)	26.81	0.54	4
4	(292.07,239.62)	37.44	0.23	6
5	(385.48,220.35)	20.06	0.40	6
6	(260.95,164.27)	25.68	0.60	8

Table 14.

Properties of tasks for the third studied case.

Task index	Task location (m)	Tactical value	Vio-execution probability	Task index	Task location (m)	Tactical value	Vio-execution probability
1	(96.670,147.46)	22.37	0.51	21	(9.3700,243.09)	56.76	0.29
2	(104.59,278.15)	33.06	0.53	22	(375.11,392.59)	60.95	0.24
3	(328.27,360.49)	65.41	0.37	23	(182.12,282.32)	66.69	0.55
4	(334.52,249.93)	41.78	0.42	24	(3.4100,380.60)	23.04	0.50
5	(37.630,348.50)	56.89	0.49	25	(351.78,331.67)	67.47	0.29
6	(377.11,351.23)	53.09	0.20	26	(387.35,165.45)	56.34	0.30
7	(276.82,311.36)	48.75	0.57	27	(256.69,384.07)	45.31	0.45
8	(378.02,143.25)	48.27	0.28	28	(162.23,260.93)	47.84	0.35
9	(56.372,226.44)	69.77	0.58	29	(281.55,52.870)	26.73	0.52
10	(217.38,326.58)	48.92	0.51	30	(28.590,45.000)	22.86	0.26
11	(258.85,259.92)	68.65	0.52	31	(365.16,284.31)	52.77	0.44
12	(275.511,34.98)	52.09	0.58	32	(298.72,246.81)	48.09	0.26
13	(208.69,378.02)	33.21	0.57	33	(179.91,268.38)	20.62	0.52
14	(371.49,198.37)	68.53	0.53	34	(224.13,171.77)	35.88	0.36
15	(214.94,1.8500)	42.40	0.31	35	(211.52,325.43)	53.40	0.45
16	(255.20,35.350)	32.40	0.28	36	(240.49,177.87)	23.83	0.50
17	(334.53,375.83)	42.74	0.59	37	(161.53,188.31)	68.27	0.20
18	(201.80,323.61)	39.88	0.36	38	(164.03,396.11)	45.58	0.36
19	(291.61,6.8100)	33.93	0.37	39	(249.02,9.4400)	50.13	0.38
20	(144.87,4.9100)	39.65	0.20	40	(106.64,55.120)	69.57	0.40

vio-execution probability: violated execution probability of the task to robots.

Table 15.

Simulation results of the allocated solutions obtained by each method for the third studied case.

	SAEGPSO	ILDWPSO	SPSO 2011	mDE	GA
Best	4.2863E+03	4.77746E+03	5.8868E+03	4.7733E+03	5.4734E+03
Worst	5.7021E+03	5.9897E+03	6.6846E+03	6.5341E+03	6.6271E+03
Mean	5.2082E+03	5.2522E+03	6.2577E+03	5.5713E+03	5.9449E+03
Std.	2.6677E+02	3.9779E+02	2.1438E+03	3.6886E+02	3.8411E+03

The best mean results with respect to the solution and the computation time among the four methods are highlighted in boldface.

Table 16.

Simulation results of the computation time (s) of each method for the third studied case.

	SAEGPSO	ILDWPSO	SPSO 2011	mDE	GA
Best	123.08	137.62	113.45	134.42	126.13
Worst	131.71	147.31	120.36	142.49	133.45
Mean	126.41	141.45	115.61	136.35	128.14
Std.	3.0728	3.6113	2.4365	2.9489	3.5124

The best mean results with respect to the solution and the computation time among the four methods are highlighted in boldface.

Table 17.

Allocation schemes of the best solutions searched by different methods for the third studied case.

Methods	Allocation scheme of individual robot
SAEGPSO	$R_1:T_{37}\to T_{34}\to T_{16}\to T_8\to T_{26}\to T_{31}\to T_{22}$
	$R_2:T_1\to T_9\to T_{21}\to T_2\to T_{27}\to T_{35}\to T_{23}\to T_{32}\to T_6$
	$R_3:T_{30}\to T_{20}\to T_{39}\to T_{40}$
	$R_4:T_{14}\to T_7\to T_{24}\to T_{38}\to T_{13}\to T_3$
	$R_5:T_5\to T_{12}\to T_{19}\to T_4\to T_{25}\to T_{18}$
	$R_6:T_{28}\to T_{33}\to T_{15}\to T_{29}\to T_{36}\to T_{11}\to T_{10}\to T_{17}$
ILDWPSO	$R_1:T_{38}\to T_{24}\to T_1\to T_{37}\to T_{29}\to T_{19}\to T_{26}$
	$R_2:T_{21}\to T_{30}\to T_{33}\to T_{28}\to T_{18}\to T_{13}\to T_{25}\to T_{12}\to T_{39}$
	$R_3:T_{40}\to T_{15}\to T_{34}\to T_3$
	$R_4:T_{31}\to T_7\to T_{11}\to T_{23}\to T_{35}\to T_{22}$
	$R_5:T_{32}\to T_{20}\to T_4\to T_6\to T_{17}\to T_{27}$
	$R_6:T_{36}\to T_5\to T_9\to T_2\to T_{16}\to T_8\to T_{14}\to T_{10}$
SPSO 2011	$R_1:T_{26}\to T_{12}\to T_{29}\to T_5\to T_{23}\to T_3\to T_{21}$
	$R_2:T_1\to T_{33}\to T_{32}\to T_{31}\to T_{14}\to T_7\to T_{18}\to T_{10}\to T_6$
	$R_3:T_4\to T_{20}\to T_{30}\to T_{36}$
	$R_4:T_{11}\to T_{25}\to T_9\to T_{40}\to T_{17}\to T_8$
	$R_5:T_{27}\to T_{13}\to T_{24}\to T_{38}\to T_{35}\to T_{22}$
	$R_6:T_2\to T_{39}\to T_{15}\to T_{16}\to T_{19}\to T_{37}\to T_{34}\to T_{28}$
mDE	$R_1:T_{33}\to T_7\to T_{22}\to T_{25}\to T_4\to T_{11}\to T_{23}$
	$R_2:T_{37}\to T_{27}\to T_{32}\to T_{31}\to T_{17}\to T_{28}\to T_{10}\to T_{38}\to T_5$
	$R_3:T_{15}\to T_{29}\to T_{36}\to T_{24}$
	$R_4:T_{12}\to T_{19}\to T_{39}\to T_{40}\to T_8\to T_6$
	$R_5:T_{21}\to T_9\to T_{30}\to T_{16}\to T_{20}\to T_3$
	$R_6:T_{26}\to T_{14}\to T_{34}\to T_{18}\to T_{35}\to T_{13}\to T_2\to T_1$
GA	$R_1:T_{35}\to T_{23}\to T_{22}\to T_{31}\to T_8\to T_7\to T_{11}$
	$R_2:T_{15}\to T_{39}\to T_{16}\to T_{26}\to T_{34}\to T_{27}\to T_2\to T_1\to T_{21}$
	$R_3:T_{17}\to T_3\to T_4\to T_{36}$
	$R_4:T_{30}\to T_{13}\to T_{38}\to T_{24}\to T_5\to T_{28}$
	$R_5:T_6\to T_{25}\to T_{32}\to T_{12}\to T_{10}\to T_{33}$
	$R_6:T_{20}\to T_9\to T_{18}\to T_{14}\to T_{37}\to T_{40}\to T_{29}\to T_{19}$

Figure 10.

The fitness curves of the best runs of different methods for the third studied case.

From Tables 15 and 16, it is apparent that the proposed SAEGPSO is followed by ILDWPSO, mDE, GA, and SPSO 2011 in terms of the average solution quality, and SPSO 2011 is followed by SAEGPSO, GA, mDE, and ILDWPSO in the average computation time. Both considering the statistical results of the solution quality and the computation time shown in Tables 15 and 16, we can conclude that SAEGPSO is highly powerful in solving the large-size MRTAP. From Table 17, it can easily noted that the best solutions obtained by all methods are feasible since all the best solutions trivially meet the constraints of the third task allocation case, which, more or less, can reflect the feasibility of each considered method in solving the large-scale MRTAP. With the same reasons depicted in ‘Case study 2’ section, as shown in Figure 10, the best fitness curves of mDE and SPSO 2011 first rise in the early stages of the evolution and then keep going down in the latter stages of the evolution.

Conclusion and future work

Focusing on enhancing the performance of PSO via overcoming two typical drawbacks of the standard PSO, this study proposes a novel enhanced SAEGPSO. To this end, we first let particles update their movements based on the moving rules defined in SPSO 2011 to avoid stagnation in SAEGPSO. Afterwards, a self-adaptive strategy determined by the evolutionary stable strategies of EGT and the iteration number of SAEGPSO is proposed for fine-tuning the control parameters of particles in order to well trade off the exploration and exploitation capabilities of SAEGPSO. To properly set the self-adaptive strategy, this paper also theoretically investigates the convergence of SAEGPSO with respect to values of the three control parameters and provides a convergence-guaranteed parameter selection principle for SAEGPSO.

Leveraging the development of SAEGPSO, this paper also develops a SAEGPSO-based allocation method for solving the MRTAP. In order to efficiently solve the MRTAP and reduce the optimization difficulties, the feasibility-based rule is adopted in the developed SAEGPSO-based allocation method to handle constraints of the MRTAP. In addition, the performance of the proposed method is verified through three different-scaled task allocation cases against four well-established evolutionary methods. The simulation results confirm that SAEGPSO outperforms its contenders in terms of the solution quality. Moreover, SAEGPSO is also highly powerful in terms of the computation time. Therefore, it allows us to conclude that SAEGPSO is highly promising in solving the MRTAP and can be considered as vital alternative in the field of MRTAP.

There are a few potential issues that deserve some future study. Firstly, the performance of SAEGPSO can be further enhanced through improving the EGT. Secondly, since the three control parameters not only affect the convergence of PSO, but also affect the convergence speed of PSO, we will study the convergence speed of SAEGPSO with respect to different control parameter settings in the future. Moreover, despite investigating the convergence to an equilibrium point of SAEGPSO, the local or global optimality of the equilibrium point has not been investigated in this study, which will be addressed in the near future. Thirdly, in order to further evaluate SAEGPSO in the MRTAP, we are also considering the possibility of comparing SAEGPSO with some other non-evolutionary methods such as the auction-based methods in the near future. Last but not least, we are considering to use the proposed PSO method to solve some other COP such as environmental/economic dispatch problem ⁴³ and TAP ⁴⁴ in the near future.