# Minimum Time Search in Real-World Scenarios Using Multiple UAVs with Onboard Orientable Cameras

Authors: Sara Pérez-Carabaza; Eva Besada-Portas; José Antonio Lopez-Orozco; Gonzalo Pajares
Journal: Journal of Sensors
 (2019)

---

## Abstract

This paper proposes a new evolutionary planner to determine the trajectories of several Unmanned Aerial Vehicles (UAVs) and the scan direction of their cameras for minimizing the expected detection time of a nondeterministically moving target of uncertain initial location. To achieve this, the planner can reorient the UAVs cameras and modify the UAVs heading, speed, and height with the purpose of making the UAV reach and the camera observe faster the areas with high probability of target presence. Besides, the planner uses a digital elevation model of the search region to capture its influence on the camera likelihood (changing the footprint dimensions and the probability of detection) and to help the operator to construct the initial belief of target presence and target motion model. The planner also lets the operator include intelligence information in the initial target belief and motion model, in order to let him/her model real-world scenarios systematically. All these characteristics let the planner adapt the UAV trajectories and sensor poses to the requirements of minimum time search operations over real-world scenarios, as the results of the paper, obtained over 3 scenarios built with the modeling aid-tools of the planner, show.


---

## Body

## 1. Introduction

The strong research interest in UAVs trajectory planning takes advantage of the UAVs capabilities to perform different types of military and civil missions, such as georeference [1], wildlife monitoring [2], or target tracking [3]. Moreover, the current developments of their onboard sensorial systems make them also ideal for performing risky long-endurance reconnaissance and surveillance operations. This work focuses on a type of Probabilistic Target Search Problem (PTSP), named minimum time search (MTS) [4], since reducing the time required by the UAVs and their onboard sensors to detect the target is a critical objective of the mission. The selected problem has several applications that include looking for survivors after natural disasters (e.g., after fires or earthquakes), search and rescue operations, or searching for military targets.Approaches that tackle PTSPs determine the trajectories of the UAVs in spite of the uncertainty associated with the target location and sensor capabilities. To do this, they probabilistically model the different uncertainty sources. In more detail, the information about the target position is modeled with an initial probability map (with the belief of target presence within the search area) and a probabilistic motion model, while the sensors uncertainty is modeled with detection likelihood functions. As an example, Figure1 schematizes a PTSP where two UAVs equipped with electrooptic sensors (cameras) look for a static lost off-road vehicle in the mountains. The colored map at the bottom shows the altitude of the search area (with green being associated with valleys and brown with mountains), the grey shadowed map in the middle displays the target initial belief (darker areas, associated with lower altitudes, indicate a higher probability of target presence), and the blue/red polygons, respectively, show the current observed areas (within the cameras footprints) by the onboard sensor of each UAV. When the PTSP is also a MTS operation, the regions with higher target probabilities (darker grey areas) should be observed (falling within the cameras footprints) as soon as possible in order to minimize the target detection time.Figure 1
Schematic representation of the problem.One of the main goals of MTS planners is to reduce the target detection time, which can be achieved by optimizing the expected time of target detection [5–10]. Other PTSP approaches optimize alternative criteria, such as maximizing the probability of target detection [11–14] or minimizing its counterpart probability of nondetection [15, 16], maximizing the information gain [17], minimizing the system entropy [18], minimizing its uncertainty (areas with intermediate belief of target presence) [19], or optimizing normalized or discounted versions of the previous criteria [4, 20–22]. A common characteristic of the different approaches is that, in scenarios with bounded resources (e.g., limited flying time or fuel), they often obtain better results than predefined search patterns (e.g., spiral, lawnmower), as they adapt the UAV trajectories to the scenario specific target initial belief and motion [6, 20]. Besides, although the approaches that optimize the previous PTSP criteria can share the same elements and probabilistic models, MTS distinctiveness is the extreme influence of the visiting order of the high probability regions in the expected time of target detection, as prioritizing flying over high probability areas first increases the chances of finding the target earlier [4].The NP-hard complexity of PTSP [23] is tackled with suboptimal algorithms and heuristics, such as gradient-based approaches [13, 15–17, 19], greedy methods [8, 12, 20], cross-entropy optimization [4, 7], Bayesian optimization algorithms [5], ant colony optimization [6, 9], or genetic algorithms [10]. Besides, streamlined formulations of the problem are typically accepted in order to further simplify the problem complexity. They range from considering static targets [8, 10, 13, 15, 16, 19, 20] instead of dynamic ones [4–7, 9, 11, 12, 14, 17, 21, 22] to modeling the sensors ideally [4–6, 8] instead of realistically (e.g., as radars [9–12] or downward-looking cameras [13, 14, 17, 19]) or to assuming that the UAVs fly following straight-lines according to the eight cardinal directions [4–8] or optimized waypoints [17, 21, 22] instead of considering the physical maneuverability constraints induced by the dynamic models of the UAVs [9–16, 19, 20]. Additionally, in some cases (e.g., in [10, 13, 16, 17, 19, 20]) the approach uses a receding horizon controller to divide the UAVs trajectory into sequentially optimized sections, narrowing the optimization search space at the expense of constructing suboptimal myopic solutions. Finally, it is worth noting that the approaches are often tested over synthetic scenarios, built by the authors, without a clear relation with a real-world problem.The previous simplifications and the lack of analysis over real-world scenarios reduce the applicability of the approaches to real-world problems. This is especially relevant within the subset of MTS approaches [4–10], as the majority of them have been developed for UAVs flying accordingly to the eight cardinal directions [4–8] or only tested over ideal sensors and synthetic unreal scenarios [4–6]. To pave the path of using MTS methods for real-world scenarios, this work extends the capabilities of the planner introduced in [10] by contributing to the following fields:(i)
MTS mission definition, by combining intelligence and geographical information to construct the target probability models that will be used in real-world missions. This feature is inspired by software tools that use geographical information to build the target initial belief or motion model [24–26], monitor search missions [26, 27], or evaluate predefined search patterns [24, 28]. Its main benefit, as the results of this paper will show, will be the substitution of the synthetic scenarios used in previous MTS planners by real-world inspired ones.(ii)
MTS planning, by optimizing simultaneously the UAV trajectory (by means of changing the UAV heading, speed, and height) and camera pose (azimuth and elevation). This new competence combines the UAV trajectory optimization capabilities of previous MTS planners (which usually manipulate only the UAV heading) with the sensor moving capabilities of only a few PTSP approaches [21, 22]. As this innovation supports a higher moving capability of the camera footprint and a quicker coverage of the high probability areas, it is especially relevant for MTS where the target has to be detected as soon as possible.(iii)
Sensor characterization, by incrementing the realism of the camera detection behavior, modeling the effects of the terrain elevation (occlusions and target-camera distance variation), camera orientation, and target and sensor size in the footprint dimensions and sensor likelihood. Not only is the realism of the likelihood model crucial to shorten the differences between the simulated and real behavior of the camera, but it will also be regarded during the mission definition presented in the Results section in order to set up the scenarios correctly.In short, this paper presents a new planner that optimizes the UAV trajectories and onboard cameras orientations in real-world dynamic MTS scenarios, where the behavior of a dynamic target is probabilistically modeled from a novel perspective by combining intelligence and geographic information, and whose new camera likelihood model takes into account the effects of the terrain elevation in the camera detection capabilities. Besides, we extend the capabilities of the MTS planner presented in [10], which only manipulated the trajectory heading of a set of UAVs equipped with fixed radars in order to optimize the detection time of static targets in synthetically built scenarios, with the simultaneous optimization of UAVs trajectories and sensor poses over realistically built scenarios with dynamic targets. Moreover, the new planner also considers the extension to heterogeneous UAVs, which can have different parametrizations and start and leave the search mission at different times. Finally, this paper describes in detail, and in an algorithmic form, the functionality of the new planner with the purpose of clarifying the interaction between the different elements during the optimization process.The remaining of this paper is organized as follows. The second section introduces the probabilistic formulation of the MTS problem, describes the novel approaches used to model the target initial state and motion behavior, and presents the new models used for the UAVs and their cameras. The third section details the new multistepped MTS planner presented in this paper, introducing its steps in an algorithmic form and analyzing its computational complexity. The fourth section compares the state-of-the-art of the closest related work with our new planner and highlights its differences with the previous planner in [10]. The fifth section analyzes its performance over three different real search scenarios and shows the benefits of letting the planner reorient the camera. And in the ending section, the conclusions are drawn and some open research questions discussed.


## 2. Minimum Time Search Definition

This section presents the probabilistic formulation of the MTS problem and describes the approaches used in this work to model the target, UAV, and sensor behavior.


### 2.1. General Problem Formulation

In the MTS problem presented in this paper, there is a set ofU UAVs overflying a search area Ω (discretized into a grid GΩ of wx·wy square cells) in order to detect the presence of a single target located within it. Besides, due to the uncertainty associated with the position and dynamics of the target and to the measurements of the onboard sensor of each UAV, the MTS problem is formulated within the probabilistic framework introduced in [5].In more detail, the information about the target initial position (represented by random variableτ0) is modeled with the initial target belief b(τ0)≜P(τ0), which is a probability mass function that represents the chances of finding the target at each τ0∈GΩ. Besides, when the target is moving, the information about its dynamics is described with the target Markovian motion model Pτt∣τt-Tτ, which states the probability that the target moves from any cell τt-Tτ∈GΩ to any cell τt∈GΩ during a time lapse Tτ. Moreover, the detection capabilities of the onboard sensor of each UAV are described with the likelihood function Pzut∣τt,sut, which returns the probability of detecting (zut=D) or not detecting (zut=D¯) a target placed at τt from the UAV and sensor pose stored in the state variable sut of the u-th UAV.For those readers familiarized with the Recursive Bayesian Filter (RBF, [12]), the previous probability models b(τ0), Pτt∣τt-Tτ, and Pzut∣τt,sut are used to obtain the target belief b(τt)≜Pτt∣z1:U0:t,s1:U0:t, which represents the chances of finding the target at time step t and at each cell τt∈GΩ, given the trajectory of the UAVs and sensors poses s1:U0:t and the sensor measurements z1:U0:t. The RBF process, stated in (1), calculates the current belief b(τt) from the previous time step belief b(τt-Tτ)≜Pτt∣z1:U0:t-Tτ,s1:U0:t-Tτ by (1) incorporating the current UAVs location and sensor poses s1:Ut and the last sensor measurements z1:Ut with the likelihood functions Pzut∣τt,sut, and by (2) displacing and redistributing the target location probability with the target motion model Pτt∣τt-Tτ. Besides, (1) assumes that the sensors provide measurements only at time steps that are multiple of Tτ and its normalization factor γ is used to ensure that ∑τt∈GΩb(τt)=1.(1)bτt=γ∏u=1:UPzut∣τt,sut∑τt-Tτ∈GΩPτt∣τt-Tτbτt-TτA MTS mission can be formulated as an optimization problem, where some criteria related to the previous probabilistic models and other mission objectives are evaluated in order to determine the best UAV trajectories and sensor posess1:U0:tend during the duration of the mission tend. A useful criterion for MTS is the Expected Time of Detection (ETD, [5–7, 9, 10]), which measures the average time of detection of the target when the UAVs and the sensors follow the trajectories and poses defined in a given s1:U0:tend. It is calculable with (2)-(4), setting b~(τt) to b(τ0) for the initial case (t=0). The first equation is similar to RBF equation (1), but as it lacks the normalization term γ and all measurements z1:U0:t are nondetection, it obtains the “unnormalized belief” b~(τt) assuming that the onboard sensors do not detect the target from the UAVs trajectories and sensor poses in s1:U0:t. The second one obtains Pnd(s1:U0:t)≜PD¯1:U0:t∣s1:U0:t, which is the probability of not detecting the target from s1:U0:t and whose value decays as the sensors make nondetection observations in regions with probability of target presence. The third expression obtains the ETD from s1:U0:tend when Pnd(s1:U0:t) becomes 0 for some t≤tend and underestimates its value when Pnd(s1:U0:tend)>0. Finally, it is worth noting that minimizing ETD(s1:U0:tend) is a better option for MTS than maximizing the probability of detection along the whole trajectory Pd(s1:U0:tend)=1-Pnd(s1:U0:tend), as the MTS objective is to collect as much probability as possible sooner [4–6].(2)b~τt=∏u=1:UPzut=D¯∣τt,sut∑τt-Tτ∈GΩPτt∣τt-Tτb~τt-Tτ(3)Pnds1:U0:t=∑τt∈GΩb~τt(4)ETDs1:U0:tend=∑t=0tendPnds1:U0:t·TτBesides, in order ensure that the UAV trajectories and sensor posess1:U0:tend calculated by our MTS planner are feasible from the maneuverability point of view, we exploit the UAV and sensor deterministic dynamic model of each UAV-sensor pair s˙ut=f(sut,aut), where aut stands for the set of control actions (e.g., commanded UAV heading, height & speed, and sensor elevation & azimuth) and s˙ut for the time derivative of the UAV location and sensor poses. In particular, our MTS planner uses this dynamic model to obtain the solutions (best UAV trajectories and sensor poses s1:U0:tend) from the initial UAVs locations s1:U0 and the sequence of sets of control actions a1:U0:tend proposed, manipulated, and evaluated by our approach in order to optimize MTS missions. Finally, it is worth highlighting that due to the deterministic nature of the UAV and sensor dynamic model, there is an unambiguous relation between the best sequence of control actions and the best UAVs trajectories and sensor poses.The realism of the four models (i.e., of the initial target beliefb(τ0), of the target motion model Pτt∣τt-Tτ, of the sensor likelihood Pzut∣τt,sut, and of the UAV and sensor motion model s˙ut=f(sut,aut)) is crucial to avoid discordances between the real ETD of a given s1:U0:t and the ETD calculated by our approach. Hence, in the rest of this section we present the new models proposed in this paper to bring realism to the definition and solution of MTS missions performed by fixed-wing UAVs equipped with orientable cameras.


### 2.2. Initial Target Belief Definition

To construct the initial target beliefb(τ0), we merge knowledge coming from different sources, using a different probability layer bl(τ0) for each information source and performing with (5) an addition of the L probability layers weighted with their reliability/importance coefficients wl. In other words, considering the first term in (5) is a normalization coefficient that ensures that b(τ0) is a probability function (i.e., ∑τ0∈GΩb(τ0)=1), our initial target belief b(τ0) is calculated as the mixture of the beliefs bl(τ0) associated with the different l-th information sources.(5)bτ0=1∑l=1Lwl∑l=1Lwlblτ0The layers can be associated with geographical information (e.g., terrain altitude, road maps) or intelligence/user-defined clues (e.g., last/habitual location areas of the target). For the examples of this work, we consider the following two layers:(1)
The terrain elevation probability layerb1(τ0), obtained with the following steps. First the digital elevation model (DEM, [29]) of Ω is automatically resampled to the size of the cells in GΩ in order to obtain the average height of each cell in GΩ. Next the user/operator is required to divide the existing elevations within the cells in GΩ into consecutive ranges and to assign a chance of target presence to each range. Finally, the method automatically determines the cells in each elevation range and the probability associated with all τ0∈GΩ, distributing the chances of target presence assigned by the operator. As this way of proceeding generates a geographical probability layer b1(τ0) where areas with similar altitude share the same initial belief, it automatically spreads uniformly the belief over different regions of the search area.(2)
The intelligence probability layerb2(τ0). For this layer the operator has to graphically define a mixture (weighted addition) of Gaussians (centered in eligible locations of the search area Ω and spread according to selectable standard deviations) and of polygonal areas (defined by their external points placed in the desired locations of the search area Ω) with uniform probabilities (assigned by the operator). The weights of each element (each Gaussian and/or each polygonal area) within the intelligence probability layer b2(τ0) are also selected by the operator according to the information gathered about the last known location of the target.As an illustrative example, we show the initial belief defined by an operator when looking for a drifting boat next to the cost in Areia Branca, Brazil. To obtain the probability elevation layerb1(τ0) shown at the bottom of Figure 2(b), where darker/lighter greys are associated with higher/lower probability cells, the operator analyzes the elevation map represented in Figure 2(a) and assigns high chances of target presence to the elevations by the coast (represented in dark green in the elevation map) to model a strong belief that the boat may have arrived a ground when the search mission starts, lower chances to the sea-level elevation (in blue) to model a moderate belief that the target is still navigating, and zero chances to higher inland elevations. In order to define the intelligence probability layer b2(τ0), represented at the top of Figure 2(b), the operator centers a moderate spread Gaussian in the last known position of the boat. Finally, the operator assigns a weight for each layer (w1=2 and w2=1) in order to produce the initial target belief b(τ0) shown in Figure 2(c).Figure 2
Initial target belief definition used for a drifting boat (in Areia Branca, Brazil).
(a)
Elevation map (m) (b)
Probability layers (b1(τ0) at bottom &  b2(τ0) at top) (c)
Initial beliefb(τ0)


### 2.3. Target Motion Model Definition

To define that target motion model, we can distinguish two cases:(1)
Scenarios with static targets. In this case, and due to the immobility of the target,Pτt∣τt-Tτ=1 if and only if τt and τt-Tτ refer to the same cell in GΩ, and Pτt∣τt-Tτ=0 otherwise. Besides, in this case, the equations that contain the term ∑τt-Tτ∈GΩPτt∣τt-Tτb(τt-Tτ) can be simplified as ∑τt-Tτ∈GΩPτt∣τt-Tτb(τt-Tτ)=b(τt-Tτ).(2)
Scenarios with dynamic targets. To construct this type of target motion model, we combine different types of information (e.g., elevation data and sea currents). In this work, the operator has to provide the following:
(i)
The elevation rangeER where the target is allowed to move within Ω. This allows the operator to indirectly assign the probabilities corresponding to the static target behavior to those cells that do not belong to the ER. In other words, ∀τt-Tτ∉moving(GΩ,ER), the target motion definition process automatically makes Pτt∣τt-Tτ=1 if τt=τt-Tτ and Pτt∣τt-Tτ=0 otherwise.(ii)
K vectors v1:K, each with 9 values that represent, for a few selected cells c1:K⊂GΩ, the chances of moving to their 8 neighbor cells N(ck) and of staying at the same cell ck. Besides, in order to force the target to stay within the search area Ω, some of the 9 values of the vectors vk that correspond to cells in the borders of GΩ are automatically set to zero. With this information, the target motion definition process makes Pτt∣τt-Tτ=0 for all τt∉N(τt-Tτ)∪τt-Tτ and computes Pτt∣τt-Tτ applying the potential field method (according to the cardinal directions that make the target move from one cell to its neighbors) for all τt∈N(τt-Tτ)∪τt-Tτ and τt-Tτ∈moving(GΩ,ER). In more detail, if adj(τt-Tτ,d) with d∈1:8 stands for the adjacent cell to τt-Tτ in the d-th cardinal direction, adj(τt-Tτ,9)=τt-Tτ represents stay in the same cell, dist(τt,ck) is the distance between cell τt and cell ck, and vk[d] is the d-th element of vk, then Pτt=adjτt-Tτ,d∣τt-Tτ∝∑k=1:Kvk[d]/dist(τt,ck) and ∑d∈1:9Pτt=adjτt-Tτ,d∣τt-Tτ=1. This way of proceeding allows the operator to define the distribution of the target probability around the neighborhood cells in only a few cells c1:K⊂GΩ and extend this definition to the remaining moving cells taking into account the distance between the selected cells and the others.(iii)
The time periodTτ that has to be used to apply Pτt∣τt-Tτ to the belief b(τt-Tτ). As the speed of the target is related to the quotient of the size of the cell and Tτ, this value is used by the planner to relate the simulation of the motion of the target to the simulation of the motion of the UAVs and sensors.As an illustrative example, we also obtain the target motion model of the boat drifting example next to the cost in Areia Branca, Brazil. In this case, the operator definesvk for three cells ck∈GΩ in the sea, making vk[9]=0.1, vk[d]=2 for one of the cardinal directions d towards the coast, and vk[r]=0 for the remaining directions r∈1:8-d. Figure 3 represents with the starting point of the red arrows the location of the selected c1:3 and with the arrows themselves the direction d where vk[d]=2 for each ck. Besides, the operator also makes the moving ER≤0 to only let the target (i.e., the drifting boat) move in the sea. To summarize the result of the target motion definition process, the green arrows in Figure 3 represent the direction obtained by weighting, for each cell τt-Tτ∈moving(GΩ,ER), the arrows in the 8-cardinal directions d with Pτt=adjτt-Tτ,d∣τt-Tτ. Besides, the lack of inland green arrows indicates the static behavior of the target in this area. In short, the green arrows show how the defined target motion makes the belief over the sea move towards the coast and remain unchanged inland. Finally, it is worth highlighting that, thanks to our target motion model definition approach, this complex model is automatically built just considering the elevation range (ER≤0) and the target movement probabilities in three cells of GΩ.Figure 3
Target motion definition used for a drifting boat (in Areia Branca, Brazil).


### 2.4. Camera Likelihood with Terrain Occlusion

The likelihood of detecting the target with the onboard camera of each UAV in a given cellτt∈Ω can be calculated scaling the Target Task Performance (TTP) metric [30] with the percentage of the selected cell within the camera footprint. This scaling, used to reduce the detection probability in the cells of the footprint border, can be modeled with (6), where τt∩footprintsut represents the area (size of the surface) of the common region of the cell τt and of the footprint of the camera (oriented and placed according to sut), τt is the total area of cell τt, and TTPF(τt,sut) is the target task performance function.(6)PDut∣τt,sut=τt∩footprintsutτtTTPFτt,sutTo determine the footprint of the camera from the UAV location and camera pose insut, we follow the process schematized in Figure 4(a) and consisting of the following two steps. First, we calculate geometrically the camera footprint at sea level taking into account the following pieces of information within sut: UAV location (xut,yut,hut) and camera azimuth, elevation, and FOV (αut,ηut,FOV). To do this, we consider that the sensor location is the same as the UAV location (since its deviation from the location of the UAV mass center is negligible for the mission) and that its orientations are measured with respect to the vehicle coordinate frame (since the camera gimbal compensates the UAV attitude). Second, we approximate the corners of the real camera footprint by the intersections (obtained using the 3D Bresenham algorithm [31]) of the terrain elevation with the 4 lines that join the camera with the 4 corners of the sea-level camera footprint.Figure 4
Camera likelihood model.
(a)
Footprint determination (b)
TTPF curve vs.VT (c)
TTPF curves vs. relative height (forStarget=5m)The value of the target task performance function TTPF for each cellτt within the sensor footprint from the UAV location and sensor pose in sut is calculated with (7), where VT(τt,sut) is the cycle size of the target and V50=2.03 (according to [32, 33]) is the critical cycle size for detecting a small target with a likelihood of 50% (since when VT=V50, TTPF=1/(1+1)=0.5). Moreover, the target cycle size is calculated with (8), where Starget is the real size of the target, β(τut) is the angle between the along-scan (vertical) and cross-scan (horizontal) directions of the footprint, and GSDi (with ∈vertical,horizontal) is the ground sample distance (corresponding size of a pixel of the camera over terrain, obtainable as the ratio of the length of the footprint in each direction to its corresponding number of pixels). Taking into account the trigonometric relations between the UAV location and camera footprint, the angle β(sut) can be obtained solving (9) and the ground sample distance in each scan direction can be calculated with (10) and (11), where hterrain(τt) is the height of the terrain at the target location τt and #pixeli is the number of pixels of the camera in the vertical or horizontal directions (i∈v,h).(7)TTPFτt,sut=VTτt,sut/V501.51+0.24VTτt,sut/V501+VTτt,sut/V501.51+0.24VTτt,sut/V50(8)VTτt,sut=StargetGSDvτt,sut·GSDhτt,sut·sin⁡βsut(9)tgβsut=tg-1ηut-FOV/2-tg-1ηut+FOV/2tgFOV/2sin-1ηut-FOV/2-sin-1ηut+FOV/2(10)GSDvτt,sut=hut-hterrainτt#pixelvtg-1ηut-FOV2-tg-1ηut+FOV2(11)GSDhτt,sut=2hut-hterrainτt#pixelhtgFOV/2sin⁡ηutAs an illustrative example of the likelihood model, Figure4(b) shows how the probability of detection, proportional to TTPF, is incremented as the cycle size of the target VT grows. Hence, and due to (8)-(11), the probability of detection grows as the target size Starget is increased, as the UAV flying altitude with respect to the terrain elevation gets smaller (through GSDi(τt,sut)), as the FOV of the camera is reduced (through β(sut) and GSDi(τt,sut)), or as the camera elevation gets closer to the vertical one (through β(sut) and GSDi(τt,sut)). In more detail, Figure 4(c) shows several TTPF curves, corresponding to different FOV and ηut values (indicated in the legend), versus the UAV flying altitude with respect to the terrain elevation (i.e., hut-hterrain(τt)) for a target of Starget=5m (e.g., the boat in Areia Branca).Finally, note that the camera likelihood function defined by (6)-(10) is able to observe partially the borders of the footprints and provide different values for targets of different sizes; cells observed from different UAV heights, camera azimuth, and elevation angles; etc. Hence, our model does not nullify the likelihood that those cells are partially within the footprint but their center falls outside the footprint as in [17, 19], and allows the planner to move the camera from the downward-looking pose assumed in other camera models tested in probabilistic target search problems [13, 14, 17].


### 2.5. UAV and Sensor Dynamic Models

The UAV dynamics corresponds to a fixed-wing UAV modeled with the upper part of the nonlinear parametrizable Simulink model represented in Figure5 or the differential equations of the appendix. The UAV motion variables within sut, highlighted in light green on the right of the figure, are the UAV 3D location (xut,yut,hut), 3D speeds (x˙ut,y˙ut,h˙ut), heading (θut), course angle (θucourse,t), air and ground velocity (vuair,t and vut), and fuel consumption (fuelut). The UAV motion variables within the command aut, highlighted in cyan on the left, are the commanded UAV velocity (vuc,t), heading (θuc,t), and height (huc,t). Additional inputs to the UAV dynamical model are the wind speed (vwt) and direction (θwt), highlighted on the left in pink. Besides, in order to let the reader identify the UAV dynamics within the Simulink model, the blocks associated with its height dynamics are colored in blue, with its velocity in grey, with its lateral displacement in white, with the wind in green, and with the fuel in magenta.Figure 5
Simulink UAV and camera motion modelThe sensor dynamic model is represented at the bottom of Figure5 and corresponds to a gimbaled camera whose scan direction can be changed commanding its elevation (ηuc,t) and azimuth (αuc,t) angles, with the model inputs highlighted in red on the left of the figure. Besides, the camera motion variables within sut are the outputs of the model highlighted in yellow on the right (camera elevation ηut and azimuth αut) and the blocks of the camera dynamics are colored in orange.The whole model includes the usual limitations related to the UAV speed, height, and heading and to the camera elevation and azimuth. Its different parameters (provided in an external input file, eligible by the operator) allow adapting the movement of the UAV and of the camera to the behavior of different real-world aircraft and camera gimbals. Finally, the model integration, performed from each UAV and camera initial statesu0 using the values of a sequence of commands au0:tend, allows obtaining the su0:tend with a given simulation resolution time step Ts.


## 2.1. General Problem Formulation

In the MTS problem presented in this paper, there is a set ofU UAVs overflying a search area Ω (discretized into a grid GΩ of wx·wy square cells) in order to detect the presence of a single target located within it. Besides, due to the uncertainty associated with the position and dynamics of the target and to the measurements of the onboard sensor of each UAV, the MTS problem is formulated within the probabilistic framework introduced in [5].In more detail, the information about the target initial position (represented by random variableτ0) is modeled with the initial target belief b(τ0)≜P(τ0), which is a probability mass function that represents the chances of finding the target at each τ0∈GΩ. Besides, when the target is moving, the information about its dynamics is described with the target Markovian motion model Pτt∣τt-Tτ, which states the probability that the target moves from any cell τt-Tτ∈GΩ to any cell τt∈GΩ during a time lapse Tτ. Moreover, the detection capabilities of the onboard sensor of each UAV are described with the likelihood function Pzut∣τt,sut, which returns the probability of detecting (zut=D) or not detecting (zut=D¯) a target placed at τt from the UAV and sensor pose stored in the state variable sut of the u-th UAV.For those readers familiarized with the Recursive Bayesian Filter (RBF, [12]), the previous probability models b(τ0), Pτt∣τt-Tτ, and Pzut∣τt,sut are used to obtain the target belief b(τt)≜Pτt∣z1:U0:t,s1:U0:t, which represents the chances of finding the target at time step t and at each cell τt∈GΩ, given the trajectory of the UAVs and sensors poses s1:U0:t and the sensor measurements z1:U0:t. The RBF process, stated in (1), calculates the current belief b(τt) from the previous time step belief b(τt-Tτ)≜Pτt∣z1:U0:t-Tτ,s1:U0:t-Tτ by (1) incorporating the current UAVs location and sensor poses s1:Ut and the last sensor measurements z1:Ut with the likelihood functions Pzut∣τt,sut, and by (2) displacing and redistributing the target location probability with the target motion model Pτt∣τt-Tτ. Besides, (1) assumes that the sensors provide measurements only at time steps that are multiple of Tτ and its normalization factor γ is used to ensure that ∑τt∈GΩb(τt)=1.(1)bτt=γ∏u=1:UPzut∣τt,sut∑τt-Tτ∈GΩPτt∣τt-Tτbτt-TτA MTS mission can be formulated as an optimization problem, where some criteria related to the previous probabilistic models and other mission objectives are evaluated in order to determine the best UAV trajectories and sensor posess1:U0:tend during the duration of the mission tend. A useful criterion for MTS is the Expected Time of Detection (ETD, [5–7, 9, 10]), which measures the average time of detection of the target when the UAVs and the sensors follow the trajectories and poses defined in a given s1:U0:tend. It is calculable with (2)-(4), setting b~(τt) to b(τ0) for the initial case (t=0). The first equation is similar to RBF equation (1), but as it lacks the normalization term γ and all measurements z1:U0:t are nondetection, it obtains the “unnormalized belief” b~(τt) assuming that the onboard sensors do not detect the target from the UAVs trajectories and sensor poses in s1:U0:t. The second one obtains Pnd(s1:U0:t)≜PD¯1:U0:t∣s1:U0:t, which is the probability of not detecting the target from s1:U0:t and whose value decays as the sensors make nondetection observations in regions with probability of target presence. The third expression obtains the ETD from s1:U0:tend when Pnd(s1:U0:t) becomes 0 for some t≤tend and underestimates its value when Pnd(s1:U0:tend)>0. Finally, it is worth noting that minimizing ETD(s1:U0:tend) is a better option for MTS than maximizing the probability of detection along the whole trajectory Pd(s1:U0:tend)=1-Pnd(s1:U0:tend), as the MTS objective is to collect as much probability as possible sooner [4–6].(2)b~τt=∏u=1:UPzut=D¯∣τt,sut∑τt-Tτ∈GΩPτt∣τt-Tτb~τt-Tτ(3)Pnds1:U0:t=∑τt∈GΩb~τt(4)ETDs1:U0:tend=∑t=0tendPnds1:U0:t·TτBesides, in order ensure that the UAV trajectories and sensor posess1:U0:tend calculated by our MTS planner are feasible from the maneuverability point of view, we exploit the UAV and sensor deterministic dynamic model of each UAV-sensor pair s˙ut=f(sut,aut), where aut stands for the set of control actions (e.g., commanded UAV heading, height & speed, and sensor elevation & azimuth) and s˙ut for the time derivative of the UAV location and sensor poses. In particular, our MTS planner uses this dynamic model to obtain the solutions (best UAV trajectories and sensor poses s1:U0:tend) from the initial UAVs locations s1:U0 and the sequence of sets of control actions a1:U0:tend proposed, manipulated, and evaluated by our approach in order to optimize MTS missions. Finally, it is worth highlighting that due to the deterministic nature of the UAV and sensor dynamic model, there is an unambiguous relation between the best sequence of control actions and the best UAVs trajectories and sensor poses.The realism of the four models (i.e., of the initial target beliefb(τ0), of the target motion model Pτt∣τt-Tτ, of the sensor likelihood Pzut∣τt,sut, and of the UAV and sensor motion model s˙ut=f(sut,aut)) is crucial to avoid discordances between the real ETD of a given s1:U0:t and the ETD calculated by our approach. Hence, in the rest of this section we present the new models proposed in this paper to bring realism to the definition and solution of MTS missions performed by fixed-wing UAVs equipped with orientable cameras.


## 2.2. Initial Target Belief Definition

To construct the initial target beliefb(τ0), we merge knowledge coming from different sources, using a different probability layer bl(τ0) for each information source and performing with (5) an addition of the L probability layers weighted with their reliability/importance coefficients wl. In other words, considering the first term in (5) is a normalization coefficient that ensures that b(τ0) is a probability function (i.e., ∑τ0∈GΩb(τ0)=1), our initial target belief b(τ0) is calculated as the mixture of the beliefs bl(τ0) associated with the different l-th information sources.(5)bτ0=1∑l=1Lwl∑l=1Lwlblτ0The layers can be associated with geographical information (e.g., terrain altitude, road maps) or intelligence/user-defined clues (e.g., last/habitual location areas of the target). For the examples of this work, we consider the following two layers:(1)
The terrain elevation probability layerb1(τ0), obtained with the following steps. First the digital elevation model (DEM, [29]) of Ω is automatically resampled to the size of the cells in GΩ in order to obtain the average height of each cell in GΩ. Next the user/operator is required to divide the existing elevations within the cells in GΩ into consecutive ranges and to assign a chance of target presence to each range. Finally, the method automatically determines the cells in each elevation range and the probability associated with all τ0∈GΩ, distributing the chances of target presence assigned by the operator. As this way of proceeding generates a geographical probability layer b1(τ0) where areas with similar altitude share the same initial belief, it automatically spreads uniformly the belief over different regions of the search area.(2)
The intelligence probability layerb2(τ0). For this layer the operator has to graphically define a mixture (weighted addition) of Gaussians (centered in eligible locations of the search area Ω and spread according to selectable standard deviations) and of polygonal areas (defined by their external points placed in the desired locations of the search area Ω) with uniform probabilities (assigned by the operator). The weights of each element (each Gaussian and/or each polygonal area) within the intelligence probability layer b2(τ0) are also selected by the operator according to the information gathered about the last known location of the target.As an illustrative example, we show the initial belief defined by an operator when looking for a drifting boat next to the cost in Areia Branca, Brazil. To obtain the probability elevation layerb1(τ0) shown at the bottom of Figure 2(b), where darker/lighter greys are associated with higher/lower probability cells, the operator analyzes the elevation map represented in Figure 2(a) and assigns high chances of target presence to the elevations by the coast (represented in dark green in the elevation map) to model a strong belief that the boat may have arrived a ground when the search mission starts, lower chances to the sea-level elevation (in blue) to model a moderate belief that the target is still navigating, and zero chances to higher inland elevations. In order to define the intelligence probability layer b2(τ0), represented at the top of Figure 2(b), the operator centers a moderate spread Gaussian in the last known position of the boat. Finally, the operator assigns a weight for each layer (w1=2 and w2=1) in order to produce the initial target belief b(τ0) shown in Figure 2(c).Figure 2
Initial target belief definition used for a drifting boat (in Areia Branca, Brazil).
(a)
Elevation map (m) (b)
Probability layers (b1(τ0) at bottom &  b2(τ0) at top) (c)
Initial beliefb(τ0)


## 2.3. Target Motion Model Definition

To define that target motion model, we can distinguish two cases:(1)
Scenarios with static targets. In this case, and due to the immobility of the target,Pτt∣τt-Tτ=1 if and only if τt and τt-Tτ refer to the same cell in GΩ, and Pτt∣τt-Tτ=0 otherwise. Besides, in this case, the equations that contain the term ∑τt-Tτ∈GΩPτt∣τt-Tτb(τt-Tτ) can be simplified as ∑τt-Tτ∈GΩPτt∣τt-Tτb(τt-Tτ)=b(τt-Tτ).(2)
Scenarios with dynamic targets. To construct this type of target motion model, we combine different types of information (e.g., elevation data and sea currents). In this work, the operator has to provide the following:
(i)
The elevation rangeER where the target is allowed to move within Ω. This allows the operator to indirectly assign the probabilities corresponding to the static target behavior to those cells that do not belong to the ER. In other words, ∀τt-Tτ∉moving(GΩ,ER), the target motion definition process automatically makes Pτt∣τt-Tτ=1 if τt=τt-Tτ and Pτt∣τt-Tτ=0 otherwise.(ii)
K vectors v1:K, each with 9 values that represent, for a few selected cells c1:K⊂GΩ, the chances of moving to their 8 neighbor cells N(ck) and of staying at the same cell ck. Besides, in order to force the target to stay within the search area Ω, some of the 9 values of the vectors vk that correspond to cells in the borders of GΩ are automatically set to zero. With this information, the target motion definition process makes Pτt∣τt-Tτ=0 for all τt∉N(τt-Tτ)∪τt-Tτ and computes Pτt∣τt-Tτ applying the potential field method (according to the cardinal directions that make the target move from one cell to its neighbors) for all τt∈N(τt-Tτ)∪τt-Tτ and τt-Tτ∈moving(GΩ,ER). In more detail, if adj(τt-Tτ,d) with d∈1:8 stands for the adjacent cell to τt-Tτ in the d-th cardinal direction, adj(τt-Tτ,9)=τt-Tτ represents stay in the same cell, dist(τt,ck) is the distance between cell τt and cell ck, and vk[d] is the d-th element of vk, then Pτt=adjτt-Tτ,d∣τt-Tτ∝∑k=1:Kvk[d]/dist(τt,ck) and ∑d∈1:9Pτt=adjτt-Tτ,d∣τt-Tτ=1. This way of proceeding allows the operator to define the distribution of the target probability around the neighborhood cells in only a few cells c1:K⊂GΩ and extend this definition to the remaining moving cells taking into account the distance between the selected cells and the others.(iii)
The time periodTτ that has to be used to apply Pτt∣τt-Tτ to the belief b(τt-Tτ). As the speed of the target is related to the quotient of the size of the cell and Tτ, this value is used by the planner to relate the simulation of the motion of the target to the simulation of the motion of the UAVs and sensors.As an illustrative example, we also obtain the target motion model of the boat drifting example next to the cost in Areia Branca, Brazil. In this case, the operator definesvk for three cells ck∈GΩ in the sea, making vk[9]=0.1, vk[d]=2 for one of the cardinal directions d towards the coast, and vk[r]=0 for the remaining directions r∈1:8-d. Figure 3 represents with the starting point of the red arrows the location of the selected c1:3 and with the arrows themselves the direction d where vk[d]=2 for each ck. Besides, the operator also makes the moving ER≤0 to only let the target (i.e., the drifting boat) move in the sea. To summarize the result of the target motion definition process, the green arrows in Figure 3 represent the direction obtained by weighting, for each cell τt-Tτ∈moving(GΩ,ER), the arrows in the 8-cardinal directions d with Pτt=adjτt-Tτ,d∣τt-Tτ. Besides, the lack of inland green arrows indicates the static behavior of the target in this area. In short, the green arrows show how the defined target motion makes the belief over the sea move towards the coast and remain unchanged inland. Finally, it is worth highlighting that, thanks to our target motion model definition approach, this complex model is automatically built just considering the elevation range (ER≤0) and the target movement probabilities in three cells of GΩ.Figure 3
Target motion definition used for a drifting boat (in Areia Branca, Brazil).


## 2.4. Camera Likelihood with Terrain Occlusion

The likelihood of detecting the target with the onboard camera of each UAV in a given cellτt∈Ω can be calculated scaling the Target Task Performance (TTP) metric [30] with the percentage of the selected cell within the camera footprint. This scaling, used to reduce the detection probability in the cells of the footprint border, can be modeled with (6), where τt∩footprintsut represents the area (size of the surface) of the common region of the cell τt and of the footprint of the camera (oriented and placed according to sut), τt is the total area of cell τt, and TTPF(τt,sut) is the target task performance function.(6)PDut∣τt,sut=τt∩footprintsutτtTTPFτt,sutTo determine the footprint of the camera from the UAV location and camera pose insut, we follow the process schematized in Figure 4(a) and consisting of the following two steps. First, we calculate geometrically the camera footprint at sea level taking into account the following pieces of information within sut: UAV location (xut,yut,hut) and camera azimuth, elevation, and FOV (αut,ηut,FOV). To do this, we consider that the sensor location is the same as the UAV location (since its deviation from the location of the UAV mass center is negligible for the mission) and that its orientations are measured with respect to the vehicle coordinate frame (since the camera gimbal compensates the UAV attitude). Second, we approximate the corners of the real camera footprint by the intersections (obtained using the 3D Bresenham algorithm [31]) of the terrain elevation with the 4 lines that join the camera with the 4 corners of the sea-level camera footprint.Figure 4
Camera likelihood model.
(a)
Footprint determination (b)
TTPF curve vs.VT (c)
TTPF curves vs. relative height (forStarget=5m)The value of the target task performance function TTPF for each cellτt within the sensor footprint from the UAV location and sensor pose in sut is calculated with (7), where VT(τt,sut) is the cycle size of the target and V50=2.03 (according to [32, 33]) is the critical cycle size for detecting a small target with a likelihood of 50% (since when VT=V50, TTPF=1/(1+1)=0.5). Moreover, the target cycle size is calculated with (8), where Starget is the real size of the target, β(τut) is the angle between the along-scan (vertical) and cross-scan (horizontal) directions of the footprint, and GSDi (with ∈vertical,horizontal) is the ground sample distance (corresponding size of a pixel of the camera over terrain, obtainable as the ratio of the length of the footprint in each direction to its corresponding number of pixels). Taking into account the trigonometric relations between the UAV location and camera footprint, the angle β(sut) can be obtained solving (9) and the ground sample distance in each scan direction can be calculated with (10) and (11), where hterrain(τt) is the height of the terrain at the target location τt and #pixeli is the number of pixels of the camera in the vertical or horizontal directions (i∈v,h).(7)TTPFτt,sut=VTτt,sut/V501.51+0.24VTτt,sut/V501+VTτt,sut/V501.51+0.24VTτt,sut/V50(8)VTτt,sut=StargetGSDvτt,sut·GSDhτt,sut·sin⁡βsut(9)tgβsut=tg-1ηut-FOV/2-tg-1ηut+FOV/2tgFOV/2sin-1ηut-FOV/2-sin-1ηut+FOV/2(10)GSDvτt,sut=hut-hterrainτt#pixelvtg-1ηut-FOV2-tg-1ηut+FOV2(11)GSDhτt,sut=2hut-hterrainτt#pixelhtgFOV/2sin⁡ηutAs an illustrative example of the likelihood model, Figure4(b) shows how the probability of detection, proportional to TTPF, is incremented as the cycle size of the target VT grows. Hence, and due to (8)-(11), the probability of detection grows as the target size Starget is increased, as the UAV flying altitude with respect to the terrain elevation gets smaller (through GSDi(τt,sut)), as the FOV of the camera is reduced (through β(sut) and GSDi(τt,sut)), or as the camera elevation gets closer to the vertical one (through β(sut) and GSDi(τt,sut)). In more detail, Figure 4(c) shows several TTPF curves, corresponding to different FOV and ηut values (indicated in the legend), versus the UAV flying altitude with respect to the terrain elevation (i.e., hut-hterrain(τt)) for a target of Starget=5m (e.g., the boat in Areia Branca).Finally, note that the camera likelihood function defined by (6)-(10) is able to observe partially the borders of the footprints and provide different values for targets of different sizes; cells observed from different UAV heights, camera azimuth, and elevation angles; etc. Hence, our model does not nullify the likelihood that those cells are partially within the footprint but their center falls outside the footprint as in [17, 19], and allows the planner to move the camera from the downward-looking pose assumed in other camera models tested in probabilistic target search problems [13, 14, 17].


## 2.5. UAV and Sensor Dynamic Models

The UAV dynamics corresponds to a fixed-wing UAV modeled with the upper part of the nonlinear parametrizable Simulink model represented in Figure5 or the differential equations of the appendix. The UAV motion variables within sut, highlighted in light green on the right of the figure, are the UAV 3D location (xut,yut,hut), 3D speeds (x˙ut,y˙ut,h˙ut), heading (θut), course angle (θucourse,t), air and ground velocity (vuair,t and vut), and fuel consumption (fuelut). The UAV motion variables within the command aut, highlighted in cyan on the left, are the commanded UAV velocity (vuc,t), heading (θuc,t), and height (huc,t). Additional inputs to the UAV dynamical model are the wind speed (vwt) and direction (θwt), highlighted on the left in pink. Besides, in order to let the reader identify the UAV dynamics within the Simulink model, the blocks associated with its height dynamics are colored in blue, with its velocity in grey, with its lateral displacement in white, with the wind in green, and with the fuel in magenta.Figure 5
Simulink UAV and camera motion modelThe sensor dynamic model is represented at the bottom of Figure5 and corresponds to a gimbaled camera whose scan direction can be changed commanding its elevation (ηuc,t) and azimuth (αuc,t) angles, with the model inputs highlighted in red on the left of the figure. Besides, the camera motion variables within sut are the outputs of the model highlighted in yellow on the right (camera elevation ηut and azimuth αut) and the blocks of the camera dynamics are colored in orange.The whole model includes the usual limitations related to the UAV speed, height, and heading and to the camera elevation and azimuth. Its different parameters (provided in an external input file, eligible by the operator) allow adapting the movement of the UAV and of the camera to the behavior of different real-world aircraft and camera gimbals. Finally, the model integration, performed from each UAV and camera initial statesu0 using the values of a sequence of commands au0:tend, allows obtaining the su0:tend with a given simulation resolution time step Ts.


## 3. MTS Planner

The MTS planner presented in this work is an extension of the Genetic Algorithm (GA) planner introduced in [10]. This section details the main elements of the new planner, introduces its algorithmic description following a top-down strategy, and discusses its computational cost.


### 3.1. Multistepped GA-Based Planner

The MTS planner obtains the UAVs trajectoriess1:U0:tend for a given mission time tend following a receding horizon controller approach [10, 13, 16, 17, 19, 20]. More concretely, the planner divides s1:U0:tend into R different sections of duration Q and optimizes each section (s1:UQ(r-1):Qr with r=1:R) sequentially, using a GA whose inputs are the final state of the last section s1:UQ(r-1) and its “unnormalized belief” b~(τQ(r-1)), and whose outputs are the new s1:UQ(r-1):Qr and b~(τQr). Besides, it incorporates the possibility of letting each UAV engage and leave the MTS mission at different time instants tustart and tuend (without loss of generality, we assume that for at least one UAV, tustart=0) and of using different time steps for the measures of each UAV Tuz and for the target update Tτ. Hence, to obtain sequentially the sections of the UAV trajectories (and of the cameras orientations), the new planner obtains the overall mission ending time, calculates the number of sections R required for the mission, and considers the UAVs and camera immobile and disabled for those time steps where they are not engaged in the MTS mission (i.e., su0:tustart=sutustart and sutuend:tend=sutuend). Besides, it also obtains the minimum time step Tb required to update the “unnormalized belief” b~(τt) during the evaluation of the ETD taking into account the time steps of target motion model Tτ and of the camera of each UAV Tuz.The algorithmic description of the new multistepped planner is sketched in the pseudocode of Algorithm1. It starts computing the mission duration tend, the number of sections R into which s1:U0:tend will be divided, and the time lapse Tb required to update the “unnormalized belief” b~(τt). Next, it initializes the “unnormalized belief” b~(τ0) and the fitness criteria vector FCV0, which stores the different values of the criteria used by the GA to decide if a solution (proposal of UAV trajectories and sensor poses in s1:U0:t) is better than another. Finally, it sequentially optimizes each of the sections of the trajectory (i.e., each st:t+Q with t=Q(r-1) and r={1:R}) using the optimizer in optimGA(·), which will be described in detail in the following section.Algorithm 1:Multistepped GA-based planner.
Require:  sc⊳Structure with scenario models and times: b(τ0), Pτt∣τt-1, Pzt∣τt,sut, f(sut,aut), t1:Ustart, t1:Uend Require:  op⊳Structure with optimizer parameters Require:  s1:Utustart⊳Initial UAVs locations and sensor poses ∀u={1:U} Require:   Q⊳ Duration of each subsequence 1:tend=max⁡sc.t1:Uend⊳ Determine the duration of the mission 2:R=tend/Q⊳ Determine the number of sections required for the multi-stepped planner 3:tend=Q⋅R⊳ Round the duration of the mission 4:sc.Tb=obtainT(sc.T1:Uz,sc.Tτ)⊳ Determine minimum time step to update b~(τt). 5:t←0⊳ Set starting time 6:b~(τt)=sc.b(τ0)⊳ Initialize b~(τt) 7:FCVt=[0,0,0,0,0,0]⊳ Initialize fitness criteria vector 8:for r=1:R do⊳ For each subsequence 9:[s1:Ut:t+Q,FCVt+Q,b~(τt+Q)]←optimGA(sc,op,s1:Ut,b~(τt),t,Q,FCVt)⊳ Obtain best subsequence 10:t←t+Q⊳ Update time for the following section 11:end for 12:return  s1:U0:tend⊳ Complete UAV trajectory and sensor poses


### 3.2. GA Optimizer

This section details the most relevant aspects of the GA used to determine the best UAV trajectory and sensor pose for each section (st:t+Q with t=Q(r-1)) of the whole solution s0:tend. We start presenting the GA genotype, continue with the main GA steps and operators, and finish with the evaluation of each possible solution s0:tend.


#### 3.2.1. GA Decision Variables

The decision variables directly manipulated by the genetic operators (i.e., the GA genotype) of the MTS planner are the subsequence of control actionsa1:UQ(r-1):Qr used to obtain, by means of the UAV motion model, the r-th subsection of the UAV trajectories s1:UQ(r-1):Qr (i.e., the GA phenotype or solution). Besides, the actions in aut are the commanded UAV heading θuc,t, speed vuc,t, and height huc,t and the camera elevation ηuc,t and azimuth αuc,t. Additionally, the number of decision variables is reduced by applying each action constantly during a time lapse Ta and by adapting the actions required for each UAV within each r-th subsection to the time steps where the UAV is engaged in the MTS mission. In other words, the GA algorithm only obtains the values of Q/Ta sets of actions (UAV heading, speed, and height; camera elevation and azimuth) for each UAV engaged in the mission during the r-th trajectory subsection.


#### 3.2.2. GA General Description

The GA that sequentially optimizes each of the trajectory sectionss1:UQ(r-1):Qr based on the set of command actions (aut=[θuc,t,vuc,t,huc,t,ηuc,t,αuc,t]) and fitness criteria (which will be explained in the following section) maintains the operators, parameters, and structure of a standard GA and includes a new step to precalculate the number of actions required by each UAV.To present it in an algorithmic form, we extend the notation of the paper with left superindexes (1:J or best) and left subindexes (par and child). The first superindex is used to indicate that there is a population or group of J action sequences 1:Ja1:Ut:t+Q and fitness criteria vectors 1:JFCVt+Q. The second one is used to identify the best action sequence besta1:Ut:t+Q, UAV trajectories and sensor poses bests1:Ut:t+Q, and fitness criteria vector bestFCVt+Q within the population of 1:J solutions. And the left subindexes are incorporated to the previous notation when the algorithm needs to distinguish between the generic population elements (1:Ja1:Ut:t+Q, 1:JFCVt+Q) and those associated with the parents selected by the GA (par1:Ja1:Ut:t+Q) and with their children (child1:Ja1:Ut:t+Q, child1:JFCVt+Q).The GA in Algorithm2 performs the following steps. First, step 3 computes the number of actions of Ta duration required by each UAV in the current subsection of the trajectory, taking into account the current time step t and each UAV starting and ending engaging time (tustart and tuend). Step 4 initializes a random population of the J action sequences 1:Ja1:ut:t+Q required for the UAVs engaged in the mission, using a parametrizable uniform distribution. Steps 5-7 simulate the UAV trajectories and sensor poses from the sequence of actions and evaluate them according to the criteria presented in the following section. In step 9, the GA enters an optimization loop which ends when a fixed computation time has passed. Next, step 10 selects J parents of action sequences from the current population using binary tournament selection; step 11 creates the children by combining pairs of parents using a single-point crossover; step 12 slightly modifies the children with an incremental Gaussian noise mutation; steps 13-15 simulate and evaluate them; and step 16 selects the survival population from the previous and the mutated children using NSGA-II recombination [34]. Once the computation time has finished, step 19 identifies the best solution in the population; step 20 obtains the FCV, trajectories and sensor poses, and “unnormalized belief” associated with it; and step 21 returns them.Algorithm 2:optimGA: NSGA-II for MTS.
1:procedure  optimGA(sc, op, b~(τt),s1:Ut,t,Q,FCVt) 2:J←op.J⊳ Get population size 3:N1:U←obtainActionNumber(t,Q,Ta,sc.t1:Ustart,sc.t1:Ustart)⊳ Obtain number of actions for each UAV 4:1:Ja1:Ut:t+Q←initialize(J,N1:U,op.init)⊳ Generate J sequences with the required actions 5:for j=1:J do 6:jFCVt+Q←sim&Eval(sc,b~(τt),s1:Ut,t,Q,FCVt,ja1:Ut:t+Q)⊳ Simulate and evaluate 7:end for 8:M←1⊳ Initialize counter of the algorithm iterations (generations) 9:while  getCompTime()≤opt.CompTime  do⊳ Main optimization loop 10:par1:Ja1:Ut:t+Q←selectParents(1:Ja1:Ut:t+Q,1:JFCVt+Q)⊳ Select pair of parents (tournament selection) 11:child1:Ja1:Ut:t+Q←crossover(par1:Ja1:Ut:t+Q,op.xover)⊳ Perform crossover of parents (single point crossover) 12:child1:Ja1:Ut:t+Q←mutate(child1:Ja1:Ut:t+Q,op.mut)⊳ Mutate children (incremental gaussian mutation) 13:for j=1:J do 14:childjFCVt+Q←sim&Eval(sc,b~(τt),s1:Ut,t,Q,FCVt,childja1:Ut:t+Q)⊳ Simulate and evaluate 15:end for 16:[1:Ja1:Ut:t+Q,1:JFCVt+Q]←recombine(1:Ja1:Ut:t+Q,1:JFCVt+Q,child1:Ja1:Ut:t+Q,child1:JFCVt+Q) 17:M←M+1⊳ Increment the GA generation number 18:end while 19:besta1:Ut:t+Q←best(1:Ja1:Ut:t+Q,1:JFCVt+Q)⊳ Select best solution 20:[bestFCVt+Q,bests1:Ut:t+Q,bestb~(τt+Q)]←sim&Eval(sc,b~(τt),s1:Ut,t,Q,FCVt,besta1:Ut:t+Q)⊳ Sim. and evaluate it 21:returnbests1:Ut:t+Q,bestFCVt+Q,bestb~(τt+Q)⊳ Return best solution 22:end procedureIt is worth noting that there are several steps whose behavior can be configured using the information in theopt configuration variable. The initialization parameters opt.init in step 4 allow, for each UAV and type of action, (1) selecting different bounds for the uniform distribution and (2) enabling/disabling its optimization. These options let our MTS planner generate trajectories with some fixed behavior (e.g., fixed height or fixed camera orientations). The stop condition parameter opt.CompTime in step 9 allows fixing the computation time of the GA algorithm. The crossover parameter opt.xover in step 11 allows selecting the percentage of parents that undergo the single-point crossover or are directly copied as children. Finally, the mutation parameter opt.mut in step 12 allows selecting the Gaussian noise that is added to all variables in a1:Ut:t+Q and the Gaussian noise that is added to a few variables uniformly selected in 1:Ja1:Ut:t+Q with a probability 1/a1:Ut:t+Q, where a1:Ut:t+Q is the number of decision variables in a1:Ut:t+Q.Besides, there are three steps of the GA that use the fitness criteria variable1:JFCVt+Q to select solutions of the population to become pairs of parents (step 10), survive in the next generation (step 16), and become the best solution (step 19). Its behavior will be further detailed in the following section, after explaining the optimization criteria used in our MTS planner.


#### 3.2.3. Simulation and Evaluation Process

In order to evaluate one of the solutions, subsequence of action commandsja1:Ut:t+Q, manipulated and proposed by the GA, first, the motion model in Figure 5 or in the appendix has to be simulated for each UAV, and, second, the fitness criteria have to be evaluated. The planner has also to take into account the UAVs mission engaging time [tustart,tuend]∀u=1:U and the possible differences among the time lapses of the problem (Tτ, Ta, Ts, Tuz). In the remainder of this section we describe the simulation process and the fitness criteria evaluation, and the steps of the sim&Eval function, presented in Algorithm 3 and used by the GA in Algorithm 2.Algorithm 3:sim&Eval: simulation and evaluation for MTS.
1:procedure  sim&Eval(sc,b~(τt),s1:Ut,t,Q,FCVt,a1:Ut:t+Q) 2:tinit←t 3:tlast←t+Q 4:⊳ Simulation of the proposed sequence of actions to obtain UAV trajectory and fuel consumption 5:for u=1:U do 6:    [sut:t+Q,fuelu]←sim(sc.f(⋅),sut,aut:t+Q,t,t+Q, ⊳ Simulate each UAV and sensor movement 7:sc.tustart,sc.tuend,sc.Ta,sc.Ts) 8:end for 9:⊳ Evaluation of the trajectory obtained from the proposed sequence of actions 10:FCVt[COL]←FCVt[COL]+∑u=1:U‍∑l=u+1:U‍collision(sc,sut:t+Q,slt:t+Q)⊳ Update collision criterion 11:FCVt[NFZ]←FCVt[NFZ]+∑u=1:U‍withinNFZ(sut:t+Q)⊳ Update fuel consumption 12:FCVt[SMOOTH]←FCVt[SMOOTH]+∑u=1:U‍smooth(aut:t+Q)⊳ Update smooth criterion 13:FCVt[FUEL]←FCVt[FUEL]+∑u=1:U‍fuelu⊳ Update fuel consumption 14:⊳ Simulation of the target and evaluation of the probabilistic criteria 15:t←tinit+sc.Tb 16:while  t≤tlast  do 17:if  motionTarget(t,sc.Tτ)  then 18:b~(τt-sc.Tτ)=b~(τt)⊳ Rename variable for the target transition 19:b~(τt)←∑τt-sc.Tτ∈GΩ‍P(τt∣τt-sc.Tτ)b~(τt-sc.Tτ)⊳ Target motion part of Equation (2) 20:end if 21:b~(τt)←∏u∈measuring(u,t,sc)‍P(zut=D¯∣τt,sut)b~(τt)⊳ Measurement update part of Equation (2) 22:Pnd(s1:U0:t)←∑τt∈GΩ‍b~(τt)⊳ Equation (3) 23:FCVt[ETD]←FCVt[ETD]+Pnd(s1:U0:t)∗sc.Tb⊳ Update ET, with equation (4) in sequential form 24:t←t+sc.Tb 25:end while 26:FCVt[MYOP]=∑τt∈GΩ‍∏u=∈engaged(u,t,sc)‍(1-λdistance(τt,sut))b~(τt)⊳ Calculate myopia reduction criterion 27:return  FCVt, s1:Utstart:tend, b~(τt) 28:end procedureThe simulation of each UAV is performed in step 6 of Algorithm3 with the sim(·) function, taking into account the Simulink or appendix model f(sut,aut), the initial UAV location sut, and the timely-equispaced set of actions in aut:t+Q. Moreover, during the simulation each set of actions aut proposed by the GA is considered constant during the time lapse Ta and the resulting UAV trajectory and sensor poses sut:t+Q are created with sampling simulation time of Ts. We use two different time lapses in the simulator (for Ts and Ta) to reduce the number of decision variables aut:t+Q manipulated by the planner, while permitting a high time resolution for evaluating certain aspects of the UAVs trajectories s1:Ut:t+Q. Finally, it is worth noting that sim(·) also obtains, from the UAV motion model, the fuel consumption of the UAV during the corresponding trajectory subsection and checks if the UAV is engaged in the mission in order to only simulate its displacement in the simulation time steps that simultaneously belong to [t,t+Q] and [tustart,tuend].The fitness criteria of the MTS planner consists in fulfilling two constraints (UAV collision and nonflying zones avoidance) and in optimizing four indexes (command smoothness, fuel consumption, expected detection time, and myopia reduction criterion), whose values are stored in the different elements of the Fitness Criterium VectorFCVt. Besides, in order to distinguish among the problem time periods Ta and Ts within the expressions of the evaluation criteria, we make [t:Δ:t+Q] stand for the time vector with elements equi-spaced by the selected Δ and k∈[t:Δ:t+Q]∧[a,b] stand for the time steps k within the set [t:Δ:t+Q] that fall within the limits [a,b]. Moreover, the tag used in Algorithm 3 within FCVt[tag] expressions permits the identification of the element of the vector FCVt associated with each criterion. Finally, the details of each criterion are presented below:(i)
UAV collision avoidance. This criterion counts the number of trajectory time instants at which the distance between the locations of all pairs of UAVs is smaller than the security collision distance dCOL. The criterion, evaluated in line 10 of Algorithm and detailed in (12), consists in checking the inter-UAV distances at time steps where both vehicles are engaged in the mission.(12)collisionsc,sut:t+Q,slt:t+Q=∑k∈t:Ts:t+Q∧sc.tustart,sc.tuend∧sc.tlstart,sc.tlendxuk-xlk2+yuk-ylk2+huk-hlk2≤sc.dCOL2(ii)
Nonflying zones (NFZ) avoidance. This criterion counts the number of trajectory time instants at which the locations of any UAV overfly a selected set of forbidden overflying cells of GΩ. Again, the criterion, evaluated in line 11 of Algorithm 3 and detailed in (13), considers locations at the time steps where the UAVs are in the mission.(13)withinNFZsc,sut:t+Q=∑k∈t:Ts:t+Q∧sc.tustart,sc.tuendinNFZxuk,yuk(iii)
Actions smoothness. This criterion, evaluated in line 12 of Algorithm 3, measures the smoothness of the commands for each UAV within aut:t+Q. It accounts, with (14) and (15), for the consecutive changes of the five types of commands only during the time interval where the UAVs are engaged in the mission. The magnitude of the accounted change, calculated with (16), makes the smooth criterion only account for those changes in action that imply a modification of sign in the signal with a value proportional to the square of the difference. Therefore, the smoothness criterion penalizes big changes of values and sign in the sequence of each type of command.(14)smoothsc,aut:t+Q=∑k∈t:Ta:t+Q-Ta∧sc.tustart,sc.tuendchangeAauk,auk+Ta(15)changeAaut,auk=chθuc,t,θuc,k+chhuc,t,huc,k+chvuc,t,vuc,k+chηuc,t,ηuc,k+chαuc,t,αuc,k(16)cha,b=a-b2sign⁡a-sign⁡b==2(iv)
Fuel consumption. Its value is directly returned by the UAV motion model and affected by the UAV speed.(v)
Expected time of detection. The main objective of the planner, described already in the general problem formulation section, is implemented in an iterative form within lines 16 and 25 of Algorithm 3. The process implements (2)-(4) after the following modifications. On the one hand, the expression in line 19 includes the target motion model transition operation for the time steps where the motion model has to be applied (i.e., when motionTarget(t,Tτ) is true, because the target is dynamic and the last target motion update was performed Tτ seconds ago). On the other hand, the “unnormalized belief” update expression in line 21 only includes the likelihood models of the UAVs at the time steps t where the UAV is engaged in the mission [tustart,tuend] and its sensor has to take a measure (according to its Tuz).(vi)
Myopia reduction criterion. The planner can create myopic solutions due to its multistepped behavior, which only considers the probability regions achievable from the starting position of each section of the trajectory. In other words, by optimizing only the previous criteria, the planner becomes greedy, ignoring during the optimization of each trajectory section the benefits of getting closer to a high probability area which is currently unreachable. To reduce this myopic behavior, the planner evaluates, in line 26 of Algorithm 3, the suitability of the final UAVs positions and sensor poses s1:UrQ of a trajectory section for the optimization of future trajectory sections by weighting the remaining unnormalized belief b~(τrQ) with a function that takes into account the distance of each cell τt in GΩ to the center of the camera footprint associated with the final location and camera pose of each UAV sut. Besides, this evaluation is only performed for the UAVs that are still engaged in the mission. In other words, the planner evaluates the suitability of the final camera footprints of the searching UAVs for optimizing the remaining sections of the trajectories.As we have previously mentioned, theFCVt values obtained after the simulation of each trajectory section are later used in 3 different steps of the GA. The way of proceeding in each step is the following:(i)
The binary tournament selection operator chooses randomly, with an uniform probability, two solutionsj of the population; among them it selects as a parent the one with lower constraint violation (i.e., with lower jFCVt[COL]+FjCVt[NFZ]), and when tied, the one that Pareto-dominates the other four optimization objectives (i.e., jFCVt[SMOOTH], jFCVt[FUEL], jFCVt[ETD], jFCVt[MYOP]).(ii)
The survival/recombination operator organizes the solutions in a similar way: first according to its constraint violation level and second according to the Pareto-dominance of the remaining four optimization objectives. Hence, as a whole, the GA tries to identify the solution that violates the constraints less times and, within them, those solutions that Pareto-dominate the other objective criteria.(iii)
The selection of a unique best solution is required at the final step of the GA to let the multistepped planner use its final position as the starting point of next trajectory section. To do this, the GA selects, among the solutions in the best Pareto-front of the last population obtained by the recombination operator, the one with smallest rounded value of myopia reduction criterion (first), smallest ETD (second, if tied in the first), smallest command smoothness (third, if tied in the previous), and small fuel consumption (fourth). The priority levels used for the selections process reflect the importance order of the objectives. Moreover, the myopia criterion is preferred to ETD, becauseFCVt[MYOP] accounts both for the probability gathered by the trajectory (through b~(τt)) and for the distance to all the cells with remaining “unnormalized belief”, while FCVt[ETD] only considers the reachable probability within each trajectory section s1:Ut:t+Q.


### 3.3. Computational Cost of the MTS Planner

Finally, we calculate the computational complexity of our planner. To do this, we will consider that all UAVs start and end the mission at the same time steps (t1:Ustart=0, t1:Uend=tend) and that the values of quotients Q/Ts, Q/Ta, Q/Tτ, Q/Tuz, tend/Q are natural numbers.We start studying Algorithm3 and estimate the computation complexity of the UAV and sensor motion simulation (ccsim) and of the evaluation of each fitness criterion (cctag). The computational cost of the UAV and sensor motion simulation is proportional to the number of UAVs U and to the steps Q/Ts in which the simulator has to return UAV and sensor poses, ccsim=O(U·Q/Ts). The computation cost of the collision avoidance criterion is proportional to the number of possible UAV combinations U(U-1) and to the steps Q/Ts in which the UAV trajectory and sensor poses have been discretized. Hence, ccCOL=O(U2·Q/Ts). For the NFZ criterion, we can follow a similar process, but accounting only for the number of UAVs. Hence ccNFZ=O(U·Q/Ts). In the case of the control smoothness, we have to consider the number of UAVs U, the number of steps Q/Ta in which the action sequence is discretized, and the number of types of actions that are being optimized A. Hence ccSMOOTH=O(U·A·Q/Ta). As the fuel consumption is calculated during the simulation of the UAV and sensor motion, its computational complexity is already included in ccsim. To estimate the complexity associated with the ETD computation, we can observe that the cost associated with the target motion is proportional to Q/Tτ and (wxwy)2, and cost associated with the sensor update is proportional to (wxwy)∑u=1UQ/Tuz. The first cost is associated to the operations in the target transition ∑τt-Tτ∈GΩPτt∣τt-Tτb~(τt-Tτ), which can be calculated, a number of Q/Tτ times, as the multiplication of two matrixes of (wxwy)x(wxwy) and (wxwy)x1 size. The second cost is associated to the update of each sensor Pzut∣τt,sutb~(τt), which is calculated a number of Q/Tuz times and is proportional to wxwy. Hence, ccETD=O((wxwy)2Q/Tτ+(wxwy)∑u=1UQ/Tuz). Finally, as the computational cost of the myopia reduction criterion is proportional to U and to wxwy, its complexity ccMYOP=O(U·wxwy). The addition of the previous complexities, after simplifying the final expression, makes the expression of the computational complexity of the simulation and evaluation function become (17).(17)ccsim&Eval=OU2QTs+U·A·QTa+wxwy2QTτ+wxwyU+∑u=1UQTuzThe expression shows that the complexity of the simulation and evaluation function ccsim&Eval grows with the increment of the duration Q of each evaluation, the number of UAVs U, the number of decision variables manipulated by the planner A, and the number of cells of the grid of the search space wxwy. Besides, ccsim&Eval also grows as the time steps of the simulation Ts, of the actions Ta, of the target Tτ, and of the camera measurements Tuz are decremented. Hence, in order to have a reasonable computational cost, it is necessary to achieve a compromise among all those parameters. Besides, the quadratic dependency of ccsim&Eval with wxwy (i.e., with the size of the grid of the search space) is worth noting. This quadratic cost is especially relevant, because a moderate size search space (e.g., of 30 x 30 cells in the Areia Branca example) has already significantly bigger values of wx·wy (e.g., 900) than those used in other parameters (e.g., U=2, A=3, Q/Ts=50, Q/Ta=Q/Tuz=5, Q/Tτ=1).To account for the computational complexity of the planner, we still have to consider the computational cost of the GA in Algorithm2 and of the multistepped approach in Algorithm 1. For the GA, we compute the cost of the main loop of Algorithm 2. Considering a population size of J solutions, U·A·Q/Ta decision variables, and M iterations of the loop, the cost of the ccGAselparent=O(M·J), ccGAxover=ccGAmut=O(M·J·U·A·Q/Ta), ccGAsim&Eval=M·J·ccsim&Eval, and ccGArecombine=O(MJ2). This happens because computational cost of the tournament selection of each pair of parents is proportional to a constant, the cost of the single-point crossover and of the incremental mutation of a solution is proportional to the number of decision variables within it, and the cost of the recombination step of the population, which includes the sorting of the solutions, is proportional to J2. Hence, the overall GA complexity can be calculated as ccGA=M·Jccsim&Eval+O(U·A·Q/Ta+J). Finally, as the multistepped GA-based planner in Algorithm 1 calls the GA R=tend/Q times, its overall computation complexity becomes ccMS=tend/Q·M·Jccsim&Eval+O(U·A·Q/Ta+J), calculable with (18), obtained by substituting within it ccsim&Eval by (17).(18)ccMS=OtendQ·M·JJ+U2QTs+U·A·QTa+wxwy2QTτ+wxwyU+∑u=1UQTuzThe overall computational complexity shows the same dependencies observed for the computational complexity of the simulation and evaluation and a (multiplying) growing behavior proportional to the population sizeJ, the number of generations in GA M, and the number of subsections in which the trajectory is divided tend/Q. Besides, as Q appears in the denominator of tend/Q that multiplies the numerator of different quotients with the time steps, there is a part of the computational cost that can be considered independent of Q and of the number of sections in R=tend/Q in which the trajectories are divided, and dependent on the quotient of the duration of the mission tend and the different time steps (Ts, Ta, Tτ, and Tuz).Finally, it is worth noting that the complexity analysis, which is useful to consider the effects of selecting certain variables of the system in the computation cost, does not provide information about (1) the time required for evaluating the solutions or (2) the convergence speed of the multistepped planner and of each GA run, whose behavior is highly dependent on the probability models (initial belief, target motion model, and sensor detection function) of each scenario. To analyze the convergence speed, in the Results section we will analyze the temporal performance of our new planner (implemented in Matlab and using its Parallel Toolbox in the evaluation loops, lines 5-7 and 13-15, of the GA algorithm) over different scenarios to be able to determine, empirically, the time needed to optimize them (in a 2.5 GHz Intel Core i7 with 6 GB RAM PC with Windows 7). To provide some insight related to the time required for evaluating the solutions within the planner, we have measured the values of the different steps of the GA algorithm in the scenarios of the following section and observed that the mean computation time of each population step of the GA (loop in line 9 of Algorithm2) is 0.46 seconds for the first scenario, 0.49 seconds for the second one, and 1.05 seconds for the third one; that the evaluation of the solutions (sim&Eval function) takes 99% of the time; and that the UAV motion simulation takes 43% of the time of the sim&Eval function, while the rest of it takes 53%.


## 3.1. Multistepped GA-Based Planner

The MTS planner obtains the UAVs trajectoriess1:U0:tend for a given mission time tend following a receding horizon controller approach [10, 13, 16, 17, 19, 20]. More concretely, the planner divides s1:U0:tend into R different sections of duration Q and optimizes each section (s1:UQ(r-1):Qr with r=1:R) sequentially, using a GA whose inputs are the final state of the last section s1:UQ(r-1) and its “unnormalized belief” b~(τQ(r-1)), and whose outputs are the new s1:UQ(r-1):Qr and b~(τQr). Besides, it incorporates the possibility of letting each UAV engage and leave the MTS mission at different time instants tustart and tuend (without loss of generality, we assume that for at least one UAV, tustart=0) and of using different time steps for the measures of each UAV Tuz and for the target update Tτ. Hence, to obtain sequentially the sections of the UAV trajectories (and of the cameras orientations), the new planner obtains the overall mission ending time, calculates the number of sections R required for the mission, and considers the UAVs and camera immobile and disabled for those time steps where they are not engaged in the MTS mission (i.e., su0:tustart=sutustart and sutuend:tend=sutuend). Besides, it also obtains the minimum time step Tb required to update the “unnormalized belief” b~(τt) during the evaluation of the ETD taking into account the time steps of target motion model Tτ and of the camera of each UAV Tuz.The algorithmic description of the new multistepped planner is sketched in the pseudocode of Algorithm1. It starts computing the mission duration tend, the number of sections R into which s1:U0:tend will be divided, and the time lapse Tb required to update the “unnormalized belief” b~(τt). Next, it initializes the “unnormalized belief” b~(τ0) and the fitness criteria vector FCV0, which stores the different values of the criteria used by the GA to decide if a solution (proposal of UAV trajectories and sensor poses in s1:U0:t) is better than another. Finally, it sequentially optimizes each of the sections of the trajectory (i.e., each st:t+Q with t=Q(r-1) and r={1:R}) using the optimizer in optimGA(·), which will be described in detail in the following section.Algorithm 1:Multistepped GA-based planner.
Require:  sc⊳Structure with scenario models and times: b(τ0), Pτt∣τt-1, Pzt∣τt,sut, f(sut,aut), t1:Ustart, t1:Uend Require:  op⊳Structure with optimizer parameters Require:  s1:Utustart⊳Initial UAVs locations and sensor poses ∀u={1:U} Require:   Q⊳ Duration of each subsequence 1:tend=max⁡sc.t1:Uend⊳ Determine the duration of the mission 2:R=tend/Q⊳ Determine the number of sections required for the multi-stepped planner 3:tend=Q⋅R⊳ Round the duration of the mission 4:sc.Tb=obtainT(sc.T1:Uz,sc.Tτ)⊳ Determine minimum time step to update b~(τt). 5:t←0⊳ Set starting time 6:b~(τt)=sc.b(τ0)⊳ Initialize b~(τt) 7:FCVt=[0,0,0,0,0,0]⊳ Initialize fitness criteria vector 8:for r=1:R do⊳ For each subsequence 9:[s1:Ut:t+Q,FCVt+Q,b~(τt+Q)]←optimGA(sc,op,s1:Ut,b~(τt),t,Q,FCVt)⊳ Obtain best subsequence 10:t←t+Q⊳ Update time for the following section 11:end for 12:return  s1:U0:tend⊳ Complete UAV trajectory and sensor poses


## 3.2. GA Optimizer

This section details the most relevant aspects of the GA used to determine the best UAV trajectory and sensor pose for each section (st:t+Q with t=Q(r-1)) of the whole solution s0:tend. We start presenting the GA genotype, continue with the main GA steps and operators, and finish with the evaluation of each possible solution s0:tend.


### 3.2.1. GA Decision Variables

The decision variables directly manipulated by the genetic operators (i.e., the GA genotype) of the MTS planner are the subsequence of control actionsa1:UQ(r-1):Qr used to obtain, by means of the UAV motion model, the r-th subsection of the UAV trajectories s1:UQ(r-1):Qr (i.e., the GA phenotype or solution). Besides, the actions in aut are the commanded UAV heading θuc,t, speed vuc,t, and height huc,t and the camera elevation ηuc,t and azimuth αuc,t. Additionally, the number of decision variables is reduced by applying each action constantly during a time lapse Ta and by adapting the actions required for each UAV within each r-th subsection to the time steps where the UAV is engaged in the MTS mission. In other words, the GA algorithm only obtains the values of Q/Ta sets of actions (UAV heading, speed, and height; camera elevation and azimuth) for each UAV engaged in the mission during the r-th trajectory subsection.


### 3.2.2. GA General Description

The GA that sequentially optimizes each of the trajectory sectionss1:UQ(r-1):Qr based on the set of command actions (aut=[θuc,t,vuc,t,huc,t,ηuc,t,αuc,t]) and fitness criteria (which will be explained in the following section) maintains the operators, parameters, and structure of a standard GA and includes a new step to precalculate the number of actions required by each UAV.To present it in an algorithmic form, we extend the notation of the paper with left superindexes (1:J or best) and left subindexes (par and child). The first superindex is used to indicate that there is a population or group of J action sequences 1:Ja1:Ut:t+Q and fitness criteria vectors 1:JFCVt+Q. The second one is used to identify the best action sequence besta1:Ut:t+Q, UAV trajectories and sensor poses bests1:Ut:t+Q, and fitness criteria vector bestFCVt+Q within the population of 1:J solutions. And the left subindexes are incorporated to the previous notation when the algorithm needs to distinguish between the generic population elements (1:Ja1:Ut:t+Q, 1:JFCVt+Q) and those associated with the parents selected by the GA (par1:Ja1:Ut:t+Q) and with their children (child1:Ja1:Ut:t+Q, child1:JFCVt+Q).The GA in Algorithm2 performs the following steps. First, step 3 computes the number of actions of Ta duration required by each UAV in the current subsection of the trajectory, taking into account the current time step t and each UAV starting and ending engaging time (tustart and tuend). Step 4 initializes a random population of the J action sequences 1:Ja1:ut:t+Q required for the UAVs engaged in the mission, using a parametrizable uniform distribution. Steps 5-7 simulate the UAV trajectories and sensor poses from the sequence of actions and evaluate them according to the criteria presented in the following section. In step 9, the GA enters an optimization loop which ends when a fixed computation time has passed. Next, step 10 selects J parents of action sequences from the current population using binary tournament selection; step 11 creates the children by combining pairs of parents using a single-point crossover; step 12 slightly modifies the children with an incremental Gaussian noise mutation; steps 13-15 simulate and evaluate them; and step 16 selects the survival population from the previous and the mutated children using NSGA-II recombination [34]. Once the computation time has finished, step 19 identifies the best solution in the population; step 20 obtains the FCV, trajectories and sensor poses, and “unnormalized belief” associated with it; and step 21 returns them.Algorithm 2:optimGA: NSGA-II for MTS.
1:procedure  optimGA(sc, op, b~(τt),s1:Ut,t,Q,FCVt) 2:J←op.J⊳ Get population size 3:N1:U←obtainActionNumber(t,Q,Ta,sc.t1:Ustart,sc.t1:Ustart)⊳ Obtain number of actions for each UAV 4:1:Ja1:Ut:t+Q←initialize(J,N1:U,op.init)⊳ Generate J sequences with the required actions 5:for j=1:J do 6:jFCVt+Q←sim&Eval(sc,b~(τt),s1:Ut,t,Q,FCVt,ja1:Ut:t+Q)⊳ Simulate and evaluate 7:end for 8:M←1⊳ Initialize counter of the algorithm iterations (generations) 9:while  getCompTime()≤opt.CompTime  do⊳ Main optimization loop 10:par1:Ja1:Ut:t+Q←selectParents(1:Ja1:Ut:t+Q,1:JFCVt+Q)⊳ Select pair of parents (tournament selection) 11:child1:Ja1:Ut:t+Q←crossover(par1:Ja1:Ut:t+Q,op.xover)⊳ Perform crossover of parents (single point crossover) 12:child1:Ja1:Ut:t+Q←mutate(child1:Ja1:Ut:t+Q,op.mut)⊳ Mutate children (incremental gaussian mutation) 13:for j=1:J do 14:childjFCVt+Q←sim&Eval(sc,b~(τt),s1:Ut,t,Q,FCVt,childja1:Ut:t+Q)⊳ Simulate and evaluate 15:end for 16:[1:Ja1:Ut:t+Q,1:JFCVt+Q]←recombine(1:Ja1:Ut:t+Q,1:JFCVt+Q,child1:Ja1:Ut:t+Q,child1:JFCVt+Q) 17:M←M+1⊳ Increment the GA generation number 18:end while 19:besta1:Ut:t+Q←best(1:Ja1:Ut:t+Q,1:JFCVt+Q)⊳ Select best solution 20:[bestFCVt+Q,bests1:Ut:t+Q,bestb~(τt+Q)]←sim&Eval(sc,b~(τt),s1:Ut,t,Q,FCVt,besta1:Ut:t+Q)⊳ Sim. and evaluate it 21:returnbests1:Ut:t+Q,bestFCVt+Q,bestb~(τt+Q)⊳ Return best solution 22:end procedureIt is worth noting that there are several steps whose behavior can be configured using the information in theopt configuration variable. The initialization parameters opt.init in step 4 allow, for each UAV and type of action, (1) selecting different bounds for the uniform distribution and (2) enabling/disabling its optimization. These options let our MTS planner generate trajectories with some fixed behavior (e.g., fixed height or fixed camera orientations). The stop condition parameter opt.CompTime in step 9 allows fixing the computation time of the GA algorithm. The crossover parameter opt.xover in step 11 allows selecting the percentage of parents that undergo the single-point crossover or are directly copied as children. Finally, the mutation parameter opt.mut in step 12 allows selecting the Gaussian noise that is added to all variables in a1:Ut:t+Q and the Gaussian noise that is added to a few variables uniformly selected in 1:Ja1:Ut:t+Q with a probability 1/a1:Ut:t+Q, where a1:Ut:t+Q is the number of decision variables in a1:Ut:t+Q.Besides, there are three steps of the GA that use the fitness criteria variable1:JFCVt+Q to select solutions of the population to become pairs of parents (step 10), survive in the next generation (step 16), and become the best solution (step 19). Its behavior will be further detailed in the following section, after explaining the optimization criteria used in our MTS planner.


### 3.2.3. Simulation and Evaluation Process

In order to evaluate one of the solutions, subsequence of action commandsja1:Ut:t+Q, manipulated and proposed by the GA, first, the motion model in Figure 5 or in the appendix has to be simulated for each UAV, and, second, the fitness criteria have to be evaluated. The planner has also to take into account the UAVs mission engaging time [tustart,tuend]∀u=1:U and the possible differences among the time lapses of the problem (Tτ, Ta, Ts, Tuz). In the remainder of this section we describe the simulation process and the fitness criteria evaluation, and the steps of the sim&Eval function, presented in Algorithm 3 and used by the GA in Algorithm 2.Algorithm 3:sim&Eval: simulation and evaluation for MTS.
1:procedure  sim&Eval(sc,b~(τt),s1:Ut,t,Q,FCVt,a1:Ut:t+Q) 2:tinit←t 3:tlast←t+Q 4:⊳ Simulation of the proposed sequence of actions to obtain UAV trajectory and fuel consumption 5:for u=1:U do 6:    [sut:t+Q,fuelu]←sim(sc.f(⋅),sut,aut:t+Q,t,t+Q, ⊳ Simulate each UAV and sensor movement 7:sc.tustart,sc.tuend,sc.Ta,sc.Ts) 8:end for 9:⊳ Evaluation of the trajectory obtained from the proposed sequence of actions 10:FCVt[COL]←FCVt[COL]+∑u=1:U‍∑l=u+1:U‍collision(sc,sut:t+Q,slt:t+Q)⊳ Update collision criterion 11:FCVt[NFZ]←FCVt[NFZ]+∑u=1:U‍withinNFZ(sut:t+Q)⊳ Update fuel consumption 12:FCVt[SMOOTH]←FCVt[SMOOTH]+∑u=1:U‍smooth(aut:t+Q)⊳ Update smooth criterion 13:FCVt[FUEL]←FCVt[FUEL]+∑u=1:U‍fuelu⊳ Update fuel consumption 14:⊳ Simulation of the target and evaluation of the probabilistic criteria 15:t←tinit+sc.Tb 16:while  t≤tlast  do 17:if  motionTarget(t,sc.Tτ)  then 18:b~(τt-sc.Tτ)=b~(τt)⊳ Rename variable for the target transition 19:b~(τt)←∑τt-sc.Tτ∈GΩ‍P(τt∣τt-sc.Tτ)b~(τt-sc.Tτ)⊳ Target motion part of Equation (2) 20:end if 21:b~(τt)←∏u∈measuring(u,t,sc)‍P(zut=D¯∣τt,sut)b~(τt)⊳ Measurement update part of Equation (2) 22:Pnd(s1:U0:t)←∑τt∈GΩ‍b~(τt)⊳ Equation (3) 23:FCVt[ETD]←FCVt[ETD]+Pnd(s1:U0:t)∗sc.Tb⊳ Update ET, with equation (4) in sequential form 24:t←t+sc.Tb 25:end while 26:FCVt[MYOP]=∑τt∈GΩ‍∏u=∈engaged(u,t,sc)‍(1-λdistance(τt,sut))b~(τt)⊳ Calculate myopia reduction criterion 27:return  FCVt, s1:Utstart:tend, b~(τt) 28:end procedureThe simulation of each UAV is performed in step 6 of Algorithm3 with the sim(·) function, taking into account the Simulink or appendix model f(sut,aut), the initial UAV location sut, and the timely-equispaced set of actions in aut:t+Q. Moreover, during the simulation each set of actions aut proposed by the GA is considered constant during the time lapse Ta and the resulting UAV trajectory and sensor poses sut:t+Q are created with sampling simulation time of Ts. We use two different time lapses in the simulator (for Ts and Ta) to reduce the number of decision variables aut:t+Q manipulated by the planner, while permitting a high time resolution for evaluating certain aspects of the UAVs trajectories s1:Ut:t+Q. Finally, it is worth noting that sim(·) also obtains, from the UAV motion model, the fuel consumption of the UAV during the corresponding trajectory subsection and checks if the UAV is engaged in the mission in order to only simulate its displacement in the simulation time steps that simultaneously belong to [t,t+Q] and [tustart,tuend].The fitness criteria of the MTS planner consists in fulfilling two constraints (UAV collision and nonflying zones avoidance) and in optimizing four indexes (command smoothness, fuel consumption, expected detection time, and myopia reduction criterion), whose values are stored in the different elements of the Fitness Criterium VectorFCVt. Besides, in order to distinguish among the problem time periods Ta and Ts within the expressions of the evaluation criteria, we make [t:Δ:t+Q] stand for the time vector with elements equi-spaced by the selected Δ and k∈[t:Δ:t+Q]∧[a,b] stand for the time steps k within the set [t:Δ:t+Q] that fall within the limits [a,b]. Moreover, the tag used in Algorithm 3 within FCVt[tag] expressions permits the identification of the element of the vector FCVt associated with each criterion. Finally, the details of each criterion are presented below:(i)
UAV collision avoidance. This criterion counts the number of trajectory time instants at which the distance between the locations of all pairs of UAVs is smaller than the security collision distance dCOL. The criterion, evaluated in line 10 of Algorithm and detailed in (12), consists in checking the inter-UAV distances at time steps where both vehicles are engaged in the mission.(12)collisionsc,sut:t+Q,slt:t+Q=∑k∈t:Ts:t+Q∧sc.tustart,sc.tuend∧sc.tlstart,sc.tlendxuk-xlk2+yuk-ylk2+huk-hlk2≤sc.dCOL2(ii)
Nonflying zones (NFZ) avoidance. This criterion counts the number of trajectory time instants at which the locations of any UAV overfly a selected set of forbidden overflying cells of GΩ. Again, the criterion, evaluated in line 11 of Algorithm 3 and detailed in (13), considers locations at the time steps where the UAVs are in the mission.(13)withinNFZsc,sut:t+Q=∑k∈t:Ts:t+Q∧sc.tustart,sc.tuendinNFZxuk,yuk(iii)
Actions smoothness. This criterion, evaluated in line 12 of Algorithm 3, measures the smoothness of the commands for each UAV within aut:t+Q. It accounts, with (14) and (15), for the consecutive changes of the five types of commands only during the time interval where the UAVs are engaged in the mission. The magnitude of the accounted change, calculated with (16), makes the smooth criterion only account for those changes in action that imply a modification of sign in the signal with a value proportional to the square of the difference. Therefore, the smoothness criterion penalizes big changes of values and sign in the sequence of each type of command.(14)smoothsc,aut:t+Q=∑k∈t:Ta:t+Q-Ta∧sc.tustart,sc.tuendchangeAauk,auk+Ta(15)changeAaut,auk=chθuc,t,θuc,k+chhuc,t,huc,k+chvuc,t,vuc,k+chηuc,t,ηuc,k+chαuc,t,αuc,k(16)cha,b=a-b2sign⁡a-sign⁡b==2(iv)
Fuel consumption. Its value is directly returned by the UAV motion model and affected by the UAV speed.(v)
Expected time of detection. The main objective of the planner, described already in the general problem formulation section, is implemented in an iterative form within lines 16 and 25 of Algorithm 3. The process implements (2)-(4) after the following modifications. On the one hand, the expression in line 19 includes the target motion model transition operation for the time steps where the motion model has to be applied (i.e., when motionTarget(t,Tτ) is true, because the target is dynamic and the last target motion update was performed Tτ seconds ago). On the other hand, the “unnormalized belief” update expression in line 21 only includes the likelihood models of the UAVs at the time steps t where the UAV is engaged in the mission [tustart,tuend] and its sensor has to take a measure (according to its Tuz).(vi)
Myopia reduction criterion. The planner can create myopic solutions due to its multistepped behavior, which only considers the probability regions achievable from the starting position of each section of the trajectory. In other words, by optimizing only the previous criteria, the planner becomes greedy, ignoring during the optimization of each trajectory section the benefits of getting closer to a high probability area which is currently unreachable. To reduce this myopic behavior, the planner evaluates, in line 26 of Algorithm 3, the suitability of the final UAVs positions and sensor poses s1:UrQ of a trajectory section for the optimization of future trajectory sections by weighting the remaining unnormalized belief b~(τrQ) with a function that takes into account the distance of each cell τt in GΩ to the center of the camera footprint associated with the final location and camera pose of each UAV sut. Besides, this evaluation is only performed for the UAVs that are still engaged in the mission. In other words, the planner evaluates the suitability of the final camera footprints of the searching UAVs for optimizing the remaining sections of the trajectories.As we have previously mentioned, theFCVt values obtained after the simulation of each trajectory section are later used in 3 different steps of the GA. The way of proceeding in each step is the following:(i)
The binary tournament selection operator chooses randomly, with an uniform probability, two solutionsj of the population; among them it selects as a parent the one with lower constraint violation (i.e., with lower jFCVt[COL]+FjCVt[NFZ]), and when tied, the one that Pareto-dominates the other four optimization objectives (i.e., jFCVt[SMOOTH], jFCVt[FUEL], jFCVt[ETD], jFCVt[MYOP]).(ii)
The survival/recombination operator organizes the solutions in a similar way: first according to its constraint violation level and second according to the Pareto-dominance of the remaining four optimization objectives. Hence, as a whole, the GA tries to identify the solution that violates the constraints less times and, within them, those solutions that Pareto-dominate the other objective criteria.(iii)
The selection of a unique best solution is required at the final step of the GA to let the multistepped planner use its final position as the starting point of next trajectory section. To do this, the GA selects, among the solutions in the best Pareto-front of the last population obtained by the recombination operator, the one with smallest rounded value of myopia reduction criterion (first), smallest ETD (second, if tied in the first), smallest command smoothness (third, if tied in the previous), and small fuel consumption (fourth). The priority levels used for the selections process reflect the importance order of the objectives. Moreover, the myopia criterion is preferred to ETD, becauseFCVt[MYOP] accounts both for the probability gathered by the trajectory (through b~(τt)) and for the distance to all the cells with remaining “unnormalized belief”, while FCVt[ETD] only considers the reachable probability within each trajectory section s1:Ut:t+Q.


## 3.2.1. GA Decision Variables

The decision variables directly manipulated by the genetic operators (i.e., the GA genotype) of the MTS planner are the subsequence of control actionsa1:UQ(r-1):Qr used to obtain, by means of the UAV motion model, the r-th subsection of the UAV trajectories s1:UQ(r-1):Qr (i.e., the GA phenotype or solution). Besides, the actions in aut are the commanded UAV heading θuc,t, speed vuc,t, and height huc,t and the camera elevation ηuc,t and azimuth αuc,t. Additionally, the number of decision variables is reduced by applying each action constantly during a time lapse Ta and by adapting the actions required for each UAV within each r-th subsection to the time steps where the UAV is engaged in the MTS mission. In other words, the GA algorithm only obtains the values of Q/Ta sets of actions (UAV heading, speed, and height; camera elevation and azimuth) for each UAV engaged in the mission during the r-th trajectory subsection.


## 3.2.2. GA General Description

The GA that sequentially optimizes each of the trajectory sectionss1:UQ(r-1):Qr based on the set of command actions (aut=[θuc,t,vuc,t,huc,t,ηuc,t,αuc,t]) and fitness criteria (which will be explained in the following section) maintains the operators, parameters, and structure of a standard GA and includes a new step to precalculate the number of actions required by each UAV.To present it in an algorithmic form, we extend the notation of the paper with left superindexes (1:J or best) and left subindexes (par and child). The first superindex is used to indicate that there is a population or group of J action sequences 1:Ja1:Ut:t+Q and fitness criteria vectors 1:JFCVt+Q. The second one is used to identify the best action sequence besta1:Ut:t+Q, UAV trajectories and sensor poses bests1:Ut:t+Q, and fitness criteria vector bestFCVt+Q within the population of 1:J solutions. And the left subindexes are incorporated to the previous notation when the algorithm needs to distinguish between the generic population elements (1:Ja1:Ut:t+Q, 1:JFCVt+Q) and those associated with the parents selected by the GA (par1:Ja1:Ut:t+Q) and with their children (child1:Ja1:Ut:t+Q, child1:JFCVt+Q).The GA in Algorithm2 performs the following steps. First, step 3 computes the number of actions of Ta duration required by each UAV in the current subsection of the trajectory, taking into account the current time step t and each UAV starting and ending engaging time (tustart and tuend). Step 4 initializes a random population of the J action sequences 1:Ja1:ut:t+Q required for the UAVs engaged in the mission, using a parametrizable uniform distribution. Steps 5-7 simulate the UAV trajectories and sensor poses from the sequence of actions and evaluate them according to the criteria presented in the following section. In step 9, the GA enters an optimization loop which ends when a fixed computation time has passed. Next, step 10 selects J parents of action sequences from the current population using binary tournament selection; step 11 creates the children by combining pairs of parents using a single-point crossover; step 12 slightly modifies the children with an incremental Gaussian noise mutation; steps 13-15 simulate and evaluate them; and step 16 selects the survival population from the previous and the mutated children using NSGA-II recombination [34]. Once the computation time has finished, step 19 identifies the best solution in the population; step 20 obtains the FCV, trajectories and sensor poses, and “unnormalized belief” associated with it; and step 21 returns them.Algorithm 2:optimGA: NSGA-II for MTS.
1:procedure  optimGA(sc, op, b~(τt),s1:Ut,t,Q,FCVt) 2:J←op.J⊳ Get population size 3:N1:U←obtainActionNumber(t,Q,Ta,sc.t1:Ustart,sc.t1:Ustart)⊳ Obtain number of actions for each UAV 4:1:Ja1:Ut:t+Q←initialize(J,N1:U,op.init)⊳ Generate J sequences with the required actions 5:for j=1:J do 6:jFCVt+Q←sim&Eval(sc,b~(τt),s1:Ut,t,Q,FCVt,ja1:Ut:t+Q)⊳ Simulate and evaluate 7:end for 8:M←1⊳ Initialize counter of the algorithm iterations (generations) 9:while  getCompTime()≤opt.CompTime  do⊳ Main optimization loop 10:par1:Ja1:Ut:t+Q←selectParents(1:Ja1:Ut:t+Q,1:JFCVt+Q)⊳ Select pair of parents (tournament selection) 11:child1:Ja1:Ut:t+Q←crossover(par1:Ja1:Ut:t+Q,op.xover)⊳ Perform crossover of parents (single point crossover) 12:child1:Ja1:Ut:t+Q←mutate(child1:Ja1:Ut:t+Q,op.mut)⊳ Mutate children (incremental gaussian mutation) 13:for j=1:J do 14:childjFCVt+Q←sim&Eval(sc,b~(τt),s1:Ut,t,Q,FCVt,childja1:Ut:t+Q)⊳ Simulate and evaluate 15:end for 16:[1:Ja1:Ut:t+Q,1:JFCVt+Q]←recombine(1:Ja1:Ut:t+Q,1:JFCVt+Q,child1:Ja1:Ut:t+Q,child1:JFCVt+Q) 17:M←M+1⊳ Increment the GA generation number 18:end while 19:besta1:Ut:t+Q←best(1:Ja1:Ut:t+Q,1:JFCVt+Q)⊳ Select best solution 20:[bestFCVt+Q,bests1:Ut:t+Q,bestb~(τt+Q)]←sim&Eval(sc,b~(τt),s1:Ut,t,Q,FCVt,besta1:Ut:t+Q)⊳ Sim. and evaluate it 21:returnbests1:Ut:t+Q,bestFCVt+Q,bestb~(τt+Q)⊳ Return best solution 22:end procedureIt is worth noting that there are several steps whose behavior can be configured using the information in theopt configuration variable. The initialization parameters opt.init in step 4 allow, for each UAV and type of action, (1) selecting different bounds for the uniform distribution and (2) enabling/disabling its optimization. These options let our MTS planner generate trajectories with some fixed behavior (e.g., fixed height or fixed camera orientations). The stop condition parameter opt.CompTime in step 9 allows fixing the computation time of the GA algorithm. The crossover parameter opt.xover in step 11 allows selecting the percentage of parents that undergo the single-point crossover or are directly copied as children. Finally, the mutation parameter opt.mut in step 12 allows selecting the Gaussian noise that is added to all variables in a1:Ut:t+Q and the Gaussian noise that is added to a few variables uniformly selected in 1:Ja1:Ut:t+Q with a probability 1/a1:Ut:t+Q, where a1:Ut:t+Q is the number of decision variables in a1:Ut:t+Q.Besides, there are three steps of the GA that use the fitness criteria variable1:JFCVt+Q to select solutions of the population to become pairs of parents (step 10), survive in the next generation (step 16), and become the best solution (step 19). Its behavior will be further detailed in the following section, after explaining the optimization criteria used in our MTS planner.


## 3.2.3. Simulation and Evaluation Process

In order to evaluate one of the solutions, subsequence of action commandsja1:Ut:t+Q, manipulated and proposed by the GA, first, the motion model in Figure 5 or in the appendix has to be simulated for each UAV, and, second, the fitness criteria have to be evaluated. The planner has also to take into account the UAVs mission engaging time [tustart,tuend]∀u=1:U and the possible differences among the time lapses of the problem (Tτ, Ta, Ts, Tuz). In the remainder of this section we describe the simulation process and the fitness criteria evaluation, and the steps of the sim&Eval function, presented in Algorithm 3 and used by the GA in Algorithm 2.Algorithm 3:sim&Eval: simulation and evaluation for MTS.
1:procedure  sim&Eval(sc,b~(τt),s1:Ut,t,Q,FCVt,a1:Ut:t+Q) 2:tinit←t 3:tlast←t+Q 4:⊳ Simulation of the proposed sequence of actions to obtain UAV trajectory and fuel consumption 5:for u=1:U do 6:    [sut:t+Q,fuelu]←sim(sc.f(⋅),sut,aut:t+Q,t,t+Q, ⊳ Simulate each UAV and sensor movement 7:sc.tustart,sc.tuend,sc.Ta,sc.Ts) 8:end for 9:⊳ Evaluation of the trajectory obtained from the proposed sequence of actions 10:FCVt[COL]←FCVt[COL]+∑u=1:U‍∑l=u+1:U‍collision(sc,sut:t+Q,slt:t+Q)⊳ Update collision criterion 11:FCVt[NFZ]←FCVt[NFZ]+∑u=1:U‍withinNFZ(sut:t+Q)⊳ Update fuel consumption 12:FCVt[SMOOTH]←FCVt[SMOOTH]+∑u=1:U‍smooth(aut:t+Q)⊳ Update smooth criterion 13:FCVt[FUEL]←FCVt[FUEL]+∑u=1:U‍fuelu⊳ Update fuel consumption 14:⊳ Simulation of the target and evaluation of the probabilistic criteria 15:t←tinit+sc.Tb 16:while  t≤tlast  do 17:if  motionTarget(t,sc.Tτ)  then 18:b~(τt-sc.Tτ)=b~(τt)⊳ Rename variable for the target transition 19:b~(τt)←∑τt-sc.Tτ∈GΩ‍P(τt∣τt-sc.Tτ)b~(τt-sc.Tτ)⊳ Target motion part of Equation (2) 20:end if 21:b~(τt)←∏u∈measuring(u,t,sc)‍P(zut=D¯∣τt,sut)b~(τt)⊳ Measurement update part of Equation (2) 22:Pnd(s1:U0:t)←∑τt∈GΩ‍b~(τt)⊳ Equation (3) 23:FCVt[ETD]←FCVt[ETD]+Pnd(s1:U0:t)∗sc.Tb⊳ Update ET, with equation (4) in sequential form 24:t←t+sc.Tb 25:end while 26:FCVt[MYOP]=∑τt∈GΩ‍∏u=∈engaged(u,t,sc)‍(1-λdistance(τt,sut))b~(τt)⊳ Calculate myopia reduction criterion 27:return  FCVt, s1:Utstart:tend, b~(τt) 28:end procedureThe simulation of each UAV is performed in step 6 of Algorithm3 with the sim(·) function, taking into account the Simulink or appendix model f(sut,aut), the initial UAV location sut, and the timely-equispaced set of actions in aut:t+Q. Moreover, during the simulation each set of actions aut proposed by the GA is considered constant during the time lapse Ta and the resulting UAV trajectory and sensor poses sut:t+Q are created with sampling simulation time of Ts. We use two different time lapses in the simulator (for Ts and Ta) to reduce the number of decision variables aut:t+Q manipulated by the planner, while permitting a high time resolution for evaluating certain aspects of the UAVs trajectories s1:Ut:t+Q. Finally, it is worth noting that sim(·) also obtains, from the UAV motion model, the fuel consumption of the UAV during the corresponding trajectory subsection and checks if the UAV is engaged in the mission in order to only simulate its displacement in the simulation time steps that simultaneously belong to [t,t+Q] and [tustart,tuend].The fitness criteria of the MTS planner consists in fulfilling two constraints (UAV collision and nonflying zones avoidance) and in optimizing four indexes (command smoothness, fuel consumption, expected detection time, and myopia reduction criterion), whose values are stored in the different elements of the Fitness Criterium VectorFCVt. Besides, in order to distinguish among the problem time periods Ta and Ts within the expressions of the evaluation criteria, we make [t:Δ:t+Q] stand for the time vector with elements equi-spaced by the selected Δ and k∈[t:Δ:t+Q]∧[a,b] stand for the time steps k within the set [t:Δ:t+Q] that fall within the limits [a,b]. Moreover, the tag used in Algorithm 3 within FCVt[tag] expressions permits the identification of the element of the vector FCVt associated with each criterion. Finally, the details of each criterion are presented below:(i)
UAV collision avoidance. This criterion counts the number of trajectory time instants at which the distance between the locations of all pairs of UAVs is smaller than the security collision distance dCOL. The criterion, evaluated in line 10 of Algorithm and detailed in (12), consists in checking the inter-UAV distances at time steps where both vehicles are engaged in the mission.(12)collisionsc,sut:t+Q,slt:t+Q=∑k∈t:Ts:t+Q∧sc.tustart,sc.tuend∧sc.tlstart,sc.tlendxuk-xlk2+yuk-ylk2+huk-hlk2≤sc.dCOL2(ii)
Nonflying zones (NFZ) avoidance. This criterion counts the number of trajectory time instants at which the locations of any UAV overfly a selected set of forbidden overflying cells of GΩ. Again, the criterion, evaluated in line 11 of Algorithm 3 and detailed in (13), considers locations at the time steps where the UAVs are in the mission.(13)withinNFZsc,sut:t+Q=∑k∈t:Ts:t+Q∧sc.tustart,sc.tuendinNFZxuk,yuk(iii)
Actions smoothness. This criterion, evaluated in line 12 of Algorithm 3, measures the smoothness of the commands for each UAV within aut:t+Q. It accounts, with (14) and (15), for the consecutive changes of the five types of commands only during the time interval where the UAVs are engaged in the mission. The magnitude of the accounted change, calculated with (16), makes the smooth criterion only account for those changes in action that imply a modification of sign in the signal with a value proportional to the square of the difference. Therefore, the smoothness criterion penalizes big changes of values and sign in the sequence of each type of command.(14)smoothsc,aut:t+Q=∑k∈t:Ta:t+Q-Ta∧sc.tustart,sc.tuendchangeAauk,auk+Ta(15)changeAaut,auk=chθuc,t,θuc,k+chhuc,t,huc,k+chvuc,t,vuc,k+chηuc,t,ηuc,k+chαuc,t,αuc,k(16)cha,b=a-b2sign⁡a-sign⁡b==2(iv)
Fuel consumption. Its value is directly returned by the UAV motion model and affected by the UAV speed.(v)
Expected time of detection. The main objective of the planner, described already in the general problem formulation section, is implemented in an iterative form within lines 16 and 25 of Algorithm 3. The process implements (2)-(4) after the following modifications. On the one hand, the expression in line 19 includes the target motion model transition operation for the time steps where the motion model has to be applied (i.e., when motionTarget(t,Tτ) is true, because the target is dynamic and the last target motion update was performed Tτ seconds ago). On the other hand, the “unnormalized belief” update expression in line 21 only includes the likelihood models of the UAVs at the time steps t where the UAV is engaged in the mission [tustart,tuend] and its sensor has to take a measure (according to its Tuz).(vi)
Myopia reduction criterion. The planner can create myopic solutions due to its multistepped behavior, which only considers the probability regions achievable from the starting position of each section of the trajectory. In other words, by optimizing only the previous criteria, the planner becomes greedy, ignoring during the optimization of each trajectory section the benefits of getting closer to a high probability area which is currently unreachable. To reduce this myopic behavior, the planner evaluates, in line 26 of Algorithm 3, the suitability of the final UAVs positions and sensor poses s1:UrQ of a trajectory section for the optimization of future trajectory sections by weighting the remaining unnormalized belief b~(τrQ) with a function that takes into account the distance of each cell τt in GΩ to the center of the camera footprint associated with the final location and camera pose of each UAV sut. Besides, this evaluation is only performed for the UAVs that are still engaged in the mission. In other words, the planner evaluates the suitability of the final camera footprints of the searching UAVs for optimizing the remaining sections of the trajectories.As we have previously mentioned, theFCVt values obtained after the simulation of each trajectory section are later used in 3 different steps of the GA. The way of proceeding in each step is the following:(i)
The binary tournament selection operator chooses randomly, with an uniform probability, two solutionsj of the population; among them it selects as a parent the one with lower constraint violation (i.e., with lower jFCVt[COL]+FjCVt[NFZ]), and when tied, the one that Pareto-dominates the other four optimization objectives (i.e., jFCVt[SMOOTH], jFCVt[FUEL], jFCVt[ETD], jFCVt[MYOP]).(ii)
The survival/recombination operator organizes the solutions in a similar way: first according to its constraint violation level and second according to the Pareto-dominance of the remaining four optimization objectives. Hence, as a whole, the GA tries to identify the solution that violates the constraints less times and, within them, those solutions that Pareto-dominate the other objective criteria.(iii)
The selection of a unique best solution is required at the final step of the GA to let the multistepped planner use its final position as the starting point of next trajectory section. To do this, the GA selects, among the solutions in the best Pareto-front of the last population obtained by the recombination operator, the one with smallest rounded value of myopia reduction criterion (first), smallest ETD (second, if tied in the first), smallest command smoothness (third, if tied in the previous), and small fuel consumption (fourth). The priority levels used for the selections process reflect the importance order of the objectives. Moreover, the myopia criterion is preferred to ETD, becauseFCVt[MYOP] accounts both for the probability gathered by the trajectory (through b~(τt)) and for the distance to all the cells with remaining “unnormalized belief”, while FCVt[ETD] only considers the reachable probability within each trajectory section s1:Ut:t+Q.


## 3.3. Computational Cost of the MTS Planner

Finally, we calculate the computational complexity of our planner. To do this, we will consider that all UAVs start and end the mission at the same time steps (t1:Ustart=0, t1:Uend=tend) and that the values of quotients Q/Ts, Q/Ta, Q/Tτ, Q/Tuz, tend/Q are natural numbers.We start studying Algorithm3 and estimate the computation complexity of the UAV and sensor motion simulation (ccsim) and of the evaluation of each fitness criterion (cctag). The computational cost of the UAV and sensor motion simulation is proportional to the number of UAVs U and to the steps Q/Ts in which the simulator has to return UAV and sensor poses, ccsim=O(U·Q/Ts). The computation cost of the collision avoidance criterion is proportional to the number of possible UAV combinations U(U-1) and to the steps Q/Ts in which the UAV trajectory and sensor poses have been discretized. Hence, ccCOL=O(U2·Q/Ts). For the NFZ criterion, we can follow a similar process, but accounting only for the number of UAVs. Hence ccNFZ=O(U·Q/Ts). In the case of the control smoothness, we have to consider the number of UAVs U, the number of steps Q/Ta in which the action sequence is discretized, and the number of types of actions that are being optimized A. Hence ccSMOOTH=O(U·A·Q/Ta). As the fuel consumption is calculated during the simulation of the UAV and sensor motion, its computational complexity is already included in ccsim. To estimate the complexity associated with the ETD computation, we can observe that the cost associated with the target motion is proportional to Q/Tτ and (wxwy)2, and cost associated with the sensor update is proportional to (wxwy)∑u=1UQ/Tuz. The first cost is associated to the operations in the target transition ∑τt-Tτ∈GΩPτt∣τt-Tτb~(τt-Tτ), which can be calculated, a number of Q/Tτ times, as the multiplication of two matrixes of (wxwy)x(wxwy) and (wxwy)x1 size. The second cost is associated to the update of each sensor Pzut∣τt,sutb~(τt), which is calculated a number of Q/Tuz times and is proportional to wxwy. Hence, ccETD=O((wxwy)2Q/Tτ+(wxwy)∑u=1UQ/Tuz). Finally, as the computational cost of the myopia reduction criterion is proportional to U and to wxwy, its complexity ccMYOP=O(U·wxwy). The addition of the previous complexities, after simplifying the final expression, makes the expression of the computational complexity of the simulation and evaluation function become (17).(17)ccsim&Eval=OU2QTs+U·A·QTa+wxwy2QTτ+wxwyU+∑u=1UQTuzThe expression shows that the complexity of the simulation and evaluation function ccsim&Eval grows with the increment of the duration Q of each evaluation, the number of UAVs U, the number of decision variables manipulated by the planner A, and the number of cells of the grid of the search space wxwy. Besides, ccsim&Eval also grows as the time steps of the simulation Ts, of the actions Ta, of the target Tτ, and of the camera measurements Tuz are decremented. Hence, in order to have a reasonable computational cost, it is necessary to achieve a compromise among all those parameters. Besides, the quadratic dependency of ccsim&Eval with wxwy (i.e., with the size of the grid of the search space) is worth noting. This quadratic cost is especially relevant, because a moderate size search space (e.g., of 30 x 30 cells in the Areia Branca example) has already significantly bigger values of wx·wy (e.g., 900) than those used in other parameters (e.g., U=2, A=3, Q/Ts=50, Q/Ta=Q/Tuz=5, Q/Tτ=1).To account for the computational complexity of the planner, we still have to consider the computational cost of the GA in Algorithm2 and of the multistepped approach in Algorithm 1. For the GA, we compute the cost of the main loop of Algorithm 2. Considering a population size of J solutions, U·A·Q/Ta decision variables, and M iterations of the loop, the cost of the ccGAselparent=O(M·J), ccGAxover=ccGAmut=O(M·J·U·A·Q/Ta), ccGAsim&Eval=M·J·ccsim&Eval, and ccGArecombine=O(MJ2). This happens because computational cost of the tournament selection of each pair of parents is proportional to a constant, the cost of the single-point crossover and of the incremental mutation of a solution is proportional to the number of decision variables within it, and the cost of the recombination step of the population, which includes the sorting of the solutions, is proportional to J2. Hence, the overall GA complexity can be calculated as ccGA=M·Jccsim&Eval+O(U·A·Q/Ta+J). Finally, as the multistepped GA-based planner in Algorithm 1 calls the GA R=tend/Q times, its overall computation complexity becomes ccMS=tend/Q·M·Jccsim&Eval+O(U·A·Q/Ta+J), calculable with (18), obtained by substituting within it ccsim&Eval by (17).(18)ccMS=OtendQ·M·JJ+U2QTs+U·A·QTa+wxwy2QTτ+wxwyU+∑u=1UQTuzThe overall computational complexity shows the same dependencies observed for the computational complexity of the simulation and evaluation and a (multiplying) growing behavior proportional to the population sizeJ, the number of generations in GA M, and the number of subsections in which the trajectory is divided tend/Q. Besides, as Q appears in the denominator of tend/Q that multiplies the numerator of different quotients with the time steps, there is a part of the computational cost that can be considered independent of Q and of the number of sections in R=tend/Q in which the trajectories are divided, and dependent on the quotient of the duration of the mission tend and the different time steps (Ts, Ta, Tτ, and Tuz).Finally, it is worth noting that the complexity analysis, which is useful to consider the effects of selecting certain variables of the system in the computation cost, does not provide information about (1) the time required for evaluating the solutions or (2) the convergence speed of the multistepped planner and of each GA run, whose behavior is highly dependent on the probability models (initial belief, target motion model, and sensor detection function) of each scenario. To analyze the convergence speed, in the Results section we will analyze the temporal performance of our new planner (implemented in Matlab and using its Parallel Toolbox in the evaluation loops, lines 5-7 and 13-15, of the GA algorithm) over different scenarios to be able to determine, empirically, the time needed to optimize them (in a 2.5 GHz Intel Core i7 with 6 GB RAM PC with Windows 7). To provide some insight related to the time required for evaluating the solutions within the planner, we have measured the values of the different steps of the GA algorithm in the scenarios of the following section and observed that the mean computation time of each population step of the GA (loop in line 9 of Algorithm2) is 0.46 seconds for the first scenario, 0.49 seconds for the second one, and 1.05 seconds for the third one; that the evaluation of the solutions (sim&Eval function) takes 99% of the time; and that the UAV motion simulation takes 43% of the time of the sim&Eval function, while the rest of it takes 53%.


## 4. Related Work

This section compares the new planner presented in this paper with existing approaches and formulations of PTSP. Instead of performing an exhaustive state-of-the-art, we present a compilation of the works that have motivated this paper. In detail, we focus on two types of approaches: methods that optimize the ETD or other criteria where the visiting order of the cells is important [4–10], or more generic PTSP approaches that obtain either smooth trajectories of the UAVs [11–17, 20] or that also optimize the sensor orientation [21, 22]. Hence, we leave out of this comparison other works such as [18, 35–42] that do not share those properties. The comparison is summarized in Table 1, where the first column identifies each work and the remaining ones show the properties highlighted in the comparison, and where the works are grouped row-wise by the type of criteria they optimize (which is shown in the tenth column). In more detail:Target column shows which approaches deal withmultiple and/ormoving targets. Although including both properties permits handling a wider range of scenarios, it is worth noting that none of the ETD strategies [5–7, 9, 10], at the bottom of the table, can still handle multiple targets.UAVs column shows which planners are developed formultiple UAVs and their type ofdynamic model, grouped under the following labels:cell-to-cell indicates that the UAVs move from one cell of the belief to another cell according to a finite set of cardinal directions,waypoint indicates that the UAVs move according to segments defined by the waypoints selected by the planner,simplified indicates that the UAV motion model s˙ut=f(sut,aut) is linearized or consists in a constant speed steering model, andcomplex indicates that the planner works with a Simulink or differential motion model as ours. The benefits of multi-UAV approaches with complex UAV motion models are the possibilities of, on the one hand, observing different regions of the space simultaneously and, on the other hand, adapting the UAVs trajectories to the actual maneuvering capabilities of the UAVs.Sensor column shows thetype of sensor, its observationrange, and itsefficiency, as well as the planner capability ofmoving it. The algorithms under comparison use mainly downward-looking (DL)cameras andradars as well as generic sensors, whose behavior is identifiable by the labels in its observation range (a singlecell, a cell of thequadtree used to store the belief, and alimited orwide group of cells) and efficiency (ideal when the likelihood is always either 1 or 0, andFN andFP if it, respectively, considers false negatives and false positives). It is also worth noting thatwide range sensors are often modeled with decaying likelihood functions with the distance [9–12, 16], while the remaining sensors are often modeled with constant likelihood functions [4–6, 8, 13, 19, 21]. Besides, although camera sensors are used in several works [13, 14, 17, 21, 22], only a few approaches (including the one presented in this paper) reorient them to observe different regions of the search space easier.Optimization column shows thecriteria and theapproach used to improve them within each planner.On the one hand, the optimization criteria are used to organize the works row-wise in four groups: optimizing the probability of detection (maximizingPd or minimizing its counterpart Pnd), maximizing the information gain (IG), considering the trajectory distance or the time (by including a time discounting factor) in previous TPSP criteria, and minimizing the expected time of detection (ETD or the Monte Carlo Simulation Time of Detection MCSTD). Besides, in the table, the notation is extended with a right superindex in Pdk and Pndk in order to represent the corresponding probability of each (k) target (each with an initial belief and target motion model independent of the others). And to evaluate MCST for a given s0:tuend, several Monte Carlo simulations are run (each one consisting in sampling the target initial position τ0 from the initial belief b(τ0) and in measuring the time step of the trajectory where s0:tuend overflies the sample τ0) and the average value of the target detection times of all the simulation is calculated. In this regard it is worth noting that optimizing either Pd, Pnd, IG, or the uncertainty by itself is not the best option for MTS, where it is necessary to collect the maximum probability of detection as soon as possible. Besides, although including the distance (and hence making the approaches prefer solutions that optimize the PTSP criterion with shorter trajectories) or the time (by including a disconting factor that gives higher importance to the optimization of the criterion in the initial steps) in the timeless objectives alleviates the problem, those solutions [4, 20–22] are biased by the strategies and the parameters used to include them. Hence, they do not optimize ETD as well as the approaches in the last group [5] do.On the other hand, the optimization approach column shows the variability of methods used to handle the high complexity of PTSP, ranging fromgreedy orlocal approaches togradient orgraph based techniques,max-sum algorithms, Sequential Quadratic Programming (QSP), Multiarmed Bandit algorithms (MAB), Cross-Entropy Optimization (CEO), Ant Colony Optimization (ACO), or Genetic Algorithms (GA). In this regard it is worth noting that the selected approach is often associated with the underlying type of decision variables (discrete when the UAV dynamic model is cell-to-cell, continuous otherwise), the type and number of criteria functions, and the authors expertise. The main benefit of the selected GA optimizer, inherited from the planner in [10], comes from its capabilities, tested by many researchers over different real-world problems, to handle many-decision multicriteria optimization problems. In this regard, it is worth noting that the selection of a GA has been especially useful to add several constraints and optimization criteria to ETD and to incorporate the five types of decision variables (actions) in the new planner.Table 1
Works comparison.
Work Targets UAVs Sensor Optimization Multiple Moving Multiple Dynamic Model Type Range Efficiency Moving Criteria Approach [11] ✓ ✓ ✓ simplified radar (DL) wide FN  ∑k‍wkPdk SQP [12]  ✓  simplified radar (DL) wide FN  P d Greedy [13]   ✓ simplified camera (DL) limited FN  Pd Gradient [14]  ✓ ✓ simplified camera (DL) limited FN  Pd Max-sum [15]   ✓ simplified generic wide FN  P n d Gradient [16] ✓  ✓ simplified generic wide FN  ⋂ k = 1 Q P n d k Gradient  [17] ✓ ✓ ✓ waypoints camera (DL) limited FN/FP  IG Gradient [19] ✓  ✓ simplified camera (DL) limited FN/FP  Uncertainty Gradient  [20]    waypoints generic quadtree FN/FP  I G + d i s t . Greedy [21] ✓ ✓ ✓ waypoints camera limited FN ✓ P d / d i s t . Graph [22] ✓ ✓ ✓ simplified camera limited FN ✓ Discounted Uncertainty MAB [4]  ✓ ✓ cell-to-cell generic cell ideal  DiscountedPd CEO  [5]  ✓ ✓ cell-to-cell generic cell ideal  ETD BOA [6]  ✓ ✓ cell-to-cell generic cell ideal  ETD A C O M M A S [7]  ✓ ✓ cell-to-cell radar, sonar, camera  ETD CEO [8]    cell-to-cell generic cell ideal  MCSTD Greedy [9]  ✓ ✓ complex radar (DL) wide FN/FP  ETD with constraints ACOR [10]   ✓ complex radar (DL) wide FN/FP  ETD and others GA This work  ✓ ✓ complex camera limited FN ✓ ETD and others GAIf we analyze the table in detail, we can see that this work is the most general within the ETD group, as it allows the planner to obtain the best camera orientations for the best UAV trajectory, deal with moving target scenarios, and optimize multiple criteria. Moreover, the optimization approach in [8], which evaluates the expected time of detection via MCSTD, is directly targeted at initial Gaussian beliefs, while ours is targeted and tested against generic initial beliefs. Regarding the UAV motion model and the quantity of optimization criteria, this work is also more general than the ones classified in the other three groups. Its main drawback, inherited from existing ETD planners, is its single-target support. However, as we will show in one of the examples, this limitation can be shortened by using as initial belief the joint belief of the targets and making the UAVs continue the search mission until all the targets are found.Besides, a more detailed comparison with its closest MTS state-of-the-art planners ([9, 10]) shows the following facts. First, the new planner optimizes five action types (UAV orientation, speed, and height; sensor azimuth and elevation) while the planners in [9, 10] only optimize the UAV orientation. Moreover, as the good performance of the ACO-based planner in [9] is supported by the use of a heuristic for generating promising UAV headings, its extension is not straightforward as it will require building heuristics for generating promising UAV speeds and heights, and promising camera orientations. Second, the camera introduced in this paper challenges the planner as its limited range likelihood enhances the stepped and multimodal behavior of ETD(s1:U0:tend), while the soft decaying wide range likelihood of the radars in [9, 10] smooths ETD(s1:U0:tend), helping those planners. Finally, although the new planner and the one in [10] share the multistepped GA-optimization support, UAV motion model, constraints, and objective functions, the new planner (and evaluation functions) extends the capabilities of the planner in [10] by allowing obtaining solutions (UAV trajectories plus sensor poses) for scenarios with dynamic targets, using five different types of actions, with different UAV motion model parameterizations, and with different UAV mission engaging and camera-capture times.An additional contribution of this work to MTS search, not included in the table due to space limitation, is the methodology followed to build realistic scenarios for testing the planner. Moreover, the next section details how to set up the scenario correctly, as there is a strong influence of some parameters (e.g., the size of the target, the camera FOV, and the height and size of the search region) on others (such as the recommended UAV flying altitude, camera likelihood, and search cell size) and, hence, on the planner performance.Finally, for comparison purposes with existing approaches, we can set up the new planner with the capabilities of the one presented in [10], as it is a direct extension of it and the new planner allows selecting which decision variables have to be optimized. We do not compare against the other works of the ETD group, as the planner in [9] is not capable of optimizing all the criteria used in this work, and the remaining ones [5–8] are designed for UAV motion models that move from cell to cell and for optimizing a single criterion.


## 5. Results

In this section the new planner is tested over three real-world inspired scenarios. We first explain their common properties and the analysis setup, and afterwards we detail the individual characteristics of each scenario and the solutions obtained by our planner under two different configurations: optimizing the UAVs trajectories and sensor poses, and only optimizing the UAV orientation (as in [10]). Therefore, this second configuration is used as a baseline to compare against and observe the benefits of this work proposal of optimizing multiple types of actions (including the sensor poses).


### 5.1. Common Scenario and Planner Setup

On the one hand, to study the benefits of optimizing multiple decision variables, we make the GA of our planner optimize: (1) only the commanded UAV headingθuc,t (as in the previous planner in [10]) or (2) the commanded UAV heading θuc,t, UAV speed suc,t, and camera azimuth αuc,t. For shortening, hereafter the first configuration will be labeled as GAθ and the second as GAθ,v,α. Besides, in this work, we prefix the values of the remaining decision variables (commanded UAV height huc,t and commanded camera elevation ηuc,t) to facilitate the analysis of the results. Finally, for comparison purposes, in GAθ we prefix the UAV speed and camera azimuth (i.e., the two extra decision variables of GAθ,v,α) to the mean value of the permitted ranges in GAθ,v,α.On the other hand, and due to the stochastic nature of a GA, we run each planner configuration 20 times and represent the evolution (against the computation time) of the mean and variance of theETD(s0:tend) obtained at different runs. Note that we only show the ETD as it is in fact the most relevant objective for this planner. Besides, as ETD(s0:tend) is not available until the last section of the trajectory is being optimized, for the computation times where the multistepped planner has not yet started the final section optimization, we approximate the value of ETD(s0:tend) with the following expression, which accounts for the ETD up to the section that is being optimized plus the worst case (associated with the case in which the UAVs do not collect more beliefs after s1:UrQ) for the remaining time (tend-(rQ+Tτ)).(19)ETDas1:U0:tend≈ETDs1:U0:rQ+tend-rQ+Tτ·Pnds1:U0:rQIt is worth noting that the approximation in (19) induces the stepped behavior of the ETD curves presented in the paper in Figures 7(d), 9(d), and 10(d). This happens because the contribution of the worst case term in the ETDa, which is associated with the belief that can/has not yet been collected by the UAVs, is usually decremented after a new section of the trajectory is being optimized, because as soon as the new probability of nondetection Pnd(s1:U0:(r+1)Q) is lower than the one of the previous section Pnd(s1:U0:rQ), there is at least a decrement of (tend-(Q(r+1)+Tτ))Pnd(s1:U0:rQ)-Pnd(s1:U0:(r+1)Q) in ETDa. Moreover, in the initial sections, where r is lower, for a similar reduction in the nondetection probability, the decrement in the represented value of the ETDa is often bigger.The results of the planner for each scenario are, respectively, schematized in Figures7, 9, and 10. In all of them, the graphic at the left-top (a) represents the initial b(τ0) and UAVs initial locations (colored airplanes), the graphic underneath (d) shows the evolution of the ETDa(s1:U0:tend) mean and standard deviation, and the second and third columns display additional information of a representative solution obtained by each planner configuration (only heading in the second column, and heading + speed + azimuth in the third one). In more detail, the top row graphics of the second and third column show the b~(τtend), ETD(s1:U0:tend), and the probability of detection Pd(s1:U0:tend) of a representative solutions of the corresponding planner configuration, while the graphics at the bottom row also represent the UAVs trajectories (with the red and blue lines), the camera footprints (with the shadowed red and blue polygons), and the camera aiming point without terrain elevation (by the end of the yellow line that starts at the UAV location of each measurement).Additionally, the elevation information used to construct the initial beliefs for each scenario is available in [43]. Besides, the wind speed uwt=0 in all the scenarios and all the cameras has a resolution of 2160x3840 pixels and a fixed elevation angle of 45°. Moreover, the camera FOV and UAV altitude huc,t are precalculated and maintained constant for each mission to ensure a high detection likelihood during the whole search. To achieve this, the FOV and huc,t are adjusted to make acceptable the value of TTPF for the worst possible case, which occurs at the highest distance between the target and camera and, hence, at the lowest terrain elevation. Additionally, the size of the cells in GΩ is calculated for each mission to ensure that the footprint covers more than one cell in the worst case (at the highest terrain elevation). Finally, all missions share the frequency of the camera captures and command changes (Tuz=Ta=10s), the simulation step (Ts=1s), the duration of the trajectory sections (Q=50s), and the collision security distance (dCOL=100m).Lastly, we inherit several values of the planner parameters from [10]. The population size J=50, the probability of crossover is 0.8, the probability of mutation 1/(∑u=1:UNu), and the high-lower additive mutation standard deviation is 1.0-0.1. Besides, the myopia parameter λ=0.95.


### 5.2. Looking for a Static Off-Road Vehicle in the Mountains

The target of the first scenario (proposed by Airbus) is an off-road vehicle broken down in Gador mountains (Spain), within a search area of 10x10km2 whose elevation varies, according to the elevation map in Figure 6(a), from 626 to 2037m. The elevation probability layer, represented at the bottom of Figure 6(b), is obtained by dividing the elevation range into four fractions and assigning lower probabilities to higher elevations. The intelligence layer, displayed at the top of Figure 6(b), consists in two Gaussians centered at the valleys oftener visited by the vehicle. The initial belief, displayed in Figure 6(c), is obtained merging both layers with w1=15 and w2=1, since this weights combination makes the two Gaussians be only 1/16 (i.e., 6.25%) of b(τ0). Besides, the only UAV engaged in this MTS operation, from t1start=0 to t1end=450s, starts the mission at the right top corner of the scenario and can fly at h1t∈[500,3000]m (with respect to the sea level) and v1t∈[30,60]m/s, with camera azimuth α1t∈[-30°,30°]. We select a camera FOV=25° and UAV h1c,t=3000m to ensure that the sensor likelihood varies, according to the elevations in Ω and target size (Starget≈5m), from 0.90 to 0.99. Besides, as the estimated smallest complete footprint, according to the UAV constant height and maximum terrain elevation, is 590x482m2, we divide Ω into a grid of 30x30 cells, each of 333x333m2.Figure 6
Scenario for searching for a lost vehicle in Gador mountains (Almeria, Spain).
(a)
Elevation (m) (b)
Probability layers (c)
Scenario (b(τ0), s10)Figure 7
Results for searching for a lost vehicle in Gador mountains (Almeria, Spain).
(a)
Scenario (b(τ0), s10) (b)
b ~ ( τ t e n d ) of GAθ (c)
b ~ ( τ t e n d ) of GAθ,v,α (d)
E T D a (e)
Solution ofGAθ (f)
Solution ofGAθ,v,αThe results of running both configurations of the planner for this scenario are presented in Figure7. The graphic of the ETD evolution in Figure 7(d), which represents the approximated ETDa up to the optimization of the last section of the trajectory, clearly shows the stepped behavior associated with the inclusion of new trajectory sections s1:UQ(r-1):Qr in the optimization of the planner. It is worth noting that the duration of the ETD steps can be easily clipped reducing the computational time of the GA stop condition, as each GA run within the multistepped planner converges to its ETD plateau long before its designated time has finished. Besides, the ETD evolution graphic also shows that the GAθ,v,α produces better results (with lower ETD) during all the optimization processes than GAθ.Furthermore, if we analyze the additional information of the representative solutions of both configurations, presented in Figures7(b), 7(c), 7(e), and 7(f), we can observe that GAθ,v,α permits the UAV + camera to collect more Pd(s1:U0:tend) than GAθ (51% vs. 43%), because the UAV adapts its speed to collect slowly the two Gaussians and arrive quicker to the southeast region of Ω. Besides, the footprints obtained with GAθ,v,α are more scattered than those obtained by GAθ, due to the possibility of changing the camera azimuth. Finally, it is worth noting how the shape of the footprints changes and gets closer to the camera aiming point at zero elevation in the valleys (at the north and at the south west) than in the mountains.


### 5.3. Looking for Survivors after an Earthquake

In the second scenario, inspired by the earthquake of Haiti in 2010, we search survivors at Carrefour (Port-au-Prince) in a rectangular area of 10x7km2, whose elevation varies, according to the elevation map in Figure 8(a), from sea level to 875m. To use our single-target search planner for a multitarget scenario, we only have to make the initial belief of our scenario the joint-probability of all the survivors and make the UAVs perform the search until all the targets are found [16]. The elevation probability layer, represented at the bottom of Figure 8(b), is defined by dividing the elevation range into three fractions and assigning 0 probability of survivors presence at the sea, high probability at the coast, and lower probability to elevations bigger than 400m. The intelligence layer, represented at the top of Figure 8(b), indicates the probability of survivors based on building damage assessment. This information, available in [44] and used in [45] for another search problem, was obtained by several public and private entities, comparing previous and a posteriori aerial imagery of the earthquake area. The initial belief, represented in Figure 8(c), is obtained merging both layers with w1=1 and w2=1, in order to make the intelligence information represent 1/2 (i.e., 50%) of b(τ0). The first UAV (in red in Figure 8(c)) arrives at the search area Ω from the East at t1start=0s, while the second UAV (in blue) does it from the south at t2start=50s. Both UAVs remain in the mission up to t1:2end=350s and can fly at h1:2t∈[500,3000]m and v1:2t∈[30,60]m/s, with camera azimuth α1:2t∈[-30°,30°]. However, due to the survivors size (Starget≈2m), both UAVs are forced to fly at lower altitudes h1:2c,t=1500m with a smaller FOV=20° than in the first scenario, in order to ensure an acceptable camera likelihood that varies from 0.82 to 0.99 (respectively, at the lowest and highest elevations of the search region). Besides, as the estimated smallest complete footprint in this scenario is 383x313m2, the scenario is divided into a grid of 40x28 cells, each of them of 250x250m2.Figure 8
Scenario of searching for survivors after an earthquake in Carrefour (Port-au-Prince, Haiti).
(a)
Elevation (m) (b)
Probability layers (c)
Scenario (b(τ0), s1:20)Figure 9
Results of searching for survivors after an earthquake in Carrefour (Port-au-Prince, Haiti).
(a)
Scenario (b(τ0), s1:20) (b)
b ~ ( τ t e n d ) of GAθ (c)
b ~ ( τ t e n d ) of GAθ,v,α (d)
E T D a (e)
Solution ofGAθ (f)
Solution ofGAθ,v,αFigure 10
Results of searching for a drifting boat (Areia Branca, Brazil).
(a)
Scenario (b(τ0), s1:20, NFZ) (b)
b ~ ( τ t e n d ) of GAθ (c)
b ~ ( τ t e n d ) of GAθ,v,α (d)
E T D a (e)
Solution ofGAθ (f)
Solution ofGAθ,v,αThe results of running both configurations of the planner for this scenario are presented in Figure9. The ETD evolution graph in Figure 9(d) shows again that the proposed GAθ,v,α outperforms GAθ regarding the ETD. This happens because both UAVs of GAθ,v,α can arrive quicker to the high probability areas (in the damaged buildings). It is worth noting also that the probabilities of detection Pd(s1:U0:tend) at the end of the mission of two representative solutions of each GA are similar (54% in GAθ vs. 55% in GAθ,v,α), because both solutions manage to get the same belief during the development of the mission. In more detail, the analysis of b~(τtend) and of the trajectories in the second and third column of Figure 9 shows that the solution by GAθ,v,α collects quicker the concentrated intelligence belief than the elevation belief, while the solution by GAθ gets quicker the elevation belief than part of the intelligence belief. However, the similarity in the values of the final probability of detection does not imply that the solutions are equivalent for our planner, as its purpose is to minimize the expected time of detection (which is better for GAθ,v,α than for GAθ) and not the probability of detection.Another interesting result of this scenario is the higher variance in the ETD ofGAθ,v,α. This is caused by the uniformity of elevation probability layer, which hardens the scenario for multistepped planners based on generic heuristics as ours. The effect is bigger in the ETD of GAθ,v,α than in the ETD of GAθ, since its capability to manipulate more decision variables allows the planner to obtain more solutions that are similar to the myopic criterium and different from the ETD perspective.


### 5.4. Looking for a Drifting Boat in the Sea

The target of the last scenario is a drifting boat (Starget≈5m) close to the coast (Areia Branca, Brazil) within a square area of 15x15km2, whose elevation varies from sea level to 147m. The construction of its initial belief and target motion models has already been described in the Minimum Time Search Definition section.Besides in this scenario, there is a NFZ (e.g., restricted military area) over the coast cells marked in red in Figure10(a), where the green arrows summarize the tendency of the target motion model. Both UAVs enter the search area Ω from the west at t1:2start=0 and stay in the mission up to t1:Uend=400s. Their flying capabilities are different: the first UAV (in red) can do it at h1t∈[500,1500]m and v1t∈[60,90]m/s, while the second (in blue) can do it at h2t∈[500,2000]m and v2t∈[30,60]m/s. The camera azimuth of both α1:2t∈[-90°,90°]. Besides, both UAVs are required to fly at their highest altitude with a FOV=30° to ensure that the sensor likelihood of the first UAV varies from 0.95 to 0.97 and of the second one from 0.89 to 0.91. Additionally, as the estimated smallest complete footprint is 990x821m2, we divide the search area Ω into a grid of 30x30 cells of 500x500m2. Finally, as the drifting boat is a dynamic target, we also have to select its motion lapse. Taking into account the size of the cells, we make Tτ=50s with the purpose of letting the target advance at approximately 10m/s.The results of running both configurations of the planner for this scenario are presented in Figure10. Again, the ETD evolution graph in Figure 10(d) shows that GAθ,v,α produces overall better ETD results than GAθ, at the expense of a bigger variability of its ETD. The initial closer mean values of the ETD results in GAθ,v,α and GAθ are due to the fact that the b(τ0) of this scenario and the initial locations of the UAVs let their cameras initially observe areas with the same probability independently of the UAVs speed and cameras azimuth. The benefit on ETD comes later, when the UAVs with changing velocity can arrive quicker to the high probability areas.For this case, thePd(s1:U0:tend) of the representative solution of GAθ,v,α is also higher than its counterpart in GAθ (74% vs. 71%). That is, the planner with GAθ,v,α lets the UAVs collect quicker more probability. Several benefits of changing vuc,t and αuc,t are clearly observed in the solutions displayed in the second and third columns of Figure 10. For GAθ,v,α configuration, not only does the red UAV arrive sooner to the moving Gaussian, but also the blue UAV orients the camera to cover the coast at its right while flying straightforward, while in GAθ the red UAV follows the track left by the Gaussian and the blue UAV has to change its orientation several times to be able to collect the belief on the coast. Hence, GAθ,v,α obtains safer straighter flights, where the UAV headings changes proposed by GAθ are substituted by camera reorientations. It is also worth noting that in both solutions the red footprint is smaller than the blue, because the first UAV flies at lower altitude than the second.


## 5.1. Common Scenario and Planner Setup

On the one hand, to study the benefits of optimizing multiple decision variables, we make the GA of our planner optimize: (1) only the commanded UAV headingθuc,t (as in the previous planner in [10]) or (2) the commanded UAV heading θuc,t, UAV speed suc,t, and camera azimuth αuc,t. For shortening, hereafter the first configuration will be labeled as GAθ and the second as GAθ,v,α. Besides, in this work, we prefix the values of the remaining decision variables (commanded UAV height huc,t and commanded camera elevation ηuc,t) to facilitate the analysis of the results. Finally, for comparison purposes, in GAθ we prefix the UAV speed and camera azimuth (i.e., the two extra decision variables of GAθ,v,α) to the mean value of the permitted ranges in GAθ,v,α.On the other hand, and due to the stochastic nature of a GA, we run each planner configuration 20 times and represent the evolution (against the computation time) of the mean and variance of theETD(s0:tend) obtained at different runs. Note that we only show the ETD as it is in fact the most relevant objective for this planner. Besides, as ETD(s0:tend) is not available until the last section of the trajectory is being optimized, for the computation times where the multistepped planner has not yet started the final section optimization, we approximate the value of ETD(s0:tend) with the following expression, which accounts for the ETD up to the section that is being optimized plus the worst case (associated with the case in which the UAVs do not collect more beliefs after s1:UrQ) for the remaining time (tend-(rQ+Tτ)).(19)ETDas1:U0:tend≈ETDs1:U0:rQ+tend-rQ+Tτ·Pnds1:U0:rQIt is worth noting that the approximation in (19) induces the stepped behavior of the ETD curves presented in the paper in Figures 7(d), 9(d), and 10(d). This happens because the contribution of the worst case term in the ETDa, which is associated with the belief that can/has not yet been collected by the UAVs, is usually decremented after a new section of the trajectory is being optimized, because as soon as the new probability of nondetection Pnd(s1:U0:(r+1)Q) is lower than the one of the previous section Pnd(s1:U0:rQ), there is at least a decrement of (tend-(Q(r+1)+Tτ))Pnd(s1:U0:rQ)-Pnd(s1:U0:(r+1)Q) in ETDa. Moreover, in the initial sections, where r is lower, for a similar reduction in the nondetection probability, the decrement in the represented value of the ETDa is often bigger.The results of the planner for each scenario are, respectively, schematized in Figures7, 9, and 10. In all of them, the graphic at the left-top (a) represents the initial b(τ0) and UAVs initial locations (colored airplanes), the graphic underneath (d) shows the evolution of the ETDa(s1:U0:tend) mean and standard deviation, and the second and third columns display additional information of a representative solution obtained by each planner configuration (only heading in the second column, and heading + speed + azimuth in the third one). In more detail, the top row graphics of the second and third column show the b~(τtend), ETD(s1:U0:tend), and the probability of detection Pd(s1:U0:tend) of a representative solutions of the corresponding planner configuration, while the graphics at the bottom row also represent the UAVs trajectories (with the red and blue lines), the camera footprints (with the shadowed red and blue polygons), and the camera aiming point without terrain elevation (by the end of the yellow line that starts at the UAV location of each measurement).Additionally, the elevation information used to construct the initial beliefs for each scenario is available in [43]. Besides, the wind speed uwt=0 in all the scenarios and all the cameras has a resolution of 2160x3840 pixels and a fixed elevation angle of 45°. Moreover, the camera FOV and UAV altitude huc,t are precalculated and maintained constant for each mission to ensure a high detection likelihood during the whole search. To achieve this, the FOV and huc,t are adjusted to make acceptable the value of TTPF for the worst possible case, which occurs at the highest distance between the target and camera and, hence, at the lowest terrain elevation. Additionally, the size of the cells in GΩ is calculated for each mission to ensure that the footprint covers more than one cell in the worst case (at the highest terrain elevation). Finally, all missions share the frequency of the camera captures and command changes (Tuz=Ta=10s), the simulation step (Ts=1s), the duration of the trajectory sections (Q=50s), and the collision security distance (dCOL=100m).Lastly, we inherit several values of the planner parameters from [10]. The population size J=50, the probability of crossover is 0.8, the probability of mutation 1/(∑u=1:UNu), and the high-lower additive mutation standard deviation is 1.0-0.1. Besides, the myopia parameter λ=0.95.


## 5.2. Looking for a Static Off-Road Vehicle in the Mountains

The target of the first scenario (proposed by Airbus) is an off-road vehicle broken down in Gador mountains (Spain), within a search area of 10x10km2 whose elevation varies, according to the elevation map in Figure 6(a), from 626 to 2037m. The elevation probability layer, represented at the bottom of Figure 6(b), is obtained by dividing the elevation range into four fractions and assigning lower probabilities to higher elevations. The intelligence layer, displayed at the top of Figure 6(b), consists in two Gaussians centered at the valleys oftener visited by the vehicle. The initial belief, displayed in Figure 6(c), is obtained merging both layers with w1=15 and w2=1, since this weights combination makes the two Gaussians be only 1/16 (i.e., 6.25%) of b(τ0). Besides, the only UAV engaged in this MTS operation, from t1start=0 to t1end=450s, starts the mission at the right top corner of the scenario and can fly at h1t∈[500,3000]m (with respect to the sea level) and v1t∈[30,60]m/s, with camera azimuth α1t∈[-30°,30°]. We select a camera FOV=25° and UAV h1c,t=3000m to ensure that the sensor likelihood varies, according to the elevations in Ω and target size (Starget≈5m), from 0.90 to 0.99. Besides, as the estimated smallest complete footprint, according to the UAV constant height and maximum terrain elevation, is 590x482m2, we divide Ω into a grid of 30x30 cells, each of 333x333m2.Figure 6
Scenario for searching for a lost vehicle in Gador mountains (Almeria, Spain).
(a)
Elevation (m) (b)
Probability layers (c)
Scenario (b(τ0), s10)Figure 7
Results for searching for a lost vehicle in Gador mountains (Almeria, Spain).
(a)
Scenario (b(τ0), s10) (b)
b ~ ( τ t e n d ) of GAθ (c)
b ~ ( τ t e n d ) of GAθ,v,α (d)
E T D a (e)
Solution ofGAθ (f)
Solution ofGAθ,v,αThe results of running both configurations of the planner for this scenario are presented in Figure7. The graphic of the ETD evolution in Figure 7(d), which represents the approximated ETDa up to the optimization of the last section of the trajectory, clearly shows the stepped behavior associated with the inclusion of new trajectory sections s1:UQ(r-1):Qr in the optimization of the planner. It is worth noting that the duration of the ETD steps can be easily clipped reducing the computational time of the GA stop condition, as each GA run within the multistepped planner converges to its ETD plateau long before its designated time has finished. Besides, the ETD evolution graphic also shows that the GAθ,v,α produces better results (with lower ETD) during all the optimization processes than GAθ.Furthermore, if we analyze the additional information of the representative solutions of both configurations, presented in Figures7(b), 7(c), 7(e), and 7(f), we can observe that GAθ,v,α permits the UAV + camera to collect more Pd(s1:U0:tend) than GAθ (51% vs. 43%), because the UAV adapts its speed to collect slowly the two Gaussians and arrive quicker to the southeast region of Ω. Besides, the footprints obtained with GAθ,v,α are more scattered than those obtained by GAθ, due to the possibility of changing the camera azimuth. Finally, it is worth noting how the shape of the footprints changes and gets closer to the camera aiming point at zero elevation in the valleys (at the north and at the south west) than in the mountains.


## 5.3. Looking for Survivors after an Earthquake

In the second scenario, inspired by the earthquake of Haiti in 2010, we search survivors at Carrefour (Port-au-Prince) in a rectangular area of 10x7km2, whose elevation varies, according to the elevation map in Figure 8(a), from sea level to 875m. To use our single-target search planner for a multitarget scenario, we only have to make the initial belief of our scenario the joint-probability of all the survivors and make the UAVs perform the search until all the targets are found [16]. The elevation probability layer, represented at the bottom of Figure 8(b), is defined by dividing the elevation range into three fractions and assigning 0 probability of survivors presence at the sea, high probability at the coast, and lower probability to elevations bigger than 400m. The intelligence layer, represented at the top of Figure 8(b), indicates the probability of survivors based on building damage assessment. This information, available in [44] and used in [45] for another search problem, was obtained by several public and private entities, comparing previous and a posteriori aerial imagery of the earthquake area. The initial belief, represented in Figure 8(c), is obtained merging both layers with w1=1 and w2=1, in order to make the intelligence information represent 1/2 (i.e., 50%) of b(τ0). The first UAV (in red in Figure 8(c)) arrives at the search area Ω from the East at t1start=0s, while the second UAV (in blue) does it from the south at t2start=50s. Both UAVs remain in the mission up to t1:2end=350s and can fly at h1:2t∈[500,3000]m and v1:2t∈[30,60]m/s, with camera azimuth α1:2t∈[-30°,30°]. However, due to the survivors size (Starget≈2m), both UAVs are forced to fly at lower altitudes h1:2c,t=1500m with a smaller FOV=20° than in the first scenario, in order to ensure an acceptable camera likelihood that varies from 0.82 to 0.99 (respectively, at the lowest and highest elevations of the search region). Besides, as the estimated smallest complete footprint in this scenario is 383x313m2, the scenario is divided into a grid of 40x28 cells, each of them of 250x250m2.Figure 8
Scenario of searching for survivors after an earthquake in Carrefour (Port-au-Prince, Haiti).
(a)
Elevation (m) (b)
Probability layers (c)
Scenario (b(τ0), s1:20)Figure 9
Results of searching for survivors after an earthquake in Carrefour (Port-au-Prince, Haiti).
(a)
Scenario (b(τ0), s1:20) (b)
b ~ ( τ t e n d ) of GAθ (c)
b ~ ( τ t e n d ) of GAθ,v,α (d)
E T D a (e)
Solution ofGAθ (f)
Solution ofGAθ,v,αFigure 10
Results of searching for a drifting boat (Areia Branca, Brazil).
(a)
Scenario (b(τ0), s1:20, NFZ) (b)
b ~ ( τ t e n d ) of GAθ (c)
b ~ ( τ t e n d ) of GAθ,v,α (d)
E T D a (e)
Solution ofGAθ (f)
Solution ofGAθ,v,αThe results of running both configurations of the planner for this scenario are presented in Figure9. The ETD evolution graph in Figure 9(d) shows again that the proposed GAθ,v,α outperforms GAθ regarding the ETD. This happens because both UAVs of GAθ,v,α can arrive quicker to the high probability areas (in the damaged buildings). It is worth noting also that the probabilities of detection Pd(s1:U0:tend) at the end of the mission of two representative solutions of each GA are similar (54% in GAθ vs. 55% in GAθ,v,α), because both solutions manage to get the same belief during the development of the mission. In more detail, the analysis of b~(τtend) and of the trajectories in the second and third column of Figure 9 shows that the solution by GAθ,v,α collects quicker the concentrated intelligence belief than the elevation belief, while the solution by GAθ gets quicker the elevation belief than part of the intelligence belief. However, the similarity in the values of the final probability of detection does not imply that the solutions are equivalent for our planner, as its purpose is to minimize the expected time of detection (which is better for GAθ,v,α than for GAθ) and not the probability of detection.Another interesting result of this scenario is the higher variance in the ETD ofGAθ,v,α. This is caused by the uniformity of elevation probability layer, which hardens the scenario for multistepped planners based on generic heuristics as ours. The effect is bigger in the ETD of GAθ,v,α than in the ETD of GAθ, since its capability to manipulate more decision variables allows the planner to obtain more solutions that are similar to the myopic criterium and different from the ETD perspective.


## 5.4. Looking for a Drifting Boat in the Sea

The target of the last scenario is a drifting boat (Starget≈5m) close to the coast (Areia Branca, Brazil) within a square area of 15x15km2, whose elevation varies from sea level to 147m. The construction of its initial belief and target motion models has already been described in the Minimum Time Search Definition section.Besides in this scenario, there is a NFZ (e.g., restricted military area) over the coast cells marked in red in Figure10(a), where the green arrows summarize the tendency of the target motion model. Both UAVs enter the search area Ω from the west at t1:2start=0 and stay in the mission up to t1:Uend=400s. Their flying capabilities are different: the first UAV (in red) can do it at h1t∈[500,1500]m and v1t∈[60,90]m/s, while the second (in blue) can do it at h2t∈[500,2000]m and v2t∈[30,60]m/s. The camera azimuth of both α1:2t∈[-90°,90°]. Besides, both UAVs are required to fly at their highest altitude with a FOV=30° to ensure that the sensor likelihood of the first UAV varies from 0.95 to 0.97 and of the second one from 0.89 to 0.91. Additionally, as the estimated smallest complete footprint is 990x821m2, we divide the search area Ω into a grid of 30x30 cells of 500x500m2. Finally, as the drifting boat is a dynamic target, we also have to select its motion lapse. Taking into account the size of the cells, we make Tτ=50s with the purpose of letting the target advance at approximately 10m/s.The results of running both configurations of the planner for this scenario are presented in Figure10. Again, the ETD evolution graph in Figure 10(d) shows that GAθ,v,α produces overall better ETD results than GAθ, at the expense of a bigger variability of its ETD. The initial closer mean values of the ETD results in GAθ,v,α and GAθ are due to the fact that the b(τ0) of this scenario and the initial locations of the UAVs let their cameras initially observe areas with the same probability independently of the UAVs speed and cameras azimuth. The benefit on ETD comes later, when the UAVs with changing velocity can arrive quicker to the high probability areas.For this case, thePd(s1:U0:tend) of the representative solution of GAθ,v,α is also higher than its counterpart in GAθ (74% vs. 71%). That is, the planner with GAθ,v,α lets the UAVs collect quicker more probability. Several benefits of changing vuc,t and αuc,t are clearly observed in the solutions displayed in the second and third columns of Figure 10. For GAθ,v,α configuration, not only does the red UAV arrive sooner to the moving Gaussian, but also the blue UAV orients the camera to cover the coast at its right while flying straightforward, while in GAθ the red UAV follows the track left by the Gaussian and the blue UAV has to change its orientation several times to be able to collect the belief on the coast. Hence, GAθ,v,α obtains safer straighter flights, where the UAV headings changes proposed by GAθ are substituted by camera reorientations. It is also worth noting that in both solutions the red footprint is smaller than the blue, because the first UAV flies at lower altitude than the second.


## 6. Conclusions

This paper presents a new planner for real-world scenarios that optimizes the UAV trajectories and poses of its onboard camera for minimum time search missions with uncertain UAV target locations and movements, and noisy camera-based target detection. This new planner extends the capabilities of the one presented in [10] by (1) permitting the change of UAV heading, speed, and height and of camera azimuth and elevation angles; (2) supporting the use of a DEM and of intelligence information of the search area to define the target behavior (initial belief and motion); (3) modeling realistically (considering the target size, sensor specifications, and DEM) the cameras used by the UAVs to search the target; and (4) handling scenarios with moving targets. Besides, we also present an algorithmic description of the planner and analyze the effects of its different elements in its computational complexity. Finally, we (1) show how to model in detail 3 real-world inspired scenarios (looking for an off-road vehicle in the mountains, survivors in an earthquake, and a drifting boat by the coast) taking into account the influence of different variables (e.g., target size, camera FOV, search area and grid size, and UAV height) and (2) analyze the performance and benefits of the new planner capabilities to tackle them.As future work, we plan to further extend the capabilities of the planner. On the one hand, we would like to let it handle the definition of the probabilistic behavior of multiple targets and the optimization of the expected detection time of all of them. On the other hand, we are considering treating each action type differently, by including weights in the smoothness functions and spending different computational resources in each type (e.g., using a different action stepTa for each type or permitting the multistepped planner to modify each of them in different runs of the GA with a different time computation limit for each type). Besides, we are also exploring the possibility of substituting the GA that supports the optimization of our planner by other metaheuristics, to include new heuristics in the planner for speeding up the optimization, or to generate solutions supported partially by predefined patterns (e.g., spirals, lawnmower) often used in general probabilistic target search problems.



---

*Source: 7673859-2019-05-12.xml*