Partial Differential Equation-Based Trajectory Planning for Multiple Unmanned Air Vehicles in Dynamic and Uncertain Environments

Unmanned air vehicles (UAVs) have been used in military operations for a number of years. Recently, UAVs have generated increased interest in civilian domains due to their lower cost and potential applications such as emergency management, law enforcement, precision agriculture, package delivery, and imaging/surveillance. The most fundamental problem for any of these applications is trajectory planning or solving for the shortest path between two points. Trajectory planning is an imperative task required for the navigation and control of UAVs, where the objective is to obtain an optimal or near-optimal flight trajectory between an initial position and the desired goal location under dynamic conditions with environmental constraints [1,2]. There are several important characteristics of an ideal trajectory planner that includes optimality and completeness. Additionally, the trajectory planners need to be computationally efficient since trajectory planning has to occur quickly due to fast vehicle dynamics for many applications using on-board computer resources [3]. The trajectory planning in a large mission area is a typical large-scale optimization problem and several methods have been proposed for solving this problem in the literature. In 1959, Dijkstra introduced his famous method for finding the shortest trajectory using a weighted directed graph [4]. Since then many other analytical solutions to the trajectory planning problem have been proposed using methods such as $A*$⁠, Dynamic Programming, Gradient Algorithms, and Newton Algorithms.

Although these algorithms promise to find the most optimal solution, they are impossible to be implemented in many real world applications due to the amount of computer memory and computational time required [5]. To address these issues, many heuristic methods for trajectory planning have been proposed such as artificial potential field (APF), Genetic Algorithm (GA) [6], fuzzy logic [7,8], Ant Colony Algorithm, Particle Swarm Optimization [9], and Gray Wolf Optimization [10]. Many of these methods implement random sampling of the many trajectory solutions in the search space and through an iterative manner improve the search process. The fundamental issue with these methods is that they cannot guarantee the optimal solution. These methods are advantageous because they require less time and computational effort to solve trajectory planning problems [5]. Real-world applications require very fast algorithms that can handle the complexity of the large-scale problem [11]. Furthermore, the real world complexities induced due to topographical features of the terrain, size of area, and time constraints make the problem even more challenging. The topographical map shown in Fig. 1 represents the typical geographical domain, size, and scale used in this optimization problem.

Partial differential equations (PDEs) are typically used to model system properties related to fluid mechanics and thermodynamics. They have also been proven useful in areas of modeling of mechanical and biological systems [12]. Several studies regarding the application of PDE equations and PDE modeling for control purposes have appeared in the recent literature. For example, Barooha et al. has studied vehicular platooning by applying decentralized bidirectional control [13]. In this study, each vehicle is modeled as a double integrator and then a PDE approximation is derived for the discrete platoon dynamics. Then, the derived PDE model is used to explain the stability of the closed-loop system while adding more vehicles in platoon. Similarly, Frihauf and Krstic [14] introduced a stable deployment method of agents using a PDE-based approach. In their method, the agents' dynamics are modeled by reaction-advection-diffusion class of PDEs. By introducing a feedback law, the agents are allowed to deploy into a family of geometric curves and the stability of the method is then guaranteed by enabling a leader approach. In Ref. [15], Sarlette and Sepulcher proposed a descritized system of agents with local interactions modeled using PDEs. In their study, they found that for coordination of subsystems through local interactions, by assuming the symmetry in translation of subsystems value, distributed control systems can benefit from the discrete approximation of one-dimensional PDEs. Pimenta et al. [16] designed feedback control laws for swarm of robots based on models from fluid dynamics. By applying the in-compressible fluid model, a solution for pattern generation task is generated. In their study, using smoothed-particle hydrodynamics technique and Kernell equations, they aimed to devise a decentralized controller. Similarly, flocking of autonomous agents using a PDE based on Cucker–Smale model is proposed in Ref. [17]. Furthermore, based on the analysis of particle flow, which is a result of PDE Cucker–Smale technique, the stability of the technique is proven. In other study [18], stochastic mapping and coverage strategy using robotic ensembles by PDE-based optimization is studied. In their study, by using an advection–diffusion–reaction PDE model, they task the agents to gather the data and identify the region of interest. Mapping task is then solved as a convex optimization problem using a spatially dependent coefficient in the PDE model. Furthermore for the coverage problem, their study suggests an optimal control problem formulation in which the generated PDE model is expressed as a bilinear control system with the agents' states as the control inputs. In another study, an advection–diffusion equation, which is originally based on Kolmogorov forward PDE equations, is studied for spatio-temporal evaluation of the swarm of agents [19]. Another application of PDE equations in multi-agent systems is introduced by Galstyan et al. [20]. In their paper, they have addressed swarm of chemotactic microscopic robotics based on the diffusion model, introduced by PDE equations and showed the result of their study in one-dimension solution. In Ref. [21], Prorok et al., represented the agents' spatial dynamics over time by PDE based Fokker–Planck diffusion model in a swarm of robots. Another study of swarm, using the same PDE based Fokker–Planck technique, inspired from the honeybee behavior in nature, is addressed in Ref. [22].

Analysis of PDEs, particularly ones representing fluid motion, has been used in the literature for finding optimal vehicle trajectories. For example, stream functions representing motion of fluid particles have been used to produce potential functions for navigation purposes [23,24]. In these studies, after formulating an APF for trajectory planning problem, for a single robot, a simple noncompressible flow equation is implemented to find the global minimum of the APF function. Implementation of a similar concept for path planning of a single UAV in three-dimensional (3D) is provided in Ref. [25].

This paper focuses on the UAV trajectory planning problem in complex environments and discusses the tradeoff between the optimality of path and time of computation. In this paper, a novel 3D trajectory planning method has been proposed which has been inspired by the use of PDEs in fluid mechanics, thermodynamics, circuits, and heat transfer. Particularly, the primary contribution of this paper is formulating and solving the trajectory planning problem in the framework of PDEs that are used to model the fluid flow in porous medium. Partial resistance of fluid flow is among the primary characteristics of the porous media considered in this paper. The fluid particle finds the trajectory from a source to a sink point based on the potential gradient while accounting for the highest permeability zones that provide the least resistance. The nonporous part in the fluid medium is considered to be analogous to threats in UAV environment where UAVs cannot pass through. A risk function is used in this method to model the threats that includes fixed obstacles, detected threats, exclusion zones, fuel limits, required way points, UAV dynamics, and time constraints. The prediction of other aircraft in the airspace is based on their translational movement assumption. The predicted points of collision with the UAV have been weighted appropriately in the risk function. A number of flight simulations were used to generate the results which indicate the robustness of the proposed solution in large and complex domains.

This paper integrates the finite element method for solving PDEs fluid flows in porous medium to motion control policy for multi-UAV system. Hence, for clarifying the structure of this paper, we add the following key to explain the organization and contribution of each of the sections:

• In Sec. 2, the clear formulation for the problem is provided. The constraints and the cost function of the problem are further explained in details.

• Section 3 provides the scheme of PDE solver and the fluid analogy which further assists to understand the proposed approach. General approach in this section is further simplified and a postprocessing section is dedicated to explain how the problem constraints would be enforced on the proposed solution. Since this section provides the solution in a numerical fashion, a theorem is then proved to support the convergence of the solver.

• This paper models the threats and risks encountered by the UAVs as nonporous part of the porous medium. Section 4 is dedicated to present the model of finding the porosity of the environment.

• Section 5 provides an algorithm for further smoothing the trajectory generated by the solver. This section is aimed to manipulate the trajectories to make the resulting paths followable by the UAVs.

• Section 6 provides the flight scenarios used to demonstrate the effectiveness of the proposed algorithm.

This paper considers a problem in which a UAV originates from point $(x0,y0,z0)$ at time t₀ and needs to travel to the goal point $(xf,yf,zf)$⁠. The objective is to reach to the goal position with minimal risk while navigating through a dynamic, uncertain, and risky environment.

The two simple UAV constraints considered in this function are maximum trajectory length and largest achievable trajectory curvature (ψ). The term t is assigned to show that the trajectory R is continuously differentiable function to the degree $t$⁠. Other constraints can be added to Eq. (2) to further formulate the optimization problem. The risk $R(x,y,z)$ is evaluated using a discrete representation of the risk function for $C$⁠.

Based on the fluid analogy, the fluid tends to flow from a higher potential to a lower potential in a field. For instance, consider a water tap. By opening the valve, the water runs from a higher potential (i.e., tap) to a lower potential (i.e., sink), in an optimized manner. Figure 2 represents a schematic of stream lines.

To make the problem more clear, we will introduce the following key for the rest of the paper.

• The UAV is considered as a fluid particle in a flow. The equations of fluid flow govern the motion of the UAV from the initial point to the final point in an optimal fashion [26]. The equations governing the fluid motions and the solution method are provided in this section.

• Each stream line represents a potential path for the UAV. The postprocessing, as explained in Sec. 5 [26], is carried out to ensure the path can be followed by the UAV under its kino-dynamic constraints.

• Potential functions f, originated from the fluid analogy, is designed to represent the higher potentials at the starting point and lower potential at the goal locations.

The variables (i, j, k) in Eq. (5) represent the cell labels of the tessellated area. Equation (5) is solved for the entire domain $C$ while accounting for the boundary conditions.

The results of Eq. (5) are multiple streamlines with varying fluid velocity u at every grid point. Since there is only one starting point and one final point in the domain, the number of solutions are reduced by applying the flux boundary conditional constraints which require that solutions start at t₀ at $(x0,y0,z0)$ and finish at t_f at $(xf,yf,zf)$⁠. The process is explained below.

Equation (6) returns all the routes that travel from the initial point to the goal position and have acceptable risk values.

Variable v represents the cells in trajectory R. The Runge–Kutta method process restricts trajectory length and eliminates undesirable trajectories.

The variable n is the unit vector along x, y, or z direction. These conditions force the generated streamlines to become parallel to the boundaries of domain $C$ when they become sufficiently close to them. This eliminates the chances of trajectories going out of the domain.

In Eq. (13), $f1(x,y,z)$ is a function of absolute values on the boundary of the domain $C$⁠. However, as a result of using Dirichlet conditions on the boundary, the flux conditions are not satisfied. Therefore, the streamlines are allowed to enter and exit from the domain, causing errors in numerical solution over the domain $C$⁠.

The computational domain is a uniform tessellated area in xy plane and is of varying size in the z direction. Another possible method for discretizing the domain is based on finite element methods for obtaining the solution in an unstructured grid. These methods, however, may generate many erroneous results and are computationally expensive [31].

Theorem 3.1. If the topographical function, defining the surface on which the fluid flows, is considered as piecewise flat in domain$C$, the error associated with the numerical solution is negligible.

Proof. The optimal trajectory is computed on a regular x, y grid where elevation is added as the z component. As a requirement of a numerical solution, the error needs to be quantified and proved to be the negligible. In order to prove this theorem, we assume that there is no pressure gradient in y direction. It may be noted that this assumption does not lose generality because the proof for y direction can be obtained in the same way as for the x direction as shown below.

The error is evaluated without losing the problem's generality by solving the equation in a two-dimensional domain.

Equation (26) shows the result of the algorithm and provides the sufficient proof to the theorem. In this case, the method gives us zero error as $g≅0$ and $g′≅0$⁠. From this theorem, it is clear that the result of implementing a nonsmooth topographical map (⁠$g′≫1$⁠) causes $Δ>1$⁠, which is undesirable. For nonflat topographical maps, a solution is presented as follows:

• Tessellation of an area must be done in a way to generate group of cells with flat and smooth surface.

• For solution, the undesired cells which cause a nonsmooth topographical map must be removed.

• A combination of unstructured cells near surface and structured cells above surface should be maintained.

It may be noted that each of the above items affects the computational time requirement of the algorithm, positively.▪

$LC∈R+$ represents a small constant and is considered in this paper to be small value near 0.

where $UC$ represents an upper bound of risk value.

In this paper, the intruder aircrafts (IAs), based on their collaborative nature with the UAV, are treated in two different ways associated with two different risk functions. The first category, identified as “cooperative IA,” is considered as the IAs equipped with Automatic Dependent Surveillance Broadcast (ADS-B) and is continuously reporting their flight data to the other air vehicles. The second category, the noncooperative aircraft, do not share their flight data with other aircraft and the UAV is responsible for predicting the movement of them based on the onboard sensory data. Here, we provide our approach separately for the cooperative and noncooperative IAs.

where ϱ and $φ$ are two positive constants, and $Rc(x,y,z,t)$ is the risk associated with the cell located at (x, y, z) at time t due to presence of cooperative IA.

where the $RH(x,y,z)$ is the risk function associated with the cell located at (x, y, z) due to the high risk area locations.

This section considers the IAs which are identified by the sensors onboard the UAV. These IAs do not report the flight data to the UAV. For every IA, the movement prediction is carried out using a Kalman filter where the position, velocity, and acceleration of the IAs are estimated using the sensory data [32].

and $RNC(x,y,z)$ is the risk function associated with the cell located at (x, y, z) due to the presence of noncooperative IA. It is important to notice that $uC$ is quantified less than $UC$⁠.

The UAV has predefined way-points to visit during its mission before reaching the goal position. The priority of the way-points is defined a priori. The way-points are handled via successively setting each way-point as the goal location, and the analysis is repeated iteratively for each consecutive way point.

The nature of the problem and the boundary conditions determine some properties of the linear set of equations. The aim is to define appropriate Dirichlet boundary conditions on $C$ so that an efficient solver can be found by just considering the symmetric and positive definite values in PDE solutions. The method implemented in this paper is multigrid method with v-cycles. We direct the reader to Ref. [33] for further details on the obtaining an efficient solution of PDEs.

The risk function (⁠$R$⁠) is adjusted based on the new waypoint parameters. In this paper, the UGV is considered as a vehicle capable of following all possible trajectories but for the UAV, the trajectory needs to have certain criteria to be considered as “achievable trajectory.” Due to the dynamics of each vehicle, the trajectories obtained using the proposed technique may be impossible for the vehicle to follow. The algorithm for smoothening the trajectory is utilized after implementing the vehicle dynamics model. Dubins trajectory Model as represented in Ref. [34] is used to generate a trajectory which is achievable by the UAV/UGV. The postprocessing of the trajectory is desirable but it does not guarantee the trajectory will visit every predefined waypoint. Algorithm 1, described ahead, moves the cells responsible for low radius curvatures in a trajectory which a fixed-wing UAV with limitations on turning radius may not be able to follow effectively, generated in domain $C$ as the solution, to achievable positions.

where z(x, y) is the elevation from the surface of ground. However, the altitude is also governed by topographical features (or obstacles on the terrain). The satellite map and processed map (that shows the risk function) of the environment have been shown in Figs. 3(a) and 3(b), respectively. Figure 3(b) shows the risk as the values along the z dimension. Furthermore it is not preferred for the UAV to fly in open areas such as roads and lakes since the chance of detection by radar increases enormously in these areas. Total area of the region in Fig. 3 is $25×108ft2$⁠. Size of the cells in this area is considered as $402ft2$⁠. For generating the risk process image from the satellite map, sketch-up software was used.

In Fig. 4, there are three areas indicated by white color. The two areas on the left-hand side are no-fly-zones over which the UAV is not allowed to fly. The area on the right-hand side is indicated as the area with high probability of UAV colliding with an IA. One way-point is set for the UAV to visit which is indicated by red color in the map. The result of trajectory planning is shown in Fig. 4.

There are several alternative approaches to solve the presented trajectory planning problem. In this paper, we include a comparative study of our proposed method with respect to: GA, Dynamic Programming, and Dijkstra algorithm. The performance is evaluated by comparing the optimality of the solutions obtained from each method. As it is known, dynamic programming represents the exact solution while the GA is a metaheuristic method that can obtain the solution faster. Table 1 represents the average results obtained from 25 runs of the code. The simulations have been performed on a computer with the following specifications:

• Operating System: Windows 8 enterprise 64-bit

• Program Language: matlab R2014a

• Processor: Intel^(R) Core(TM)i7 -2500 CPU @ 3.30 GHz

• (2 CPUs)

• Memory: 8.00 GB RAM

In Table 1, the normalized cost is defined as the average of the cost, calculated by each of the algorithm, running for 25 times, divided by the exact (optimal) cost, found by dynamic programming which is 356927.2.

where $Γ>0$ is the penalty weight and $0<ς≪1$ is a constant. $0≪M$ is the maximum acceleration of the UAV and $a(t)=x¨(t)2+y¨(t)2+z¨(t)2$ represents the acceleration of the UAV at time t. Parameter $b$ is the binary variable which would be equal to 1 if $a(t)≤M$⁠. Otherwise, it is considered as 0. The cost for acceleration exhibits a behavior shown in Fig. 5. Since dealing with $∞$ is a problem in PDEs, ς is set to be very small and then the cost of acceleration passing from $M$ would be extremely high.

The result for a two-dimensional formation scenario is shown in Fig. 6. In this scenario, three UAVs and two IAs are flying in an environment full of uncertainties. The map indicating the unidentified regions has been provided. There is a cost which linearly adds toward the cost function upon entering the uncertain environment. This scenario is based on the real-world scenarios in which the UAVs operate in fully identified regions rather than unidentified regions. The UAVs start from a triangular formation in positions of $(−50,−75), (−50,−25)$⁠, and $(0,−75)$ and constrain to remain in the triangular formation.

The results obtained from the proposed PDE-based method have been verified using the MILP method [32,35] that provides the exact solutions. The results indicate less than 1% error in two-dimensional (2D) analysis.

In order to carry out extensive evaluation and to demonstrate the effectiveness, the proposed technique was applied to several cases in 2D with scenarios of increasing complexity. The results obtained from this study have been provided in Table 2. Here, case IDs represent the scenario descriptions where the first three characters represent the number of obstacles, and the last four characters represent number of UAVs. For instance, “O02U004” means there are two fixed obstacles and four UAVs. Each of the cases was run for 100 instances, and the results shown in the table represent the mean obtained from those runs.

As the results demonstrate, the proposed PDE-based method fares better in computational time as the scenarios get bigger while the optimality of solution does not degrade as much (over 98% optimality accuracy).

This section presents a comparison in results obtained using two cost functions (Eqs. (36) and (37)) for trajectory planning in 3D. The geographical area of the study is Sevilleta National Wildlife Refuge, Albuquerque, NM. The trajectory has been generated by cost functions (36) and (37) have been shown in Fig. 7 by white and black color, respectively. It may be noted that since cost function represented by Eq. (37) represents fuel cost, the trajectory obtained is shorter. While the trajectory obtained by optimizing cost function (36) represents a safer trajectory, it is longer. The result represents the robust performance of the proposed method.

Overall, the results indicate the benefits of utilizing the proposed PDE-based method. Comparisons provided in this paper suggest that the PDE method is capable of finding the acceptable near-optimal solution for large-scale, realistic environments. Previous research indicated the ability of MILP solutions in finding the optimal trajectories for UAVs [32,36], but as it can be seen in Table 2, as the complexity of the problem become larger, the efficiency of the proposed numerical PDE method becomes more pronounced.

As explained in this paper, the proposed PDE method relies on solving PDE equations in an iterative manner. This results in more computational time requirement for cases involving smaller number of UAVs as compared to the other existing methods, such as MILP, for solving similar problems. As it was presented in Table 2, for complex problems, the proposed PDE method becomes more computationally efficient as compared to the MILP.

The paper presents the formulation of UAV trajectory planning within the framework of PDEs and obtains the near-optimal solution via solving the PDE using an efficient numerical method. The presented method utilizes the theory of fluid flow in porous media to solve the problem of UAV trajectory planning in region filled with moving and stationary obstacles. The proposed method has proven desirable as it is numerically solved and shown to provide acceptable solutions that require less computational time. The method is compared to several different methods and is shown to provide computational advantage without much compromise in the optimality of the solution. Specifically, the proposed method is applied to actual flight scenarios in large-scale problems and the results show that the method has better performance than the GA and Dijkstra in terms of optimality and computational time. Furthermore, in comparison with the Dynamic Programming and MILP, the method provides solution in much less computational time at a small sacrifice in optimality.

The problem studied in this paper involves a dynamic environment with moving obstacles and intruder aircraft. As a part of future work, a number of steps can be performed to make the computations faster. These include parallelizing the solver [37] and multikernel learning [38]. Future works also include generating the decentralized law based on the dynamical behavior of fluid particles for the purpose of computational efficiency. This would allow such methods to be used for more complex scenarios. Another future direction of research would be to incorporate a human-in-the-loop study for the scenarios such as area coverage problems.