Coordinated Target Tracking via a Hybrid Optimization Approach

Recent advances in computer science and electronics have greatly expanded the capabilities of unmanned aerial vehicles (UAV) in both defense and civil applications, such as moving ground object tracking. Due to the uncertainties of the application environments and objects’ motion, it is difficult to maintain the tracked object always within the sensor coverage area by using a single UAV. Hence, it is necessary to deploy a group of UAVs to improve the robustness of the tracking. This paper investigates the problem of tracking ground moving objects with a group of UAVs using gimbaled sensors under flight dynamic and collision-free constraints. The optimal cooperative tracking path planning problem is solved using an evolutionary optimization technique based on the framework of chemical reaction optimization (CRO). The efficiency of the proposed method was demonstrated through a series of comparative simulations. The results show that the cooperative tracking paths determined by the newly developed method allows for longer sensor coverage time under flight dynamic restrictions and safety conditions.


Introduction
Recently, unmanned aerial vehicles (UAV) have been widely used in both defense and civilian applications. In many of these applications, the UAV is required to continuously track a moving target, such as in the surveillance and tracking of ground objects. This task becomes challenging when the target is moving in a complex and dynamic environment, where the target may be occluded by ground features, such as buildings and mountains. In the case where the tracking cannot be achieved by using a single UAV, multiple UAVs may be deployed to improve the tracking robustness. The methods to coordinate the movements of multiple UAVs can be generally grouped into two categories, namely centralized and distributed algorithms, respectively [1]. In centralized methods, the mission of each UAV is dispatched from a central station or agent based on global information collected from the all of the aircraft. These methods generally tend to produce better solutions since the global information is available. In contrast, the distributed approaches determine the motion of each UAV by using the local information obtained from the UAV itself and its neighboring UAVs, which enjoys the advantage of computational efficiency and robustness when compared with their centralized counterparts. In a moving target tracking problem, besides generating the tracking path for each UAV, the motion of the target should be estimated since such information cannot be obtained from prior knowledge. As the purpose of this paper is to determine the optimal path of each UAV for cooperative tracking tasks, the motion pattern of the target is assumed to be available when it can be observed by the tracking UAVs.
The aim of cooperative tracking is to keep the target under the coverage of the UAV sensors at all time. To address this problem, Yang et al. adopted the Lyapunov Vector Field Control (LVC) coordinate system, where XI, YI are the Cartesian inertial reference frame. va represents the velocity of the aircraft and denotes the heading of the UAV.
where xi(t) and yi(t) represent the position of the i-th UAV at time instant t. i(t) denotes the roll angle of the UAV, which is bounded to the range withinmax and max. Assuming that each UAV is equipped with an autopilot so that it is able to follow the given roll angle, let Ts denote the sampling rate and we apply zero-order hold to the input signal, then the discrete kinematics of the UAV can be found as: The sensors carried by the UAV can be categorized into to two groups, in terms of the installation conditions, namely gimbaled sensors and body fixed sensors. As shown in Figure 2a, gimbaled installed sensors are capable of pointing to a fixed direction regardless of the attitude of the UAV. The sensor converge area is a cone-shaped region under the aircraft. The pointing orientation of body fixed sensors, on the other hand, are usually affected by the state of the UAV, such as the roll and In this paper, the angle of attack and the angle of side slip are neglected, then the kinematics equations of each aircraft defined in the inertial frame can be represented as: (1) where x i (t) and y i (t) represent the position of the i-th UAV at time instant t. φ i (t) denotes the roll angle of the UAV, which is bounded to the range within -φ max and φ max . Assuming that each UAV is equipped with an autopilot so that it is able to follow the given roll angle, let T s denote the sampling rate and we apply zero-order hold to the input signal, then the discrete kinematics of the UAV can be found as: If the input roll angle is zero, Equation (2) can be rewritten as:

Sensor Visible Regions in Complex Environment
The sensors carried by the UAV can be categorized into to two groups, in terms of the installation conditions, namely gimbaled sensors and body fixed sensors. As shown in Figure 2a, gimbaled installed sensors are capable of pointing to a fixed direction regardless of the attitude of the UAV. The sensor converge area is a cone-shaped region under the aircraft. The pointing orientation of body fixed sensors, on the other hand, are usually affected by the state of the UAV, such as the roll and pitch angles, as shown in Figure 2b. To simplify the problem, in this work UAVs are assumed to be equipped with gimbaled sensors. Hence the sensor convergence on the ground would not be affected by the changes of UAV's attitude.
In many applications, the target to be tracked is moving in a complex environment, which may be occluded by terrain features, such as buildings in an urban environment. Figure 3 illustrates a simulated urban environment, where the 3D rectangular structures represent buildings. Due to the occlusion of buildings, not all of the convergence area of the sensor is visible. The actual visible area depends on the height and the field of view of the sensor. As shown in Figure 3, assuming the sensor is located at (sx, sy, sz), the visible regions on the ground can be determined calculated by using the line of slight (LoS) methods. Rather than using conventional geometrical analysis approaches, such as LoS algorithms, to define the visible areas of the sensors, we employ a fast yet robust spatial analysis method developed in [17], which greatly improves the efficiency of this procedure. pitch angles, as shown in Figure 2b. To simplify the problem, in this work UAVs are assumed to be equipped with gimbaled sensors. Hence the sensor convergence on the ground would not be affected by the changes of UAV's attitude. In many applications, the target to be tracked is moving in a complex environment, which may be occluded by terrain features, such as buildings in an urban environment. Figure 3 illustrates a simulated urban environment, where the 3D rectangular structures represent buildings. Due to the occlusion of buildings, not all of the convergence area of the sensor is visible. The actual visible area depends on the height and the field of view of the sensor. As shown in Figure 3, assuming the sensor is located at (sx, sy, sz), the visible regions on the ground can be determined calculated by using the line of slight (LoS) methods. Rather than using conventional geometrical analysis approaches, such as LoS algorithms, to define the visible areas of the sensors, we employ a fast yet robust spatial analysis method developed in [17], which greatly improves the efficiency of this procedure.    In many applications, the target to be tracked is moving in a complex environment, which may be occluded by terrain features, such as buildings in an urban environment. Figure 3 illustrates a simulated urban environment, where the 3D rectangular structures represent buildings. Due to the occlusion of buildings, not all of the convergence area of the sensor is visible. The actual visible area depends on the height and the field of view of the sensor. As shown in Figure 3, assuming the sensor is located at (sx, sy, sz), the visible regions on the ground can be determined calculated by using the line of slight (LoS) methods. Rather than using conventional geometrical analysis approaches, such as LoS algorithms, to define the visible areas of the sensors, we employ a fast yet robust spatial analysis method developed in [17], which greatly improves the efficiency of this procedure.     Figure 3b depicts the visible areas looked from the point A, the blue regions are invisible areas due to the occlusion of buildings or terrains. Hence, the detectable regions of the sensor can be defined as: where R con is the convergence region of the sensor positioned at p, and R in denotes all of the possible visible areas visible from p. Given the field of view angle of the gimbaled sensor θ, and the height of the sensor h, the convergence region on the ground can be calculated as:

The Motion of the Ground Target
The ground object to be tracked is assumed to move within the 2-dimensional plane, and the associated motion can be described by a state vector x t = [x t , y t , .
y t ]. Let F() represent the state-transition function of the target dynamic, the motion of the target can be described as follows: where W() denotes the state noise and is usually assumed to be subject to a zero-mean Gaussian distribution, i.e., p(W)~N(0, Q), Q is the covariance matrix. The possible motion of the target can be predicted by filtering approaches [13].

The Proposed Method
In this paper, we present a solution to the problem of tracking by multiple coordinated UAVs based on the concept of the Receding Horizon Control (RHC) method, which minimizes the uncertainties arising from the application environment and the motion of the target. The associated optimal control sequence estimation problem is then solved by using a novel bio-inspired optimization approach, based on the framework of chemical reaction optimization.

The Receding Horizon Control in Coordinated Tracking
Since the motion of the target is not available as prior knowledge, it requires predicting the trajectory of the target based on previous observation. Receding horizon control (RHC) approach is an advanced model predictive controller allows for determining the optimal solution at current time instant while taking the future states of the target into account. Let u[k:k + L] denote the control sequence starts from the time k, and L is a constant indicating the window size of the RHC sequence. X[k] represents the state of the dynamic system, including the positions of the UAV and the target, respectively. The main steps of RHC method for solving the target tracking problem are as follows: Step 1: Determine the optimal control sequence based on the state X [k]: where u* denotes the optimal control sequence minimizing the cost function J. Since the dynamic model of the UAV is obtained from prior knowledge, the trajectory of the UAV can be updated in an iterative fashion based on the input control signal and its current state. Because the motion of the target is unknown, its future positions can only be estimated based on previous observation. The cost function J returns the goodness of the tracking based on the relative position between the target and the UAV, which would be discussed in detail in Section 3.2.
Step 2: Apply the optimal input sequence u* to the system, and update the states of the system.

Dynamic Constraint of the UAV
In practice, the roll angle of the UAV should be within a reasonable scale when changing its orientation, due to safety considerations. Hence, the input roll angle should satisfy the following constraint: where φ max is a constant representing the maximal roll angle of the UAV during turning process.

The Target Observation Time
The objective of the coordinate target tracking is to maintain the target located within the detectable region of at least one UAV at any time instant. In addition, it favors the targets located next to the center of the sensor visible area. Let (rx, ry) denote the relative position of the target with respect to the center of the visible area of the sensor: 2σ 2 , the target is within the visible region of at least one sensor P, the target is out of the visible region for all sensors where the function Sen() quantifies whether the target can be observed by the sensors of UAVs based on the relative position between the UAVs and the target, This term reaches its minimum if the target is located within the center of the visible area of all of the sensors, and approaches zero if the target is moving to the boundaries of the sensor visible regions. P is a positive constant used to penalize the case where the target is out of the visible region of all of the sensors.

The Anti-Collision Constraint
In this paper, we assume the UAVs applied to track the target are distributed on the same flight level, thus it is necessary to prevent the UAV from colliding with each other during tracking: where k is a constant, d indicates the desired distance between UAVs, x i and x j represent the positions of i-th and j-th UAVs in the inertia reference frame, respectively. erfc denotes the Gaussian error function which is defined as: This term would reach its minimum when the positions of the UAVs are evenly distributed around the target.

The Regulation Term
In order to ensure the planned path is smooth, the change of the input signal is used as the regulator to penalize a jagged path: This term approaches zero when a straight path is generated, and it would be large in magnitude if the resulting path is highly curved. The total cost can be obtained by summing the aforementioned terms together, leading to the following function: where w 1 , w 2 and w 3 are constants controlling the influence of each term on the total energy.

Chemical Reaction Optimization-Based Coordinated Tracking
The Chemical Reaction Optimization (CRO) algorithm, first reported by Lam and Li [16], is a nature-inspired optimization algorithm based on the principle of chemical reactions. It mimics the process of chemical reactions occuring in a closed container. CRO is a type of population-based algorithm where the basic individuals are molecules. Each molecule is described by some attributes, including the molecular position, i.e., a possible solution of the optimization problem, the potential energy (PE) and the kinetic energy (KE). The molecule structure represents a possible solution of the optimization problem, the potential energy (PE) is associated with the cost of the resulting solution based on the objective function. In general, a small value of the potential energy may indicate a better solution of the optimization problem. Kinetic energy (KE) defines the tolerance of the system for accepting a worse solution than the existing one. In addition to molecules, the container (buffer), where the chemical reaction takes place, is another key component of the CRO algorithm. It is able to supply or absorb the energy resulting during the process of the molecular reactions. Different from many existing population-based optimization algorithms, the total number of molecules may not be the same in each iteration of the optimization process, but the entire amount of energy (the energies of molecules and the energy of the buffer) of the system stays constant throughout the reaction process.

Coding Scheme
As shown in Figure 4, the molecular structure is defined as a sequence of bounded input signals, i.e., the desired roll angle to the autopilot in the time from t k to t k+1 for each UAV. Thus, the dimensions of each molecule would be N × L, where N represents the total number of UAVs employed for coordinated tracking, and L is the number of steps in the predicted ahead current state.
where w1, w2 and w3 are constants controlling the influence of each term on the total energy.

Chemical Reaction Optimization-Based Coordinated Tracking
The Chemical Reaction Optimization (CRO) algorithm, first reported by Lam and Li [16], is a nature-inspired optimization algorithm based on the principle of chemical reactions. It mimics the process of chemical reactions occuring in a closed container. CRO is a type of population-based algorithm where the basic individuals are molecules. Each molecule is described by some attributes, including the molecular position, i.e., a possible solution of the optimization problem, the potential energy (PE) and the kinetic energy (KE). The molecule structure represents a possible solution of the optimization problem, the potential energy (PE) is associated with the cost of the resulting solution based on the objective function. In general, a small value of the potential energy may indicate a better solution of the optimization problem. Kinetic energy (KE) defines the tolerance of the system for accepting a worse solution than the existing one. In addition to molecules, the container (buffer), where the chemical reaction takes place, is another key component of the CRO algorithm. It is able to supply or absorb the energy resulting during the process of the molecular reactions. Different from many existing population-based optimization algorithms, the total number of molecules may not be the same in each iteration of the optimization process, but the entire amount of energy (the energies of molecules and the energy of the buffer) of the system stays constant throughout the reaction process.

Coding Scheme
As shown in Figure 4, the molecular structure is defined as a sequence of bounded input signals, i.e., the desired roll angle to the autopilot in the time from tk to tk+1 for each UAV. Thus, the dimensions of each molecule would be N × L, where N represents the total number of UAVs employed for coordinated tracking, and L is the number of steps in the predicted ahead current state.

On-Wall Ineffective Collision
In this elementary process, the new path is determined by slightly deforming the current path, which can be considered as a local search around the current solution. Let KEω and PEω represent the kinetic and potential energy of a molecule before it hitting the wall of the container, respectively. The potential energy can be calculated based on Equation (12), and the kinetic energy for each molecule is initialized randomly. Due to the collisions, a certain amount of the energy may transfer to the wall of the container (i.e., the buffer), and thus the potential PEω′ and kinetic energy KEω′ of such molecule at the new statue can be computed as:

On-Wall Ineffective Collision
In this elementary process, the new path is determined by slightly deforming the current path, which can be considered as a local search around the current solution. Let KE ω and PE ω represent the kinetic and potential energy of a molecule before it hitting the wall of the container, respectively. The potential energy can be calculated based on Equation (12), and the kinetic energy for each molecule is initialized randomly. Due to the collisions, a certain amount of the energy may transfer to the wall of the container (i.e., the buffer), and thus the potential PE ω and kinetic energy KE ω of such molecule at the new statue can be computed as: Let KElossRate (KElossRate [0, 1]) be a constant defining the percentage of energy that lost before and after the molecule hits the wall, α is a random number evenly distributed from KElossRate to 1. Based on the law of conversion of energy, the energy absorbed by the buffer can be expressed as: If the kinetic energy of the molecule before wall collision is large enough, the after collision potential energy may larger than PE ω , leading to a worse solution when compared with the existing one. It allows the local search algorithm the ability of 'hill-climbing', which increases the probability of avoiding local optimal solutions.

Decomposition
Unlike the process of the 'on-wall ineffective collision', the molecule breaks into parts after it hits the wall of the container (ω→ω 1 + ω 2 ). For simplicity we only consider two parts. Because more solutions are created through the decomposition process, the sum of KE and PE energies of the original molecule may be less than the sum of potential energies of the resulting molecules: In this case, such decomposition is not valid since it disobeys the law of energy conservation. In order to increase the possibility of a correct decomposition process, the buffer may supply a certain amount of the energy to the decomposition process so that the following condition can be satisfied: where δ 1 and δ 2 are two independently random variables uniformly distributed in the range of [0, 1]. The remaining energy is translated as the kinetic energies for two new molecules as: where δ 3 is another independent random variable range from 0 to 1. The decomposition process provides the molecule the ability of exploring more solution regions with respect to the initial position ω.

Inter-Molecular Ineffective Collision
The reaction process is very similar to its single molecule reaction counterpart, where the collision takes place between molecules. The number of the molecules is assumed to be unchanged before and after the collision, and the molecular energies of the associated reaction process should satisfy the condition: The kinetic energies of the new molecules are: The inter-molecule collision allows a median range local search, thus increasing the probability of detecting a better solution in a local area.

Synthesis
Synthesis is the opposite reaction process of decomposition, where two (multiple) molecules merge into a single molecule. As the new molecule is produced from multiple molecules, it is more likely to meet the following requirement for such a process: The remaining energy of the new molecule is: In this reaction process, the new molecule is a combination of multiple old molecules, thus the newly produced molecule may exhibit a large deviation from the original ones in the solution space.
The workflow of the CRO based optimization framework is shown in Figure 5.

Synthesis
Synthesis is the opposite reaction process of decomposition, where two (multiple) molecules merge into a single molecule. As the new molecule is produced from multiple molecules, it is more likely to meet the following requirement for such a process: The remaining energy of the new molecule is: In this reaction process, the new molecule is a combination of multiple old molecules, thus the newly produced molecule may exhibit a large deviation from the original ones in the solution space.
The workflow of the CRO based optimization framework is shown in Figure 5.

Simulation Results and Analysis
This section presents comparative simulation results to demonstrate the efficiency and accuracy of the proposed technique in the coordinated UAV tracking problem. The simulations were conducted using 3D synthetic environmental data, which allows for generating a variety of virtual environments with different resolutions and urban features. The simulations were carried out in MATLAB (R2010b, MathWorks, Natick, MA, USA) on a standard specification PC (Dell Precision T3500, DELL, Round Rock, TX, USA, Inter(R) Xeon(R) CPU operating at 2.67 GHz).
The dimensions of the virtual environment are defined within the range of 1.5 km × 1.5 km, and 30 buildings with different dimensions are randomly generated. Three UAVs are employed to track one moving target, and the motion of the UAVs and target are described using the dynamic model discussed in Section 2. The flight speed vg for each UAV is set to be 30 m/s, the flight level is 200 m,

Simulation Results and Analysis
This section presents comparative simulation results to demonstrate the efficiency and accuracy of the proposed technique in the coordinated UAV tracking problem. The simulations were conducted using 3D synthetic environmental data, which allows for generating a variety of virtual environments with different resolutions and urban features. The simulations were carried out in MATLAB (R2010b, MathWorks, Natick, MA, USA) on a standard specification PC (Dell Precision T3500, DELL, Round Rock, TX, USA, Inter(R) Xeon(R) CPU operating at 2.67 GHz).
The dimensions of the virtual environment are defined within the range of 1.5 km × 1.5 km, and 30 buildings with different dimensions are randomly generated. Three UAVs are employed to track one moving target, and the motion of the UAVs and target are described using the dynamic model discussed in Section 2. The flight speed v g for each UAV is set to be 30 m/s, the flight level is 200 m, the maximum roll angle φ max is set to 45 • and the minimal distance between each UAV is 50 m. The trajectory of the target is randomly generated in the flexible region on the ground, with a constant speed at 15 m/s. The RHC window size L is set as 5 s, and the total simulation time is 100 s, and the sampling time is 1 s. The parameter settings of the CRO algorithm are shown in Table 1.  Figures 6e and 7e, respectively), which approach the total flight time (100 s). This indicates that the proposed method is able to generate an optimal tracking path for each UAV. The costs of the proposed method at each sampling time are depicted in Figures 6f and 7f, respectively. The cost grows rapidly when the target is out of the visibility region for all of the sensors, and for the rest the time the costs are almost constant, implying that the proposed algorithm can determine the optimal path within the sampling time (1 s). As shown in Figures 6d and 7d, the control input for each UAVs, i.e., the roll angle, was changing frequently during the tracking process, since the UAVs are required to maneuver to track the moving target while avoid the occlusions of the buildings. The control inputs determined by the proposed method are generally smooth and changing slowly, which can be tracked by the autopilot of the UAV. the maximum roll angle ϕmax is set to 45° and the minimal distance between each UAV is 50 m. The trajectory of the target is randomly generated in the flexible region on the ground, with a constant speed at 15 m/s. The RHC window size L is set as 5 s, and the total simulation time is 100 s, and the sampling time is 1 s. The parameter settings of the CRO algorithm are shown in Table 1. Figures 6 and 7 illustrate the performance of the proposed method for tracking the moving target in two scenarios. The total visibility time for these two scenarios are 95 s and 94 s (as shown in Figures  6e and 7e, respectively), which approach the total flight time (100 s). This indicates that the proposed method is able to generate an optimal tracking path for each UAV. The costs of the proposed method at each sampling time are depicted in Figures 6f and 7f, respectively. The cost grows rapidly when the target is out of the visibility region for all of the sensors, and for the rest the time the costs are almost constant, implying that the proposed algorithm can determine the optimal path within the sampling time (1 s). As shown in Figures 6d and 7d, the control input for each UAVs, i.e., the roll angle, was changing frequently during the tracking process, since the UAVs are required to maneuver to track the moving target while avoid the occlusions of the buildings. The control inputs determined by the proposed method are generally smooth and changing slowly, which can be tracked by the autopilot of the UAV. Figures 6c and 7c illustrate the relative distances between each UAV during the tracking process for two scenarios, respectively. It can be seen that the generated path by the proposed method always satisfies the safety constraint (i.e., the minimal distance between each UAV). In order to demonstrate the efficiency of the proposed CRO based cooperative tracking path planning algorithm, this section firstly compares the proposed approach with commonly used evolutionary algorithms based optimizer, i.e., the genetic algorithm (GA) [13] and particle swarm optimization (PSO) [14], for determination of the tracking paths. Figures 8 and 9 illustrate the performances of the proposed method, GA-based method and PSObased algorithm for the tracking task, respectively. Figures 8a,b and 9a,b illustrate the convergence of the proposed method and the other evolutionary-based optimizers from randomly selected sampling times for two scenarios. It can been seen that the proposed approach outperforms to GAand PSO-based optimizer in terms of convergence speed and optimality (in terms of the cost value). Figures 8c and 9c depict the visibility time of the ground target of the aforementioned algorithms obtained from 10 repeat simulations. The overall performance of the proposed method is better than GA-and PSO-based planners in terms of the visibility time. This is because the CRO is a variable population-based metaheuristic approach, which allows for generating more solutions near the optimality and thus increasing the convergence speed. Figures 8d and 9d are the average execution times for each decision point (i.e., the sampling point), where the proposed method requires more time than the GA-and PSO-based planners. However, this execution time (around 0.55 s for the proposed method) is acceptable when compared the sampling time (1 s).
Next, the proposed method is compared with a recently developed Lyapunov guidance vector field (LGVF)-based approach for planning the cooperative tracking path [18]. Figures 10 and 11 illustrate the results of using the LGVF-based approach for two scenarios. It can be observed that the proposed method achieves better solution to the cooperative path planning problem in terms of the visibility time. The minimal distances between UAVs are not always satisfactory (i.e., less than the safety distance, i.e., 50 m), as shown in Figures 10c and 11c. The statistics of the proposed method and the LGVF-based approach obtained from 10 repeat simulation tests are shown in Tables 2 and 3, respectively. The mean, maximum and minimum visibility times of the proposed method are better than those of the LGVF-based planner. This is because it is difficult for the LGVF-based optimizer to incorporate nonlinear and non-convex constraints, such as the visibility of sensors and anti-collision constraints in the tracking optimization problem, and thus cannot always reach the optimal solution. Figures 6c and 7c illustrate the relative distances between each UAV during the tracking process for two scenarios, respectively. It can be seen that the generated path by the proposed method always satisfies the safety constraint (i.e., the minimal distance between each UAV). In order to demonstrate the efficiency of the proposed CRO based cooperative tracking path planning algorithm, this section firstly compares the proposed approach with commonly used evolutionary algorithms based optimizer, i.e., the genetic algorithm (GA) [13] and particle swarm optimization (PSO) [14], for determination of the tracking paths. Figures 8 and 9 illustrate the performances of the proposed method, GA-based method and PSO-based algorithm for the tracking task, respectively. Figure 8a,b and Figure 9a,b illustrate the convergence of the proposed method and the other evolutionary-based optimizers from randomly selected sampling times for two scenarios. It can been seen that the proposed approach outperforms to GA-and PSO-based optimizer in terms of convergence speed and optimality (in terms of the cost value). Figures 8c and 9c depict the visibility time of the ground target of the aforementioned algorithms obtained from 10 repeat simulations. The overall performance of the proposed method is better than GA-and PSO-based planners in terms of the visibility time. This is because the CRO is a variable population-based metaheuristic approach, which allows for generating more solutions near the optimality and thus increasing the convergence speed. Figures 8d and 9d are the average execution times for each decision point (i.e., the sampling point), where the proposed method requires more time than the GA-and PSO-based planners. However, this execution time (around 0.55 s for the proposed method) is acceptable when compared the sampling time (1 s).
Next, the proposed method is compared with a recently developed Lyapunov guidance vector field (LGVF)-based approach for planning the cooperative tracking path [18]. Figures 10 and 11 illustrate the results of using the LGVF-based approach for two scenarios. It can be observed that the proposed method achieves better solution to the cooperative path planning problem in terms of the visibility time. The minimal distances between UAVs are not always satisfactory (i.e., less than the safety distance, i.e., 50 m), as shown in Figures 10c and 11c. The statistics of the proposed method and the LGVF-based approach obtained from 10 repeat simulation tests are shown in Tables 2 and 3, respectively. The mean, maximum and minimum visibility times of the proposed method are better than those of the LGVF-based planner. This is because it is difficult for the LGVF-based optimizer to incorporate nonlinear and non-convex constraints, such as the visibility of sensors and anti-collision constraints in the tracking optimization problem, and thus cannot always reach the optimal solution.

Conclusions
Tracking moving targets using UAVs is a challenging task due to the uncertainties arising from the motion of the target and the environmental features. This paper presents a novel method, coupling the concept of model predictive control into the framework of chemical reaction optimization, to solve the coordinated tracking path planning problem in complex urban

Conclusions
Tracking moving targets using UAVs is a challenging task due to the uncertainties arising from the motion of the target and the environmental features. This paper presents a novel method, coupling the concept of model predictive control into the framework of chemical reaction optimization, to solve the coordinated tracking path planning problem in complex urban environments. The sensor visible regions under urban environment, flight dynamics and anti-collision constraints are considered. Particularly, the LoS occlusion caused by dense buildings is discussed, and the FOV for a gimbaled sensor is also analyzed. A series of comparative simulations have shown that the proposed method outperforms other metaheuristic and mathematical methods in terms of tracking ability and safety, while achieving similar computational efficiency. The simulation results imply that the proposed method is capable of improving the performance of target tracking in urban environments, with more visible time, higher security and stronger robustness compared to other methods.