Cooperative search for moving targets with the ability to perceive and evade using multiple UAVs

2.1. Motion model of UAVs

Mission area $$ E $$ is divided into $$ N_x\text{×}N_y $$ grid cells with $$ L_U $$ as the side length of a single cell, which are used as mission grid or UAV motion grid and denoted as $$ \Omega^U $$. It is agreed that the UAV has limited movement between mission grid cells, where the cell with m-th row and n-th column in $$ \Omega^U $$ is marked as $$ g_{mn}^U $$. Figure 1 shows the specific representation of the mission grid.

Figure 1. Definition of mission grid.

Based on this definition, the multi-UAV system executing cooperative search tasks can be considered as a complete control system. The state variables of the UAVs at time $$ t_k $$ can be represented as $$ X(t_k)=\{X_1(t_k),X_2(t_k),…,X_i(t_k),…,X_{N_U}(t_k)\} $$, where $$ X_i(t_k),i=1,2,…,N_U $$ represents the position at time $$ t_k $$ in $$ \Omega^U $$ of i-th UAV. Each UAV can move to the adjacent cells or stay in the current cell every motion cycle of $$ {\Delta}T $$ time. Figure 2 shows the definition of flight directions of UAVs.

Figure 2. Definition of flight directions of UAV.

The control input of i-th UAV at time $$ t_k $$ can be represented as $$ u_i(t_k)\in\{1,2,3,4,5,6,7,8,9\} $$, where $$ i=1,2,…,N_U $$; flight direction "9" represents staying in the current cell, and flight directions 1-8 represent moving one step towards adjacent cells. Therefore, the control input of the multi-UAV system at time $$ t_k $$ can be represented as $$ U(t_k)=\{u_1(t_k),u_2(t_k),…,u_{N_U}(t_k)\} $$, and the state model of the multi-UAV system can be recorded as $$ X(t_{k+1})=\mathrm{F}(X(t_k),U(t_k )) $$, where $$ \mathrm{F}(\cdot) $$ is the state transition function of the system.

2.2. Motion model of moving targets

As the motion law of the moving targets is unknown, it is assumed that the moving target has a maximum movement distance of $$ L_{Tmax} $$, and the motion of each target can be represented by a Gauss-Markov motion model^[9]. The movement of each target in every motion cycle can be considered as a combination of random movement and directional evasive movement. Figure 3 shows the movement of the moving target with different evasion intensities.

Figure 3. Combination of movement with different evasion intensities.

In order to effectively evade the search of UAVs, a virtual evasive force is defined in this paper, which includes two directions of evasion: away from the current position of the UAV and away from the forward path of the UAV. Figure 4 shows the virtual evasive force of a moving target when facing the search for i-th UAV.

Figure 4. Virtual evasive force of moving targets.

In Figure 4, $$ {\bf{F}}_{ia} $$ represents the virtual evasive force A to keep away from the current position of i-th UAV, $$ {\bf{F}}_{ib} $$ represents the virtual evasive force B to keep away from the forward path of i-th UAV, and $$ {\bf{F}}_i $$ represents the comprehensive virtual evasive force of i-th UAV. Considering the distance from the target to the UAV, the virtual evasive force A is defined as:

where $$ {\bf{e}}_{ia} $$ is the unit direction vector from i-th UAV to the target, $$ A_a $$ is the response amplitude of virtual evasive force A, and $$ \delta_a $$ is the corresponding straight response distance. Considering that the target cannot effectively determine the flight direction and threat level of the UAV at a long distance, it is reasonable to assume that the evasive force is negatively correlated with distance. The virtual evasive force B is defined as:

where $$ {\bf{e}}_i $$ is the unit direction vector of the flight direction of i-th UAV, $$ {\bf{e}}_{ib} $$ is the unit direction vector pointing vertically to the target from the direction of $$ {\bf{e}}_i $$, $$ A_b $$ is the response amplitude of virtual evasive force B, and $$ \delta_{b1} $$ and $$ \delta_{b2} $$ represent the corresponding response distances in straight and vertical directions, respectively. This assumption reasonably estimates the target's evasive behavior toward the UAV's pursuit. The final evasive virtual force of the moving target can be represented as:

In summary, based on the Gauss-Markov motion model, it is assumed that the motions in the x and y axes are independent of each other, the transition probability density function of the target from the position $$ (x_k,y_k) $$ to the position $$ (x_{k+1},y_{k+1}) $$ after one motion cycle $$ {\Delta}T $$ can be represented as:

In order to appropriately allocate the proportion of random movement and the directional evasive movement, the hyperbolic tangent (tanh) function is used to saturate the virtual evasive force of the target in this paper. Therefore, the means of $$ \mu_x $$ and $$ \mu_y $$ and variance of $$ \sigma_{xy}^2 $$ can be defined as:

where $$ \omega_{sat}\in[0,1] $$ is the saturated partition coefficient of evasive movement. According to the common definition of the $$ 3\sigma $$ rule, $$ L_{Tmax} $$ is appropriately allocated by the guidance of virtual evasive force, and consequently, the evasive maneuver of targets can be realized.

2.3. Detection model of UAV

The UAV detection range is defined as a square area with a side length of $$ L_D $$ centered on the UAV's current position, which can be represented as the mission grid cell and its adjacent cells where the UAV is located. Considering the real-time changes in UAV position, denote the detection grid of i-th UAV as $$ \Omega^{Di} $$, which can be regarded as a certain range in $$ \Omega^U $$, and its side length is $$ \gamma_D $$ times that of a single mission grid cell. Same as above, the detection grid cell of i-th UAV with m-th row and n-th column is marked as $$ g_{mn}^{Di} $$. Figure 5 shows the specific representation of the detection grid of i-th UAV.

Figure 5. Definition of detection grid of i-th UAV.

The i-th UAV can only detect the area within $$ \Omega^{Di} $$. The detection is carried out in every single detection grid cell, and the detection probability corresponding to $$ g_{mn}^{Di} $$ is $$ p_{mn}^D $$, which can be defined as a general form^[25] as:

where $$ p_{max}^D $$ and $$ p_{min}^D $$ are the maximum and minimum detection probabilities, and $$ \delta_{Dmin} $$ and $$ \delta_{Dmax} $$ are their corresponding boundary detection distances. In order to comply with the actual detection situation, it is assumed that the probability of a false alarm occurring during the detection process is $$ p^F $$; that is, during a motion cycle $$ {\Delta}T $$, UAV may generate a false alarm signal within the detection range with a probability of $$ p^F $$.

Moreover, the additional conventions are given here:

(ⅰ) To avoid collisions between UAVs, at most one UAV is allowed to search within each mission grid cell at any time instance.

(ⅱ) Communications between UAVs are ideal, and there is no communication delay.

(ⅲ) There is at most one target in a single mission grid cell.

(ⅳ) After discovering the target, the UAV will detach from the search sequence, enter a tracking state, and will not participate in subsequent search missions.

3.1. Probability map

A probability map is a search information map that reflects the existence probability of targets. Considering the difference in speed between a UAV and a moving target, the mission grid $$ \Omega^U $$ is refined and divided as the region division of the probability map. After division, the mission area is divided into $$ M_x\text{×}M_y $$ grid cells with $$ L_T $$ as the side length of a single grid cell where $$ L_{Tmax}{\le}L_T $$, and the side length of a mission grid cell is $$ \gamma_T $$ times of $$ L_T $$. Denote the grid of probability map as $$ \Omega^T $$, in which the cell with m-th row and n-th column is marked as $$ g_{mn}^T $$. Figure 6 shows the specific representation of the grid definition of a probability map.

Figure 6. Grid definition of a probability map.

Denote the existence probability of the target in $$ g_{mn}^T $$ as $$ p_{mn} $$. Considering that the detection is carried out in a single detection grid cell, i.e., a mission grid cell, if the existence of a target in the mission grid cell with r-th row and c-th column is denoted as $$ T_{rc} $$, then the probability of $$ T_{rc} $$ is:

where $$ g_{mn}^T{\in}g_{rc}^U $$ denotes that $$ g_{mn}^T $$ is within the area of $$ g_{rc}^U $$ in terms of regional division, and the definition of such notations in the rest of the paper is the same as here. It can be determined that the target is found in $$ g_{rc}^U $$ when $$ P(T_{rc}) $$ is not less than a certain threshold $$ p_{thold}^D $$.

The transition of the existence probability based on target motion prediction under the influence of the movement of UAVs in any grid cell of the probability map can be represented by the probability diffusion coefficient of one cell to another. Denote the probability diffusion coefficient for transferring the existence probability of the target from $$ g_{rc}^T $$ to $$ g_{mn}^T $$ at time $$ t_k $$ as $$ \Phi_{rc}^{mn}(t_k) $$, which can be expressed based on target motion prediction with Equations (4) and (5) as follows:

where $$ (Y_r,X_c) $$ and $$ (Y_m,X_n) $$ are the centers of $$ g_{rc}^T $$ and $$ g_{mn}^T $$. Obviously, the sum of the probability diffusion coefficient from any grid cell to all cells remains constant at 1. According to the definition of maximum movement distance, it can be further constrained to have a sum of 1 from any grid cell to the adjacent cells and itself:

Based on the definition above, the updated existence probability of the target based on target motion prediction in $$ g_{mn}^T $$ after time $$ t_k $$ is denoted as $$ p_{mn}^-(t_{k+1}) $$, which is as follows:

After the update based on target motion prediction, the probability map will be updated again based on the real-time detection results of airborne detection sensors of every UAV. In this paper, Bayesian criteria are used to dynamically update the probability map^[8]. The update of existence probability after time $$ t_k $$ can be classified into the following three situations:

(ⅰ) i-th UAV finds a target in $$ g_{rc}^{Di}(t_{k+1}) $$ with $$ g_{mn}^T{\in}g_{rc}^{Di}(t_{k+1}) $$:

(ⅱ) i-th UAV finds nothing in $$ g_{rc}^{Di}(t_{k+1}) $$ with $$ g_{mn}^T{\in}g_{rc}^{Di}(t_{k+1}) $$:

(ⅲ) $$ g_{mn}^T $$ is outside the area of the detection range of any UAV:

3.2. Certainty map

A certainty map is a search information map that reflects the degree of UAV exploration of the mission area. Considering that UAVs move in mission grid cells, the region division of the certainty map is the same as $$ \Omega^U $$, and the certainty map varies with the UAVs' detectability and access situation. Denote the grid of a certainty map as $$ \Omega^H $$, in which the cell with m-th row and n-th column is denoted as $$ g_{mn}^H $$ and the corresponding certainty is denoted as $$ \chi_{mn} $$. The update of a certainty map after time $$ t_k $$ can be classified into the following two situations based on the access and the detection range of every UAV:

(ⅰ) i-th UAV flies past with the situation that $$ g_{mn}^H $$ has the same regional definition as $$ g_{rc}^{Di}(t_{k+1}) $$ in $$ \Omega^{Di}(t_{k+1}) $$:

As in Equation (14), the certainty of the grid cell increases with the detectability and access of UAVs.

(ⅱ) $$ g_{mn}^H $$ is outside the area of the detection range of any UAV:

where $$ \eta_H\in(0,1) $$ is the attenuation coefficient of certainty.

Detection response map

A detection response map is a search information map defined in this paper to provide feedback on the target detection results of UAVs, which solves the problem of the lack of direct search results in conventional search information maps. The region division of the detection response map is the same as $$ \Omega^U $$, and the detection response map varies with the UAVs' detection results. Denote the grid of detection response map as $$ \Omega^R $$, in which the cell with m-th row and n-th column is denoted as $$ g_{mn}^R $$ and the corresponding response value is denoted as $$ \psi_{mn} $$. The update of the detection response map after time $$ t_k $$ can be classified into the following three situations based on the detection results of every UAV:

(ⅰ) i-th UAV finds a target in $$ g_{rc}^{Di}(t_{k+1}) $$ with the situation that $$ g_{mn}^R $$ has the same regional definition as $$ g_{rc}^{Di}(t_{k+1}) $$ in $$ \Omega^{Di}(t_{k+1}) $$:

where $$ \eta_{R1}\in(0,1) $$ is the sensitivity coefficient of detection response. In order to avoid missed detections, the response value is sensitive to detected targets and insensitive to undetected targets. The increase in response value is related to the detection probability of UAVs.

(ⅱ) i-th UAV finds nothing in $$ g_{rc}^{Di}(t_{k+1}) $$ with the same situation above:

(ⅲ) $$ g_{mn}^R $$ is outside the area of the detection range of any UAV:

where $$ \eta_{R2}\in(0,1) $$ is the attenuation coefficient of response value.

4.1. MPC method

According to the idea of MPC, by predicting and evaluating the search behavior in a limited period of time under different search decisions, the control input required by the current UAV swarm can be determined based on the current environment information. Figure 7 shows the cooperative search decision-making process for multi-UAV systems based on the idea of MPC.

7. Cooperative search decision-making process for multi-UAV systems based on the idea of MPC.

It is assumed that an optimization process will consider the cooperative search process of $$ N_K $$ steps in future and evaluate the search benefits associated with different control inputs. Denote the search decision set of certain control inputs in future for multi-UAV systems at time $$ t_k $$ as $$ \mathbb{U}(t_k{\vert}t_k)=\{\mathbb{U}_1(t_k{\vert}t_k),$$$$\mathbb{U}_2(t_k{\vert}t_k),…,\mathbb{U}_i(t_k{\vert}t_k),…,\mathbb{U}_{N_U}(t_k{\vert}t_k)\} $$, $$ 1{\le}i{\le}N_U $$, in which $$ \mathbb{U}_i(t_k{\vert}t_k)=\{u_i(t_k{\vert}t_k),$$$$ u_i(t_{k+1}{\vert}t_k),…,u_i(t_{k+N_K-1}{\vert}t_k)\} $$ is the set of the control input sequence in future for i-th UAV and $$ u_i(t_{k+q}{\vert}t_k) $$ is the corresponding control input at time $$ t_{k+q} $$.

Obviously, the quantitative search effectiveness of the search path can be uniquely obtained through an objective function $$ \mathrm{J}(\cdot) $$, which is based on the state $$ X(t_k) $$ and search decision set $$ \mathbb{U}(t_k{\vert}t_k) $$ of the multi-UAV system at time $$ t_k $$, and the prediction evaluation process of the search can be achieved. Finally, the optimal search decision set $$ \mathbb{U}^*(t_k{\vert}t_k) $$ can be obtained by the optimization model of MPC as follows:

where $$ \mathbb{U}_i^*(t_k{\vert}t_k)=\{u_i^*(t_k{\vert}t_k),$$$$u_i^*(t_{k+1}{\vert}t_k),…,u_i^*(t_{k+N_K-1}{\vert}t_k)\} $$ is the optimal set of the control input sequence in future for i-th UAV, and $$ u_i^*(t_{k+q}{\vert}t_k) $$ is the corresponding optimal control input at time $$ t_{k+q} $$. Finally, by using the first term of the optimal decision set as the system optimal control input $$ U^*(t_k{\vert}t_k)=\{u_1^*(t_k{\vert}t_k),u_2^*(t_k{\vert}t_k),…,u_{N_U}^*(t_k{\vert}t_k)\} $$ at time $$ t_k $$ and guiding the UAV swarm to perform a one-step search, the cooperative search decision-making process for the current step can be completed.

4.2. Objective function

Considering that no actual search actions occurred during the estimation of the search results in prediction, the virtual search process of future steps is based on the premise that no targets are found in each step.

5.1. Priority-encoded IGAFA algorithm

5.2. Complete search process

Based on the previous sections, the procedures of a complete search process for moving targets using CSMTPE, are as shown in Table 2.

6.1. Motion prediction of moving targets

Assuming that there are two targets to be searched in the mission area, the initial position of the target is subject to the two-dimensional normal distribution, where the means of position are $$ [2500 \quad 1500]^T $$ and $$ [2500 \quad 2600]^T $$, and the variances of position are $$ [160000 \quad 0; \quad 0 \quad 160000] $$ and $$ [160000 \quad -128000; \quad -128000 \quad 160000] $$.

To test the update process of the probability map, three experimental scenarios were used: (Ⅰ) target prediction without the search interference of UAVs, (Ⅱ) under the search interference of static UAVs, and (Ⅲ) moving UAVs. Figure 8 shows the updating process of the probability map in different scenarios at simulation steps k = 0, k = 10 and k = 20. In scenario Ⅰ, as the targets are not affected by the UAV search, the target follows a Gauss-Markov motion model. The result shows that the existence probability of the target spreads towards adjacent grid cells, and the uncertainty of the target position increases with the simulation steps. In Scenario Ⅱ, the target will be subject to the evasive virtual force A generated by the UAVs in addition to random movement, showing a trend of evading the UAVs. The result shows that there is a significant reduction in the existence probability of the target around the UAVs. In Scenario Ⅲ, the target is affected by the evasive virtual forces A and B, and the result shows that the target in the UAVs' forward direction will additionally evade both sides of the forward direction, while the target behind the UAVs may return to the previous search area by random movement after the UAVs have moved away.

Figure 8. Updating process of the probability map in different scenarios.

Compared with previous methods^[10,11], the changes in the probability map in the simulation results are consistent with the actual experimental scenario rather than purely numerical assumptions, reflecting the reasonable evasive maneuvers of moving targets and improving the efficiency of the search process.

6.2. Analysis on the influence of different factors

In addition to the configuration in the first experiment, set $$ N_U=3 $$ with the initial position of $$ [850 \quad 850]^T $$, $$ [1850 \quad 3550]^T $$ and $$ [1250 \quad 2250]^T $$. To avoid a lengthy coefficient tuning process, a set of coefficients that have performed well in multiple experiments is directly presented here, and the factor tunning of the objective function and IGAFA is mainly discussed in the following. For the coefficients in the search information map, set $$ N_K=20 $$, $$ \eta_H=0.95 $$, $$ \eta_{R1}=0.67 $$, $$ \eta_{R2}=0.9 $$ as a correspondence to the update method of search information map mentioned earlier. For the coefficients in the basic parameters in IGAFA, set $$ P_{c1}=0.4 $$, $$ P_{c2}=0.3 $$ and $$ P_m=0.3 $$ as the basic configuration of the following experiments.

In this experiment, the coefficients in the objective function and the factors in IGAFA will be determined. Firstly, considering the sensitivity to various indices and the balance between benefits and penalties, set $$ \lambda_E=1.0 $$ as the basis of the main search benefit and set $$ \lambda_C=0.1 $$ as the basis of the penalty for UAV motion cost. To test the influence of the configuration with different ratios of $$ \lambda_H $$ and $$ \lambda_R $$, we use a traditional GA as the basic algorithm to search the moving target with the same probability distribution of the initial position as that in the previous experiment. Table 3 shows the search result of the number of targets found at simulation step $$ k=100 $$.

Table 3

Average number of targets found with different search coefficients

	$$ \lambda_R=0.2 $$	$$ \lambda_R=0.5 $$	$$ \lambda_R=0.8 $$
$$ \lambda_H=0.2 $$	1.15	1.48	1.62
$$ \lambda_H=0.5 $$	1.28	1.55	1.68
$$ \lambda_H=0.8 $$	1.25	1.53	1.58

Table 3 indicates that the group of $$ \lambda_H=0.5 $$ and $$ \lambda_R=0.8 $$ has the best optimization performance. Therefore, the coefficients are set accordingly for better detection-evasion effectiveness of adversarial games in the search process.

Secondly, the factors in IGAFA will be determined by calculating the optimal search path for UAVs at their initial position. First of all, without considering the weight matrix during the mutation process, test the influence of different $$ A_{IT} $$ and $$ B_{IT} $$ on the result of optimization. Figure 9 shows the best fitness under different iteration steps with different $$ A_{IT} $$ and $$ B_{IT} $$. It indicates that the group of $$ A_{IT}=1 $$ and $$ B_{IT}=1 $$ has the best optimization performance. Therefore, the factors are set accordingly for a better evolution process of the probability distribution of selected genes.

Figure 9. Best fitness under different iteration steps with different $$ A_{IT} $$ and $$ B_{IT} $$.

Next, considering the weight matrix during the mutation process, we test the influence of different $$ \eta_W $$ on the result of optimization. Figure 10 shows the best fitness under different iteration steps with different $$ \eta_W $$. It indicates that $$ \eta_W=0.5 $$ has the best optimization performance. Therefore, the factor is set accordingly for a better process of mutation.

Figure 10. Best fitness under different iteration steps with different $$ \eta_W $$

6.3. Process and result analysis of a complete search simulation

Based on the configuration in the second experiment, randomly arrange moving targets with pre-set position distribution and simulate the complete cooperative search process using CSMTPE with the procedures in Table 2. Figure 11 shows the search process of a complete search simulation, including the current search path, the search path for the next $$ N_K $$ steps, and the probability map, certainty map, and detection response map at the simulation steps $$ k=1 $$, $$ k=19 $$, and $$ k=40 $$, where UAV 1 finds and locks on Target 1 at step $$ k=19 $$ and UAV 2 finds and locks on Target 2 at step $$ k=40 $$.

Figure 11. Search process under different simulation steps.

In the probability map in Figure 11(a), the existence probability of the target dynamically changes with the search process, and the UAV swarm converges and aggregates the probability to the actual position of the targets, achieving an effective search process. The certainty map in Figure 11(b) reflects the currently searched area of the UAV swarm, avoiding repeated searches of the same area in a short time, while the lower certainty after a long time will guide the UAVs to revisit the area to enhance the search order in the later stage of the search. The detection response map in Figure 11(c) reflects the current target detection situation of the UAV swarm, which includes both the real alarm signal caused by the target and the false alarm signal of the sensor, and finally guides the UAVs to focus on searching the alarm area. Figure 11(a) also reflects that the search path guided by the CSMTPE can effectively explore the area with the initial existence probability of the target, and the priority of accessing the area with higher probabilities is greater. Meanwhile, the allocation of search areas among UAVs is also reasonable.

Furthermore, the above initial conditions are used for multiple simulation experiments, where the initial positions of the targets are randomly set differently in each simulation. The target search results within the maximum simulation step of 100 are recorded. Figure 12 shows the result of 40 Monte Carlo experiments of the search for two moving targets. The result indicates that the number of targets found gradually increases with the simulation steps, with an average of 1.83 targets found and a search completion rate of 92% within 100 simulation steps.

Figure 12. Numbers of targets found under different simulation steps.

In order to simulate more search scenarios, five cases of configuration were designed here, with the number of UAVs and targets being 3:1, 3:2, 6:2, 6:3, and 6:4, respectively. Figure 13 shows the search process in different cases of UAV-target configuration. In each case, the background of each figure shows the probability distribution of the targets' initial position, while the positions of UAVs and moving targets are plotted at the moment when all targets were discovered. The result of the average number of targets found is shown in Table 4.

Figure 13. Search process in different cases of UAV-target configuration.

The result in Table 4 indicates that the larger the ratio of the target number to the UAV number, the higher the corresponding search completion rate, and the faster the target can be found, and it also proves that the CSMTPE proposed in this paper can effectively cope with the multiple target search in different target distribution situations.

In order to verify the effectiveness of the update method of probability map based on target motion prediction, the search result was tested under the same simulation conditions without using Equation (10) for the update of the probability map. Based on the search result in Table 4, Figure 14 shows an intuitive comparison of the search result with and without the motion prediction, and it clearly indicates that when target motion prediction is not used, the average number of targets found will significantly decrease, which is because probability map cannot effectively predict the random and evasive movement of the target. Moreover, when the UAV approaches the target, the existence probability of the target shifts between adjacent grid cells without motion prediction due to the rapid movement of the target, and the probability of losing the target will be significantly increased.

Figure 14. Comparison of the search results with and without the motion prediction.

6.4. Comparison of different optimization algorithms

In order to compare the optimization effectiveness of the proposed optimization algorithm, a set of optimization algorithms will be used as a comparison in the following experiment. First of all, a motion-encoded particle swarm optimization^[10] (ME-PSO) and a motion-encoded genetic algorithm^[15] (ME-GA) are chosen as the comparison group of the direct motion-encoding method as these are commonly used to solve similar problems in recent research. Secondly, to test the proposed priority-encoding method, these optimization algorithms have been rewritten with corresponding encoding methods and denoted as PE-PSO and PE-GA, respectively. In addition, to compare with more optimization algorithms, a differential evolution (DE) algorithm^[28] and an artificial bee colony (ABC) algorithm^[29] are chosen and rewritten with priority-encoding methods, which are PE-DE and PE-ABC.

In order to match the parameters used by IGAFA in the second experiment, the parameters in GA are set to $$ P_c=0.7 $$ and $$ P_m=0.3 $$, and the PSO, DE and ABC have the same population size as IGAFA. Figure 15 shows the best fitness under different iteration steps using different optimization algorithms. The result shows that the fitness optimized by ME-PSO or ME-GA is significantly higher than the fitness using the indirect priority-encoding method proposed in this paper, which indicates that the direct motion-encoding method is not appropriate when dealing with the proposed optimization problem. The fundamental reason is that the direct motion-encoding method cannot guarantee the legitimacy of its descendants, and it is difficult to achieve convergence of the optimal solution during the optimization. Each bit in the encoding has a significant influence on the results, and the optimization is inefficient under a large solution space.

Figure 15. Best fitness under different iteration steps using different optimization algorithms.

Among the algorithms that use the priority-encoding method, traditional PE-GA has the slowest optimization speed. By comparison, PE-PSO, PE-DE and PE-ABC have faster optimization speeds in the early stage due to their characteristics being close to random-search. However, in the later stage, considering the continuity problem of the solution, the optimization speed significantly decreases. As an improved algorithm for GA, priority-encoded IGAFA optimizes the evolution direction while retaining the excellent genes of its parents, making the evolution process closer to the encoding characteristics, thus maintaining good optimization efficiency in the middle and later stages of the optimization process. It proves that the fine-adjustment mechanism proposed in this paper has efficient optimization performance when dealing with priority-encoding methods.

In order to further compare the influence on the target search process when using different optimization algorithms, PE-GA and ME-PSO are selected as typical optimization algorithms representing priority-encoding method and motion-encoding method, respectively.

Table 5 shows the average execution time with different $$ N_K $$ and $$ N_U $$. The result indicates that ME-PSO has a shorter execution time because of the defect of illegal descendant, and IGAFA has a longer execution time because of the added improvement methods, while the execution time of the three algorithms is roughly the same.

We apply different algorithms to the search scenarios set in the third experiment and simulate the cooperative search using CSMTPE. The results of the average number of targets found using different optimization algorithms at simulation step $$ k=100 $$ are shown in Table 6.

Table 6 indicates that the ME-PSO has a significantly smaller number of targets found compared to other algorithms, while the search result using PE-GA is slightly worse than that using IGAFA. This is consistent with the conclusion obtained in Figure 15, proving that IGAFA is more suitable for solving such problems.

Comparison of different search methods

We now compare the CSMTPE proposed in this paper with traditional greedy search method and parallel search method. The greedy search method calculates each grid cell within the motion range of each UAV in each simulation step according to Equation (10) with $$ N_K=1 $$, and selects the cell with the minimum objective function value as the cell to be searched in the next simulation step. The parallel search method refers to the parallel search of targets by each UAV in the mission area, achieving a coverage search. Table 7 shows the search results using different search methods at simulation step $$ k=100 $$ under the same configuration in the third experiment, and Figure 16 shows the combination result of targets found at simulation steps $$ k=25 $$, $$ k=50 $$, and $$ k=100 $$.

Figure 16. Comparison of the search results using different search methods.

Table 7 indicates that CSMTPE has more targets in comparing search methods. Parallel search does not consider the initial distribution information of the target, and the search range is too rough, resulting in poor search efficiency. Greedy search has a similar search rate to CSMTPE in Case I and Case III, but overall, the number of targets found is lower than CSMTPE.

As shown in Figure 16, when using the greedy search method, the number of targets found is significantly less than that of CSMTPE at simulation step $$ k=25 $$, reflecting that it wasted more search steps due to blind search in the early stage of the search, while the number of targets found in parallel search method is related to the target distribution, resulting the poor search efficiency. In contrast, the CSMTPE proposed in this paper continuously updates various search information maps to obtain more target information and takes into account both short-term and long-term search efficiency when planning search paths.

Considering the existence of situations where the prior information of the target is unknown, the original problem will degenerate into a covering search problem. Under this circumstance, another similar experiment is conducted, and the search results are shown in Table 8.

When the prior information of the target is unknown, the search results reflect the conclusion that the search efficiency of CSMTPE will degenerate into the level of that using the parallel search method. Meanwhile, it also shows that the greedy search method performs significantly worse compared to other search methods above due to the lack of long-term predictability. In summary, it proves that CSMTPE can adaptively and effectively solve the multi-UAV regional cooperative search problem.

Authors' contributions

Made substantial contributions to the conception and design of the study and performed simulation and interpretation: Wang Z

Provided administrative, technical, and material support: Zou W, Li S

Availability of data and materials

Not applicable.

Financial support and sponsorship

None.

Conflicts of interest

All authors declared that there are no conflicts of interest.

Ethical approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Copyright

© The Author(s) 2023.