A Distributed Algorithm for Real-Time Multi-Drone Collision-Free Trajectory Replanning

In this paper, we present a distributed algorithm to generate collision-free trajectories for a group of quadrotors flying through a common workspace. In the setup adopted, each vehicle replans its trajectory, in a receding horizon manner, by solving a small-scale optimization problem that only involves its own individual variables. We adopt the Voronoi partitioning of space to derive local constraints that guarantee collision avoidance with all neighbors for a certain time horizon. The obtained set of collision avoidance constraints explicitly takes into account the vehicle’s orientation to avoid infeasiblity issues caused by ignoring the quadrotor’s rotational motion. Moreover, the resulting constraints can be expressed as Bézier curves, and thus can be evaluated efficiently, without discretization, to ensure that collision avoidance requirements are satisfied at any time instant, even for an extended planning horizon. The proposed approach is validated through extensive simulations with up to 100 drones. The results show that the proposed method has a higher success rate at finding collision-free trajectories for large groups of drones compared to other Voronoi diagram-based methods.


Introduction
Trajectory generation is a key element for the execution of complex autonomous vehicle missions. It can be defined as the computational problem of finding a valid trajectory that guides a vehicle from an initial state to a given final state in an environment with static and/or moving obstacles. In most applications, the main concern, rather than just finding a feasible trajectory between the initial and final states, is to obtain the optimal trajectory with respect to a certain objective function. In such a setting, trajectory generation is formulated as an optimization problem with a cost function that quantifies the accomplishment of mission goals and objectives, and different types of constraints to ensure safety and feasibility of resulting trajectories.
With rapid advances in communication and computational technology, autonomous vehicles continue to take part in more complex missions, and even engage in teams of collaborating vehicles to take on increasingly demanding tasks. This necessitates incorporating intervehicle collision avoidance constraints in the optimization problem to guarantee that generated trajectories for a group of vehicles sharing a common workspace are collisionfree. Therefore, for a large group of vehicles, the optimization problem would involve a large number of constraints and decision variables, and the computational cost of solving it centrally can be prohibitively high. To reduce the computational complexity, a multitude of distributed schemes, reviewed below, have been proposed for decomposing the optimization problem into smaller sub-problems that can be solved locally by each vehicle. The major challenge is to ensure that local decisions do also satisfy the coupling collision avoidance constraints. This is mainly addressed by exchanging information among the vehicles on their current states, future input sequences, etc. Depending on the communication strategy, the sub-problems might be solved sequentially or concurrently, with possibly several iterations of optimization and communication to achieve the required performance.
In [1], the collision avoidance constraint, usually expressed in terms of the two-norm of a relative position vector, is approximated by a set of linear constraints. The sub-problem for each vehicle is then formulated as a mixed-integer linear programming (MILP) that includes the vehicle's individual variables as well as variables of a subset of neighbors. This enables cooperation among vehicles by allowing a vehicle to make feasible perturbations to neighboring vehicles' decisions. The sub-problems are solved sequentially by each vehicle, and the algorithm iterates over the group of vehicles until convergence, during each cycle of a model predictive control (MPC) scheme.
Sequential convex programming (SCP)-based methods have also been used for solving distributed multiple vehicle trajectory generation problems [2,3]. The work in [4] addresses the infeasiblity of intermediate problems in decoupled-SCP methods, arising from convex approximation of collision avoidance constraints, i.e., linearizing them, and proposes incremental SCP (iSCP), which tightens collision constraints incrementally. Compared to sequential approaches in [5][6][7], that cast the trajectory of anterior vehicles as dynamic obstacles for a posterior vehicle, the methods proposed in [1,4] result in less constrained intermediate problem and faster convergence rate, yet, similar to most MPC-SCP-based methods, they would require the vehicles to exchange a full representation of their decisions to neighboring vehicles over a communication network.
The synchronous approach in [8] extends the distributed MPC (dMPC) scheme in [9] for formation control, based on alternating direction method of multipliers (ADMM), to problems with intervehicle collision avoidance constraints. These constraints are decoupled using separating hyperplanes, which enforces each vehicle to stay within one half-space of a time-varying plane over a certain time horizon. The resulting sub-problems are solved simultaneously by vehicles, while the normal vector and offset shared between a vehicle and a neighbor, for characterizing their separating hyperplane, are updated at each cycle of the dMPC, using the interchanged information about generated trajectories at the previous cycle.
In the decentralized trajectory planner proposed in [10], vehicles replan their trajectories asynchronously, independent of the planning status of other vehicles. At each iteration, a vehicle considers trajectories assigned to neighboring vehicles as constraints, and solves an optimization problem including as decision variables the normal and offset of planes that separate the outer polyhedral representation of its trajectory and those of its neighbors. A check-recheck scheme is then performed to ensure that the generated trajectory does not collide with trajectories other vehicles have committed to during the optimization time. Therefore, to guarantee deconfliction between vehicles, the planner requires a vehicle to broadcast its computed trajectory to its neighboring vehicles at the end of each replanning iteration.
The on-demand approach to local collision avoidance, proposed in [11], imposes constraints only at specific time instances when collisions between a vehicle and its neighbors are predicted. Predicting collisions along a time horizon, however, relies on an accurate knowledge of the neighbors' future actions which must be communicated at every sampling time. The dMPC scheme in [12] for distributed trajectory generation is based on this predict-avoid paradigm and an event-triggered replanning strategy, and has been shown to result in less conservative trajectories, but at the cost of voiding the collision avoidance guarantees for all time instances over the horizon. To capture the downwash effect of quadrotor's propellers, the collision avoidance constraint in [12] is modified with a diagonal scaling matrix, which approximates the quadrotor body with a translating ellipsoid elongated along the vertical axis, yet it ignores the quadrotor's rotational motion.
Reciprocal velocity obstacle (RVO) and its variants have been widely used in distributed collision avoidance [13][14][15][16][17]. At each time step, RVO [13] builds the set of all relative velocities, leading to a collision between a vehicle and its neighbors, and chooses a new constant velocity outside this set, and closest to the desired value, to avoid collisions. Therefore, RVO requires the position and velocity information to be communicated, or sensed, between nearby neighbors. Other variants, such as acceleration velocity obstacle (AVO), which addresses the instantaneous change of velocity in RVO by taking into account acceleration constraints, need further information such as acceleration to be interchanged. Reciprocally-rotating velocity obstacle (RRVO) [18] uses rotation information to mitigate deadlocks caused by symmetries of representing vehicles with translating discs in RVO. It relies on the assumption that neighbors may rotate equally (or equally opposite), bounded by a maximum value, to compute an approximation of swept areas for rotating polygonshaped vehicles, and uses them for constructing velocity obstacles. A new velocity and rotation is then selected at each time step to avoid collisions.
Another approach to distributed collision avoidance is to construct the Voronoi diagram of the group of vehicles and generate the trajectory for each vehicle so that it is entirely within the vehicle's Voronoi cell [19][20][21][22]. Since Voronoi cells do not overlap, it can be guaranteed that the generated trajectories are collision-free. To consider the physical size of a vehicle, the modified Voronoi cell used in [23,24] retracts boundary hyperplanes of the general Voronoi cell by a safety radius for disc-shaped vehicles. At each sampling time, upon receiving the relative position information, trajectories are replanned to conform to the updated Voronoi diagram. The resulting sub-problems can be solved simultaneously, in a receding horizon manner, until the vehicles reach their final positions. The Voronoi-based approaches only require the vehicles to know relative positions to neighboring vehicles, and therefore are well suited to applications where vehicles only have relative position sensing and no communication network [25].
In this paper we develop a distributed trajectory generation framework, with low computation and communication demands, for multiple quadrotors flying in (relatively) close proximity to each other. We specifically address the shortcomings of approximating the drone body with a disc (or sphere) for generating feasible collision-free trajectories for large groups of drones. A sphere model, used in most existing distributed collision avoidance schemes, may be overly conservative in confined spaces since it invalidates trajectories whose feasibility depends on the consideration of the flight attitude. Instead, we model the drone body with an ellipsoid, and employ the Voronoi partitioning of space to derive local collision avoidance constraints that take into account the drone's real size and orientation. The same approach can be integrated into other distributed schemes that utilize separating hyperplanes for decoupling collision avoidance constraints. Yet the main reason for adopting Voronoi diagram is that using time-invariant boundary hyperplanes determined prior to solving a sub-problem, despite being more conservative, can significantly reduce communication and computational load, allowing for higher replanning rates. Incorporating the resulting set of constraints into sub-problems, solved by each vehicle, allows finding collision-free trajectories for guiding a group of drones through confined spaces by proper adjustment of attitude angles. In addition, the obtained set of constraints can be expressed as Bézier curves, and hence can be efficiently evaluated to guarantee that intervehicle collision avoidance requirements are met at any instant of time even over a long planning horizon.
In the proposed synchronous distributed scheme, each vehicle uses the position information of its neighbors, updated at each sampling time, and solves a sub-problem to generate its trajectory inside (a subset of) its Voronoi cell towards the closest point (in the cell) to its goal position. We present an efficient method to compute this point, which is needed to appropriately define the terminal constraint and cost in the sub-problem. A sequence of sub-problems are then solved in a receding-horizon manner until the vehicles reach their goal positions. The simulation results show that the proposed method has a higher success rate at finding collision-free trajectories for larger groups of quadrotors compared to other Voronoi diagram-based methods. In addition, it can effectively reduce the total flight time required to perform point-to-point maneuvers. Furthermore, the computation time of generating those trajectories satisfies timing constraints imposed by real-time applications.
The rest of this paper is organized as follows: In Section 2 we formulate the optimization sub-problem solved by each vehicle. In Section 2.1 we study the differentially flat system describing the drone equations of motion, and parameterize trajectories with Bézier curves. We derive the set of local collision avoidance constraints in Section 2.2, and present an efficient algorithm for finding the closest point in a Voronoi polytope to a goal position in Section 2.3. In Section 3 we obtain the continuity conditions between two adjacent Bézier curve segments, and present a method for evaluating inequalities in Bézier form without discretization. Finally, we provide simulation results in Section 4.

Problem Formulation
The multiple vehicle trajectory generation problem addressed in this paper can be defined as finding optimal trajectories that act as references to guide a group of vehicles from their initial positions to some desired final positions. The generated trajectories should jointly minimize a cost function, corresponding to the accomplishment of mission goals and objectives, and satisfy a set of local and coupling constraints, so that they are dynamically-feasible and collision-free. For N v vehicles, this problem can be formulated as the following optimal control problem.
The cost to be minimized is the sum of the vehicles' individual costs, J, given by the functional, where x i (t) ∈ R n x and u i (t) ∈ R n u are the state and the input vectors of the vehicle's model described by an ODE, and x i,0 and x i, f are the initial and final values of the state of the i-th vehicle, respectively. X i and U i denote the set of admissible states and inputs for the i-th vehicle derived from limits imposed by vehicle dynamics and the surrounding environment. In order to reduce the computational complexity of solving (1) for larger N v with increased numbers of constraints and variables, one can divide the problem into a set of small-scale sub-problems. Here, the sub-problems are formulated such that each involves only a vehicle's individual costs and constraints, and hence can be solved independently by the vehicle. The sub-problems must include constraints to ensure that the trajectory generated locally by a vehicle does satisfy the coupling collision avoidance constraints.
The key idea to ensure intervehicle collision avoidance is to decompose the space into non-overlapping regions, provided by a Voronoi diagram, and generate the trajectory for each vehicle such that it is entirely within its partition. The Voronoi diagram is updated at each sampling time according to the relative positions of vehicles, and a sequence of sub-problems is solved in a receding horizon manner until the vehicles reach their final positions. For the i-th vehicle, the problem that has to be solved at the time instant t k can be formulated as min where x i,k (t) and u i,k (t) are the state and the input profiles of the vehicle over the time interval [t k , t k + t h ], with t h being the planning horizon, andx i,k denotes its state at the time instant t k . The cost function in the above sub-problem is modified as where the second term is added to penalize the distance, at t k + t h , to the point in the Voronoi partition that is closest to the goal position.
In the optimization problem (3), C i,k may denote the Voronoi partition assigned to the i-th vehicle. The Voronoi diagram is updated for each sub-problem according to the vehicles' configuration at each time instant t k , i.e., Since Voronoi partitions are disjoint and the assigned trajectory to each vehicle for the time horizon t h is contained within its partition, it can be guaranteed that there is no collision between the trajectories over the time The distributed trajectory generation framework is summarized in Algorithm 1. In Section 2.2, we study the Voronoi diagram for a group of vehicles and modify C i,k to explicitly take into account the orientation while generating collision-free trajectories for multiple drones. Solve the optimization sub-problem 11:

Quadrotor Model
In this paper, the simplified quadrotor equations of motion are described by where p ∈ R 3 is the position and m is the mass of the quadrotor. In addition, g = 9.8 m s 2 is the gravitational acceleration, and e 3 = [0 0 1] T . The first term on the right-hand side of (5) is gravity in the z I direction, and the second term, f ∈ R 3 , is the thrust force aligned with the body's z-axis.
where T ∈ R is the net thrust, I z B = R B z B = Re 3 is the body frame z-axis expressed in {I}, and R ≡ I B R ∈ SO(3) is the rotation matrix from the body frame {B}, centered at the quadrotor's center of gravity, to the fixed inertial frame {I}. For simplicity, we drop the superscript I and consider z B = I z B . Figure 1 is a graphical representation of the quadrotor and the associated reference frames.

Trajectory Parametrization
The quadrotor dynamics (5) with the four inputs is differentially flat [26], and therefore the state and the input of the system can be expressed as functions of the flat outputs and a finite number of its derivatives. The position vector together with the yaw angle can be selected as flat outputs of the system. Here, p ∈ R 3 is parameterized as a Bézier curve, given by wherep l ∈ R 3 are the control points, B l,n are Bernstein basis polynomials of degree n, and τ ∈ [0, 1] is defined as The linear velocity, v =ṗ, and linear acceleration, a =p, can be expressed as parametric Bézier curves through the first and second derivative of p with respect to time, yielding where the control pointsv l andā l are obtained as The thrust T and rotation matrix R can also be expressed as functions of the flat output and its derivatives. The net thrust T can be written as Assuming that the rotation matrix then the columns of the rotation matrix are extracted from where the unit vector r is defined as The above equations declare that the vehicle's orientation can be fully determined from the second derivative of the trajectory and the yaw angle. As mentioned before, the yaw angle ψ is a component of the flat output, and therefore it can be controlled independently without affecting the trajectory generation. Using the differential flatness property of the system, trajectories consistent with dynamics can be planned in the space of flat outputs, where (5) is trivially satisfied and the original input and state constraints are transformed into constraints on the flat output and its derivatives.

Collision Avoidance
In this section, we present a Voronoi diagram-based approach to decoupling intervehicle collision avoidance constraints. Although the presence of obstacles, interpreted as non-decision-making agents, is not explicitly considered here, incorporating vehicleobstacle collision avoidance constraints into the problem simply amounts to taking into account the obstacles' position when updating the Voronoi partition (step 6 of Algorithm 1).
The widely used approach in the literature to avoiding collisions with obstacles in the environment is to model the drone body as a sphere with radius r D , and then simply building the collision-free space, C f ree , by inflating the obstacles with a factor r D . As a result, collision-free trajectories can be obtained by enforcing the vehicle, which is now treated as a point in the space, to be inside C f ree [27]. Considering now the collision avoidance between the i-th and j-th drones, the corresponding constraint can be derived similarly by where . denotes the Euclidean distance. Ignoring the real shape and orientation of the drone, and approximating its body with a sphere, invalidates trajectories that are feasible upon considering the flight attitude. For this reason, the above approach might be too conservative for trajectory generation in confined spaces and can even fail to find feasible collision-free trajectories when multiple drones are involved. Approximating the drone body with an ellipsoid, whose principal axes are aligned with the body frame axes, allows considering the drone orientation while inspecting for collisions against other vehicles. For the i-th drone, the ellipsoid, E i , centered at the drone position p i , is given by where Λ is a 3 × 3 diagonal matrix of the form with r D and h D being the lengths of the principle semi-major and semi-minor axes, respectively. (See Figure 2). As proposed in [28], collision avoidance constraints for two ellipsoid-shaped drones can be derived using separating hyperplanes. Denoting by a ∈ R 3 and b ∈ R the normal vector and offset of a hyperplane, respectively, the separating hyperplane for E i and E j , defined as H ≡ {p|a T p − b = 0}, must satisfy The set of inequalities (18) holds if, and only if, Satisfying the set of constraints (20) will guarantee that there is no collision between the two ellipsoids E i and E j associated with the i-th and j-th drones, respectively.
For multiple vehicle trajectory generation, collision avoidance constraints, either in the form of the inequality constraint (15), for spheres, or the set of constraints (20), for ellipsoids, must be incorporated in the optimization problem for each pair of vehicles. As the number of vehicles involved in a mission grows, the resulting increase in the number of constraints would inevitably exacerbate the computational issues of finding collision-free trajectories in a centralized manner.

Distributed Collision Avoidance
Here, we propose a distributed approach to collision avoidance which takes into account the shape and orientation of a drone. The presented approach uses Voronoi partitioning of space and generates the trajectory of each vehicle such that it is entirely within (a subset of) the vehicle's Voronoi cell for a time interval t h . Since Voronoi cells are pairwise disjoint, and each vehicle only moves inside its Voronoi cell, intervehicle collision avoidance is guaranteed for all future time before t h .
Each Voronoi cell in an n-dimensional space is a convex polytope bounded by a number of (n − 1)-dimensional convex polytopes. For a group of vehicles in three-dimensional space, the general Voronoi cell of the i-th vehicle is defined as where p ij = p j − p i , and p i and p j are the position of the i-th and j-th vehicles at the current time instant. Note that V i is the intersection of half-spaces corresponding to hyperplanes with a = p ij and b = 1 2 p T ij (p i + p j ). An arbitrary point in V i is closer to the i-th vehicle than any other vehicle [22], i.e., The boundary of the Voronoi cell, ∂V i , is the union of multiple faces, each of which include points in the space that are equidistant to the i-th vehicle and a neighboring vehicle.
In order to account for the size of a vehicle, the buffered Voronoi cell (BVC) proposed in [24] retracts the boundary of the general Voronoi cell by a safety radius, so that if the vehicle's center is inside the BVC, its body, approximated by a sphere of radius r D , will be entirely within its Voronoi cell. The BVC of the i-th vehicle, denoted byV i , is defined as It can be easily shown that BVC is a subset of the general Voronoi cell, i.e.,V i ⊂ V i . In addition, for any two points p ∈V i and q ∈V j , p − q ≥ 2r D holds. Therefore, the vehicles are guaranteed to avoid collisions due to the buffer region of r D along ∂V i . Figure 3 shows the Voronoi diagram for 10 drones in a collision-free configuration and the BVC for two adjacent drones. The BVC is defined based on a symmetrical approximation of the vehicle's body with a translating disc. In order to reduce the conservatism and avoid infeasibility issues due to ignoring the real shape and orientation of the vehicle, we approximate the drone with an ellipsoid (16), and bearing in mind that we propose C i in problem (3) to be defined as If the trajectory of the i-th drone p i (t) satisfies the above set of local collision avoidance constraints for all t ∈ [t k , t k + t h ], then the ellipsoid representing the drone body is within the Voronoi cell for the entire time horizon; that is, the ellipsoid centered at p i and aligned with the columns of R i does not intersect the Voronoi boundary, stated mathematically Noting that it can be induced thatV where proj XYZ C i is the projection of C i onto the three-dimensional subspace spanned by e 1 , e 2 , and e 3 . Therefore, incorporating (25) into the optimization problem (3) will ensure that the generated trajectories are collision-free while alleviating infeasibility problems by taking orientations into account. Also, since z B is fully obtained fromp (13), the above set of local collision avoidance constraints can be expressed as constraints imposed on Bézier curves. Later, we present an efficient method for evaluating inequalities in Bézier form.

Finding the Closest Point to the Goal Position
As explained above, at each time instant the Voronoi cell V i is updated according to the relative position of the i-th vehicle to other vehicles. The optimization problem (3) is then solved to generate a trajectory, for a time horizon t h , that guides the vehicle towards the point in the cell closest to the goal position. This process is repeated until the vehicle reaches its final position. At each sampling time, the closest point must be found prior to solving the trajectory generation problem. Therefore, having an efficient scheme for finding the closest point is critically important to avoid long computational delays between updating the Voronoi cell and replanning the trajectory.
The point in a convex polytope that is closest to a query point q is either q itself or a point on the boundary of the polytope. A naive way to find the closest point in a convex polytope in a three-dimensional space, represented by P = (F, E, V), where F is the set of faces, E is the set of edges, and V is the set of vertices, is to check the distance between q ∈ R 3 to all faces, edges, and vertices for finding the minimum. However, for complex polytopes, the computation time is not negligible.
The geometric approach proposed in [24] for a polygon in a two-dimensional space calculates the barycentric coordinates and an angle from the query point to the two vertices of each edge to find the closest point. Since this approach iterates over all edges, its computational complexity can significantly increase as the number of Voronoi neighbors of a vehicle increases. Here, we make use of the Gilbert-Johnson-Keerthi (GJK) distance algorithm and devise an approach that can efficiently determine whether the query point is inside the polytope, i.e., q ∈ P, in which case the closest point is q itself. Otherwise, the presented algorithm returns the closest feature of P to q, and the closest point can be obtained by projecting q onto it. The proposed approach is not limited to distance queries between a point and a polytope, and can also be used when the final constraint in (3) is relaxed to a small box or sphere around the goal position, and used in conjunction with a terminal cost term.
The GJK distance algorithm or simply GJK algorithm is an iterative algorithm that relies on a support mapping function to incrementally build simplices that are closer to the query point [29]. The algorithm has been extensively used for collision detection between generic convex shapes [30,31]. The original algorithm, however, can be used to compute the minimum distance, and also the respective pair of (closest) points, between two convex shapes [32].
In order to obtain the minimum distance between two general convex sets A and B, GJK approximates the closest point in the Minkowski difference of the two sets, C = A − B to the origin, denoted by ν(C), with an iterative search. At each iteration, a simplex in C is constructed such that it is closer to the origin O than the simplex in the previous iteration. In R 3 , a simplex can be a point, a line, a triangle, or a tetrahedron with 1, 2, 3, and 4 affinely independent vertices.
GJK relies on the so-called support mappings to construct a new simplex. A support mapping function s C (d) of the convex set C maps the vector d to a point in the set, called the support point, according to At each iteration, a support point w k = s C (−ν k ) is added as a vertex to the current simplex, indicated with V k , where ν k is the closest point of V k to the origin, i.e., ν k = ν(conv(V k )). V k+1 is then updated such that it only contains the smallest set of vertices that supports ν k+1 = ν(conv(V k ∪ w k )), and earlier vertices that do not back ν k+1 are discarded [30].
It is proved in [30] that in each iteration the new ν is closer to the origin than the previous one, and thus the sequence of {ν k } converges to the closest point of C to the origin. In addition, it is shown that which is used to construct the terminating condition of the GJK algorithm for general convex shapes. The GJK algorithm, as explained above, depends heavily on the computation of ν k to test the termination condition and to determine the search direction for finding the support point. In each iteration of the algorithm, ν k must be computed with the Johnson Distance Subalgorithm [29] or more robust alternatives such as the signed volume method in [33]. Here, we exploit unique features of polytopes and propose a faster way to evolve simplices in the GJK algorithm without computing ν k in each iteration.
For polytopes, GJK arrives at the actual ν(C) in a finite number of iterations [30]. The pseudocode in Algorithm 2 describes the GJK distance algorithm for a polygon P. In order to find the support point w k , we employ a search direction d k ↑↓ ν k , which is updated in each iteration of the algorithm with S1D, S2D, or S3D subroutines. To update d and V, these subroutines, summarized in Algorithms 3-5, inspect the Voronoi regions of the simplex for the one that contains the origin. Figure 4 shows the Voronoi regions of an m-simplex (m = 1, 2, 3) where the origin possibly lies. Once the Voronoi region containing the origin is found, d is determined as a vector from the associated vertex, edge, or face (of the simplex) pointing towards the origin.
The stop criterion for the conditional loop in Algorithm 1 is also constructed using the search direction, offered as where v 1 ∈ V k . Considering that d T k (ν k − v 1 ) = 0, we can conclude that the above criterion, for deciding whether V k represents the closest feature of P to the origin, is equivalent to the stop criterion in [30] for determining whether ν k s the closest point, that is If the inequality (31) holds, then the closest feature of P to the origin can be determined from V k . Figure 5 shows all possible outputs of Algorithm 2, which can be a vertex, an edge, or a face of P or a tetrahedron inside it, and ν(P) for each of the cases. A support point of a convex polytope can also be computed efficiently. For a polytope P, the support point is a vertex of P, i.e., s P (d) ∈ vert(P), and we can take w = s P (d) Therefore, for polytopes, the support point can be uniquely determined by simply scanning through the list of vertices for the vertex that is the most extreme in the search direction d. Therefore, the computation time is linear in the number of vertices of P. For complex polytopes, the vertices adjacency information and the coherence between consecutive calls to support mapping functions can be exploited to find the support point with almost constant time complexity [30].
end if 16: end if 17: end if 18: end function Figure 5. Examples of the closest feature of a polyhedron to a query point are shown above. Once the closest feature is obtained from Algorithm 1, the closest point, i.e., ν(P), can be determined, as shown above.
else 10: end if 16: end if 17: end if 18: end function

Continuity Conditions
As explained before, at each sampling time, a trajectory, expressed as a parametric Bézier curve, is generated for the time horizon [t k , t k + t h ], and the trajectory for the entire flight time [0, T] is formed by joining segments of these Bézier curves end-to-end. The smoothness of the resulting composite trajectory must be guaranteed by enforcing continuity at the joining points of two consecutive segments up to a certain derivative. In the following, in order to derive conditions that address parameter continuity between consecutive curves, we assume that the time horizon is equal to ∆t k = t k+1 − t k , which is not necessarily the same for all sub-problems. In practice, however, the time horizon is greater than ∆t k , in which case a Bézier curve describing the segment over the time interval [t k , t k+1 ] can be obtained by subdividing p k (.) at t k+1 with the de Casteljau's algorithm. For simplicity we drop the subscript i ∈ [N v ].
The Bézier curve describing the trajectory over the time interval [t k , t k+1 ] is defined as Assuming that the global parameter t runs over the interval [t k , t k+1 ], the local parameter τ k is related to t by The parametric continuity condition, C r , requires the r-th derivative and all lower derivatives of two consecutive segments to be equal at the joining point. In other words, Zero-order parametric continuity, C 0 , requires the endpoints of two consecutive curves, p k (.) and p k+1 (.), to intersect at one endpoint, that is, Since a Bézier curve is coincident with its control points at the two ends, i.e., the position continuity condition (37) translates intō The first-order parametric continuity condition, C 1 , for the k-th and k + 1-th Bézier curves, can be obtained as ∆t k+1 n k (p n k ,k −p n k −1,k ) = ∆t k n k+1 (p 1,k+1 −p 0,k+1 ) Finally, the k-th and k + 1-th Bézier curves are C 2 -continuous if Higher-order parametric continuity conditions can be obtained likewise.

Evaluating Inequalities in Bézier Form
Parameterizing the trajectory with a Bézier curve converts the original infinite dimensional problem (3) into a semi-infinite optimization problem with a finite number of decision variables and an infinite number of constraints. The commonly used approach to obtaining a standard optimization problem is time gridding, which inspects satisfaction of constraints on a finite number of points only. Although this method is straightforward, it cannot guarantee that constraints are satisfied over the entire time interval. Using fine discretization can remedy this issue, but, it will increase the number of constraints as well as the computation time. Since all constraints involved in the trajectory generation problem addressed above can be expressed as Bézier curves, we can employ the method proposed in [33] to recast the semi-infinite optimization problem into one that is computationally tractable. As explained below this method exploits unique features of Bézier curves to efficiently evaluate constraints while avoiding problems associated with time gridding.
If h(τ) can be expressed as a Bézier curve, then any inequality constraint of the form h(τ) ≤ 0, τ ∈ [0, 1] can be replaced by a finite set of constraints on its control points. More specifically, if h(τ) is defined as then from the convex hull property of Bézier curves we know that where CH(H) = {α 0h0 + · · · + α n hh n h |α 0 + . . . , α n h = 1, α l ≥ 0} is the convex hull defined by the set of control points [34]. Thus, the inequality constraint h(τ) ≤ 0 holds if This finite set of inequality constraints can ensure that the original inequality constraint is satisfied over the entire interval [0, 1]. However, Ineqs. (44) might be conservative due to the existing gap between the control pointsh l and the actual curve h(τ). This problem can be alleviated by refining the control polygon and finding closer control points to the curve using recursive subdivision of h(τ) with the de Casteljau's algorithm. The sequence of control polygons generated with repeated subdivision converges to the underlying Bézier curve [35]. Figure 6 shows a threefold subdivision of a cubic Bézier curve. Furthermore, the de Casteljau's algortithm allows refining the control polygon locally. Using recursive subdivsion of h(τ) to reduce the conservatism in the finite set of constraints results in an increase in the number of constraints; hence, a trade-off has to be made between the computational effort and the conservatism. Nevertheless, the optimization variables remain the same [36].

Simulation Results
In this section, the efficacy of the proposed method for generating feasible and collisionfree trajectories in (vehicle-) dense environments are assessed through different simulation examples. We compare the resulting trajectories to those generated with the well-studied BVC approach. We specifically test the capability of the two methods to generate trajectories that ensure all drones involved in a simulation example reach their final positions, and compare the flight time, obtained with each of them, to complete point-to-point transition missions. We also present the recorded computation time for executing the proposed algorithm in this paper to emphasize its suitability for real-time applications.
In the simulations presented below, we assume all drones have the same size, and their BVC (23) (4) i,k (τ) 2 dτ. The time horizon and the replanning step are also considered to be the same for both methods. The obtained solution at the previous replanning step is used to set the initial guess for the current sub-problem. We use FORCES Pro [37] to generate solvers for the resulting sub-problems. The sub-problems, involving the set of control pointsp i,k as decision variables, can be reformulated to match the supported classes of problem in FORCES Pro. Here, all computations are executed on a single desktop computer, with 2.60 GHz i7-4510U CPU and 6.00 GB RAM; however, in practice, the resulting independent sub-problems can be solved in parallel.
As mentioned before in the paper, in Voronoi-based methods, a vehicle only requires the position information from its neighboring vehicles to generate its trajectory. Therefore, they are more suitable for implementation when vehicles have limited communication capability, and have to rely solely on onboard sensing. In reality, the position sensor noise can impact the planner performance, yet this is more pronounced when estimating other information, such as velocity, from noise-corrupted measurements is needed. Therefore, Voronoi-based planners are more robust when there is no communication network. Nevertheless, in the following simulations, we assume that accurate position information is available with no delay at the replanning time.
In the first example, we consider five drones flying from their initial positions to given final positions. This example is similar to one in [24] where a random offset is added to break the symmetry in the drones' initial and final configurations. Figure 7 (right) shows collision-free trajectories generated with the distributed scheme described above, with a replanning rate of 20 Hz. For this particular example, the resulting trajectories match those generated with BVC with a flight time of 11.6782 s. Figure 7 (left) shows collision-free trajectories obtained from the centralized solution, which delivers a total flight time of 9.4347 s, yet, while the central solution is obtained in 601 ms, the average computation time for solving the sub-problems in the decentralized scheme is only 49 ms.
In the next example, we consider 18 drones switching positions in a 3 m × 5 m × 2 m space, with a maximum speed and acceleration of ±4.7 m s and ±9.8 m s 2 , respectively. Figure 8a shows the initial and final configurations, and Figure 8b displays collision-free trajectories generated with the proposed distributed scheme in the paper implemented at 10 Hz. While both methods could find collision-free trajectories for guiding the team of drones from their initial positions to their goal positions, the flight time achieved with the proposed method is markedly shorter than the time obtained with BVC. We also performed a trial simulation with 34 drones in a similar configuration. Table 1 compares the success rate and the flight time to complete the transition using BVC and the proposed method.
In the third example, we consider 100 drones flying in an 8 m × 8 m × 3.5 m space. The initial and final positions for the drones are displayed with dot and square markers in Figure 9. We test both methods in 30 different trials. In each trial, final positions are randomly assigned to drones. A trial is considered successful if all drones could reach their final positions within the stipulated time. The proposed method with 23 successful trials and an average total flight of 1s outperforms the BVC with only 16 completed trials. It should be noted that using well-devised deadlock prevention strategies or loosening time constraints can improve the success rate of both methods. Figure 9 shows collision-free trajectories generated with the proposed method for one of the trials at different time steps. The average computation time for solving sub-problems in this example was around 115 milliseconds. In addition, compared to the geometric algorithm in [24], the closest point in a Voronoi cell to the goal position was obtained at least 10 times faster with the proposed algorithm in Section 2.3. The computation time for finding the closest point, and solving the optimization problem, depends on the number of neighboring drones (See Table 2 for recorded computation times in simulation examples with 18, 34, and 100 drones). In most applications, with typical Voronoi diagrams, the number of boundary planes, i.e., the number of Voronoi neighbors, is small. Thus, the proposed distributed algorithm is scalable to arbitrary numbers of vehicles.

Conclusions
In this paper, we introduce an efficient distributed algorithm for generating collisionfree trajectories for multiple drones, taking into account their orientation. In order to avoid substantial communication between drones, we adopt Voronoi-based space partitioning and derive local constraints that guarantee collision avoidance with neighboring vehicles for an entire time horizon. We leverage Bézier curve properties to ensure that the set of collision avoidance constraints are satisfied at any time instant of the planning horizon. The same approach can be employed to obtain local collision avoidance constraints for the cases where the normal vector and offset of separating planes are time-varying parameters or decision variables of sub-problems. Yet, adopting Voronoi diagram with fixed planes for an entire planning horizon, though being conservative, results in simple, small sub-problems allowing for the trajectories to be replanned at a higher rate. We present different simulation results to highlight the scalability of the algorithm to large numbers of drones, and also its capability to generate less conservative trajectories with notably shorter flight times, compared to other Voronoi-based methods.
Our future work includes implementation and experimental validation of the algorithm for teams of drones. As we explain in the paper, at each time sample, upon receiving (or sensing) the new position information, a vehicle must find the closest point in its Voronoi cell to the goal position, and solve an optimization problem that uses the current state of the vehicle as the initial condition, to generate its trajectory for a certain time horizon. Although the time to compute the closest point is mainly negligible, the computation time to find the optimal solution can lead to a (significant) delay between updating the position information and executing the trajectory. Therefore, in practice, the computational delay must be explicitly considered to avoid performance degradation (or even failure) of the planner.