In recent years, robots have gradually found their way into the entertainment industry. Due to the advancements in both software and hardware, autonomous quadrotors, for example, have become active participants in art work and entertainment activities, showing off their exceptional agility in navigating in the three-dimensional space. Companies such as Intel and SKYMAGIC have delivered dazzling drone shows, where hundreds or even thousands of drones equipped with dedicated LEDs were coordinated to display enormous animations in the sky. Other companies, such as Verity Studios and ElevenPlay, have taken their shows to a different level, involving coordinated movements of drones with human performers, light effects, and music.
Drone shows where visual effects dominate the audience’s experience have been pushed to an unprecedented scale and flown as many as 2018 drones 
. However, shows that primarily rely on the effect of highly dynamic motions and fully leverage the motion capabilities of quadrotors remain at a much smaller scale. One reason is that drone swarms performing highly dynamical motions demand fast and accurate state estimates, motion planners and controllers, as well as an efficient and reliable overall flight control system architecture (e.g. with little or predictable time delays). Moreover, from a producer’s perspective, the design and choreography of attractive, highly dynamic, swarm-like motion patterns for a large troupe of drones becomes more challenging due to inter-agent collision constraints and spatial constraints, as well as aerodynamic effects. Therefore, our work aims to provide motion planning and control strategies that enable a quadrotor swarm to act as an integrated entity and display highly dynamic motion patterns, as shown in Fig. 1.
In the context of motion planning for complex task specifications, motion primitives have proven to be a successful framework. Motivated by creating theatrical performances, we aim to design a library of motion primitives, from which choreographers can then choose and combine motion primitives to form a coherent story. First, individual motion primitives must be designed. Rhythmic motion primitives for individual vehicles were proposed in . A set of representative behaviour descriptors in  was designed to enable online interaction with a small group of robots. However, a unified framework for encoding a library of motion primitives for swarms is lacking.
Next, to concatenate swarm motion primitives, we develop an algorithm for smooth and safe transitions. Many techniques in the robot planning and control literature already exist to transition vehicles from one configuration to another. Of particular interest is a simultaneous goal assignment and trajectory planning algorithm with a simple collision avoidance scheme proposed in , which has been applied in the context of drone performances . Despite its scalability to large swarms and its consideration of vehicle dynamics at an early stage in the planning process, this method does not immediately apply to our situation. Namely, it does not guarantee a common start and end time for all vehicles, and it assumes the vehicles are at rest at the start and end, which makes designing smooth trajectories difficult for our dynamic motions.
Finally, once motion primitives and their transitions are designed, they need to be executed reliably by the real system. To this end, a phase comparator along with a correction algorithm waswere proposed in  to eliminate the phase error in tracking the desired position trajectories. Moreover, a correction method was taken to compensate for the predictable delay in sensing and communication. In this work, we directly apply the correction method to quadrotors with standard position and attitude controllers. We demonstrate that this method can reliably minimize tracking errors of periodic signals for quadrotor swarms.
The contributions of this paper are as follows. First, we present a generic formulation for swarm motion primitives that is suitable for describing periodic, highly dynamic motion patterns, where drones appear as a coordinated unit embodying moving and deforming objects. Second, we adapt and combine various state-of-the-art trajectory planning methods to safely and smoothly connect arbitrary swarm motion primitives and, additionally, preserve patterns in the swarm behaviour that could emerge from strategic goal assignment. Third, we show that a simple correction algorithm is capable of compensating for amplitude and phase errors that arise when trying to synchronize a swarm’s motion to a desired periodic motion pattern. Finally, results are demonstrated experimentally on a large swarm of drones.
Ii Dynamic, Periodic Motion Patterns
This section presents our design of dynamic, periodic motion primitives for drone swarms. These motion primitives embody moving and deforming objects and are inspired by natural particle phenomena such as wave motions or rigid body rotations. The goal is that vehicles appear as an integrated entity.
Ii-a Generic Formulation
We define a motion primitive for a swarm of drones as a tuple
which consists of the start time and end time
of the motion primitive, a unique characteristic configuration vectorfor each drone , and a desired position trajectory generator . Then, the desired position evolution for the -th drone over the interval is given as .
As in , we parameterize the trajectory generator as a finite sum of periodic functions to encode rhythmic motions. Moreover, we use the configuration vectors to emphasize how each drone’s motion is related to the overall display. Specifically, trajectory function is defined as
with frequencies as well as parameter functions , and that define the centre and amplitude of the sinusoidal functions in 3D. Here indexes motions with different frequencies and parameter functions.
The parameters in (2) reflect desired rhythmic and visual effects. The temporal components describe intervals of repeated patterns. The spatial components , , , and are responsible for the spatial pattern collectively demonstrated by the swarm. In particular, the parameter functions determine the overall picture presented to the audience, where each drone’s contribution is reflected by their unique characteristic configuration vector.
Some possible choices for spatial patterns are given as below. In our experiments, we demonstrate two specific representations of (2), where vehicles represent the particle movement seen in waves and rigid body rotations.
Ii-B Wave Patterns
Wave motions are ubiquitous in nature. Typical examples include surface water waves, sound traveling through air, and electromagnetic field propagation. They can be modeled by the wave propagation equation . Consider a two-dimensional rectangular elastic surface with bounded edges as shown in Fig. 2a. Suppose the equilibrium location of each point on the surface is parameterized by . The disturbance to a point at time , , by a three-dimensional wave propagating through the surface at speed is governed by
where is the spatial Laplacian. Its solution can be expressed as a product of spatial and temporal components:
where the amplitudes can be determined from the surface’s initial position and velocity, and frequencies for are dictated by the dispersion relation .
Finally, we can use a finite approximation to this solution to generate reference trajectories for drones situated on
the surface. Suppose that we select pairs of terms in the finite approximation and assign a unique point to each drone on the surface at some desired height . Then the desired reference is for . Comparing to (1) and (2), we can design a swarm motion primitive, where the indices correspond to a pair , the frequencies are , and parameters are , and .
The surface wave (4) has desirable geometric properties, such as symmetry and periodicity in its spatial components , . However, the parameters , , , and must be carefully selected when representing a pattern on the continuous rectangular surface with a finite number of drones. To be specific, and determine the spatial frequency of the oscillation amplitude and thus, the symmetry and periodicity in the overall wave pattern. The pattern gets more interesting as the axes or points of symmetry increase until their spacing is smaller than that between drones. Although any distribution of the drones within the surface is valid and should not influence the overall pattern, a selection of the vectors that shares similar symmetry as , may offer a better visual experience to the audience, such as a radial or rectangular mesh. Examples illustrating the role of each parameter are shown in Fig. 2b and Fig. 2c.
Ii-C Rigid Body Rotation
A rigid body rotation can also be expressed in the form (2). Consider a rigid body rotating in an inertial frame at a constant angular velocity, for example, the tilted cone in Fig. 3. First we attach a body frame with origin relative to the inertial frame . Next, we define an inertially fixed frame which overlaps with at time .
Let be the position of a point expressed in . Then this point is expressed in the inertial frame as
where is the rotation matrix from to and, without loss of generality, is a principle rotation along the axis of given by . Expanding (5), we obtain the trajectory generator of a rigid body rotation:
This may now be compared with (2) to obtain the parameters and frequencies.
Since determines the periodicity and (, ) simply denotes the pose of the rotating objectis simply the translation of the rotating object, the characteristic vectors are the only parameters to be designed. Given the availability of only drones as well as the constraints on the minimum distance among them, should be strategically chosen in order to have the audience recognize the shape of the body. For instance, it is easier to identify a cone if the drones outline a helix lying on the surface of the cone, as opposed to spacing them uniformly.
Iii Transition Trajectory Planning
This section presents our transition trajectory planner that coordinates the swarm to make smooth and safe transitions and preserves geometric features in this process. We formulate the problem in Sec. III-A followed by our proposed method in Sec. III-B and III-C and justification in Sec. III-D.
Iii-a Problem Statement
Consider two consecutive motion primitives and defined as
In what follows, we assume a given common start time and end time for all drone actors, where and for some small numbers and . We aim to compute:
A common start time and end time for all drone actors that satisfy and for some small numbers and .
An assignment that assigns each drone identified by in a unique characteristic configuration vector in .
For each drone identified as in , a smooth trajectory from to that respects the quadrotors’ motion constraints, the flight space boundary and inter-agent collision constraints.
Our preliminary results showed that it is possible to optimize and for energy consumption. In this work, however, we assume and are given and focus on the goal assignment and trajectory generation problem.
We address the trajectory planning problem in two steps: (i) goal assignment and (ii) trajectory generation. As in [4, 5], we employ a collision resolution strategy in the second step, but we additionally enforce a common start and arrival time, as well as the smoothness of the resulting collision-free trajectories to achieve smooth and coordinated transitions.
Iii-B Goal Assignment
We formulated the goal assignment problem as a combinatorial optimization problem. In our application, the assignmentaims to maximize the smoothness of trajectories generated in the next step. More complicated swarm transition objectives can be defined, such as minimizing the likelihood of collisions. However, in this case, the cost function of one assignment depends on how other assignments are made, making it a challenging nonlinear optimization problem.
If we denote the cost of assigning to as , the mapping can be found by solving the linear assignment problem
using the Hungarian algorithm . To obtain the assignment cost , we solve a simplified minimum snap trajectory generation problem as in , for which only state continuity constraints are enforced. We then take the optimal cost function value of that problem as the assignment cost, namely, T_α,β(⋅)∫_t_s ^t_e (T_α,β^(4)(τ))^2 dτ J_a(r_1^α,r_2^β) = T_α,β^(p)(t_s) = T_d,1^(p)(r_1^α, t_s) T_α,β^(p)(t_e)=T_d,2^(p)(r_2^β, t_e), with . Note that is the transition trajectory assigning to parametrized with a single polynomial curve of order in each direction
where is the coefficient vector.
Following , the optimization problem in (7) can be written in quadratic form, given as x_α, β^n 12 x_α, β^n^T H x_α, β^n H(x_α, β^n, r_1^α, r_2^β)= 0, where is the hessian of the minimum snap cost function with respect to the decision variable and denotes the state continuity constraints in (7). Note that this problem is fast to solve, making it a suitable cost function for the linear assignment problem.
Iii-C Collision-Free Smooth Trajectory Generation
Given , and the assignment , we aim to find a smooth and collision-free trajectory for each drone. Inspired by , we decouple the problem into sub-problems to avoid accounting for the collision constraints in a large joint space. However, our proposed algorithm consists of two steps: (i) find a dynamically feasible candidate trajectory for each drone denoted as and (ii) iteratively resolve collisions in (if any) in a sequential manner to obtain . In what follows, we parametrize both candidate and collision-free trajectory as in (8) denoted as and respectively.
Iii-C1 Generating Candidate Trajectories
To generate candidate trajectories for the assignment , we solve the full minimum snap problem in  for each drone. This problem is the same as (III-B) if we let and but with an additional set of state constraints that bound each drone’s position and higher order states at each discrete time step . with , and ,
Iii-C2 Iterative Collision Resolution
Given the candidate trajectories , we construct a directed collision graph . A vertex represents the drone identified by , and an edge , points from to , indicating a collision between and . For any edge , we force the collision avoidance maneuver to be executed solely by drone , while drone follows its intended path. Our goal is to remove all edges by finding a collision-free trajectory for drone .
In particular, we find by penalizing its difference from with additional ellipsoid collision constraints, given as. Define and to be the set of equality constraints and inequality constraints in (LABEL:eq:SmthTrajGen), respectively. Then, the collision avoidance problem for drone is given as
x_s^n ∑_k=1^K e_k^n^TWe_k^n H(x_s^n, r_1^n,M(r_1^n) )= 0 F_k(x_s^n)⩽0, k=1,2,…,K. ——E^-1(T_s^n&(t_k) - T_c^m(t_k))——_2^2⩾2, where , specifies the ellipsoid collision boundary as in , and (, ) are as previously defined. The deviation of drone at time is given by . The weighting matrix is a positive definite diagonal matrix trading off the deviation in each direction. Note that we can write for any polynomial trajectory with a suitable matrix that depends only on . Therefore, the cost function and collision constraints in (III-C2) can be written in quadratic form in .
We solve (III-C2) sequentially for drones in in decreasing order of the amount of outbound edges. If an optimal solution is found for , we remove all of its outbound edges; otherwise, in cases where the problem is temporarily infeasible or hard to find a solution (e.g., other drones are not in favourable positions), we skip to the next drone. This procedure is repeated until either is empty or the maximum iteration is exceeded.
Iii-D Design Considerations
In order to generate smooth transition trajectories, a few key design decisions were made that differentiate our method from those in previous work. First, the vehicle dynamics and transition time are incorporated into the assignment cost function to facilitate smooth trajectory generation in the following steps. Second, trajectories are parameterized with a single polynomial curve instead of concatenated polynomials to minimize unnecessary curvature. Finally, in contrast to [4, 5], continuity in trajectories at and is guaranteed in the collision resolution step. We illustrate the first design feature in the top two figures of Fig. 4c, which highlights the difference between two different choice of , the minimum snap function proposed in this work and the Euclidean distance used in . The former assignment cost function induces smooth, fluid, and energy-saving candidate trajectories as compared to the latter.
Our choice of the cost function in (III-C2) allows us to preserve geometric features in the transition process. As an example, theThe bottom figure in Fig. 4c shows a human-designed goal assignment , where rotational symmetry of the initial and final motion primitives were incorporated. Consequently, it introduces nice geometric properties, for example rotation symmetry in this particular example, in the candidate trajectories that enable a coordinated transition process. However, incorporating these strategic assignments into our framework is non-trivial and would result in a nonlinear assignment problem. Nonetheless, if we do have particular features in the candidate trajectories, we are able to preserve them during the collision resolution step by penalizing the difference between the final and candidate trajectories in (III-C2). Moreover, we introduced the weighting matrix to optionally preserve the swarm transition patterns in some dimensions while relaxing them in the others. An example is shown in Fig. 4b, where patterns in the plane are preserved during collision avoidance by encouraging motion in the direction.
Iv Motion Synchronization
To execute the proposed highly-dynamic motion patterns in tight formations and in sync with the desired periodic motion pattern, a high-accuracy controller for periodic motions is required. In practice, one major issue that was observed is the phase shift and amplitude amplification or attenuation  due to the delay in communicating vehicle commands, the inherent delay in dynamical systems, as well as the delay in sensor measurements. Although it is possible to reduce delay and amplitude error by adding feed-forward terms to the drones’ position and attitude controllers, the communication overhead it incurs is undesirable for a large swarm. Similarly to , we actively adjust the desired trajectories in (2) by scaling their amplitude and shifting them in time to obtain a new reference trajectory.
If we approximate the quadrotor as a linear system in each spatial direction with transfer functions , , we can estimate the amplitude attenuation and phase shift at each oscillating frequency from the system’s frequency response to sinusoidal reference signals. The compensated reference trajectory in the direction is given as
where the amplitude scalings and phase shifts can be determined from the closed-loop system’s transfer function observed in experiments, . That is, we compute , . We only compensate for the motion primitives but not for the transition trajectories since the latter has an infinitely wide frequency spectrum while the former’s is finite.
V Experimental Results
We demonstrate our algorithms using the small-sized quadrotor platform Crazyflie 2.0 in a multi-agent testbed at the Dynamic Systems Lab (inspired by the Crazyswarm ). From a central computer, we gathered position data from all the drones in the fleet using an overhead motion capture system. The state information was sent via radios to each drone’s onboard computer, along with the desired position, which is tracked using an onboard position controller.
V-a Motion Synchronization
We empirically determined the Crazyflie’s closed-loop transfer functions , , by commanding sinusoidal position trajectories at different frequencies in the and directions. Based on the phase shift and amplitude attenuation between the desired and actual trajectories, we found the system’s frequency response characterized by amplitude attenuation and phase shift as shown in Fig. 5. Based on the data, we constructed a look-up table from which the scaling factors for motion primitives, and
, are determined by linearly interpolating points in the magnitude and phase mapping, respectively. We constructed the look-up table using data from one drone but evaluated it on another 8 drones. As shown in Fig.6
, the tracking performance was greatly improved with negligible variances among the 8 drones being testedFig.6 summarizes the improvements on tracking performance of the 8 drones compared to when there is no synchronization.
V-B Transition Trajectory Planning
The optimization problems in trajectory planner are solved using the nonlinear optimizer IPOPT . We used as the polynomial order and as a starting point for partition intervals, which may be doubled in the next iteration if a feasible solution is found but a collision occurs in between two time steps and . To avoid numerical issues, we nondimensionalized the minimum snap cost function as in .
We evaluated the transition trajectory planner in both simulation and experiments. In simulation, we tested with 25 drones transitioning in a volume with a collision radius of in the - plane and in the direction. Our algorithm found a feasible solution in of 1800 randomly generated motion primitive pairs. The main reason for failure are numerical difficulties in solving (III-C2) when the transition time is long. One repair strategy is to parameterize the trajectories with a few concatenated polynomials, instead of just one polynomial.
V-C Choreographed Drone Performance
The video at http://tiny.cc/fast-periodic presents a one-minute drone performance where 25 drones are choreographed to fly in tight formations and perform fast and dynamic motions. We seamlessly concatenate five wave motions and one rigid body rotation to create a cohesive performance. The peak velocity and pitch angle reached and degrees, respectively. Since the motion primitives are aggressive, we used a larger collision ellipsoid than in , with in and in , which are the radii of the smallest ellipse in - plane that two drones can trace at a moderate angular velocity with phase shift. These parameters were coarsely estimated and worked well in our experiments. However, to fully exploit the spatial and the vehicles’ physical limits, more sophisticated methods should be used to explicitly model the rotors’ aerodynamics.
In this paper, we provide guidelines for creating performances with quadrotor swarms that fully leverage their motion capabilities to create appealing visual effects. The swarm motion primitives are formulated as coupled periodic motions, which are described by a single equation indicating the overall motion pattern and the relationship between the individual actors of the swarm. The geometric properties of parameter functions in this formulation were discussed. Moreover, we provided a hierarchical transition trajectory planner to seamlessly connect motion primitives together and preserve geometric characteristics. Lastly, a correction algorithm is proposed to improve the quadrotors’ tracking performance of the periodic motions, which allows the swarm to perform synchronously and in close proximity to each other. The method is validated with a swarm performance of 25 drones in a compact space.
-  Intel. (2018) Intel breaks guinness world records title for drone light shows in celebration of 50th anniversary. [Online]. Available: https://newsroom.intel.com/news/intel-breaks-guinness-world-records-title-drone-light-shows-celebration-50th-anniversary/.
-  F. Augugliaro, A. P. Schoellig, and R. D’Andrea, “Methods for designing and executing an aerial dance choreography,” IEEE Robot. Autom. Mag, vol. 20, no. 4, pp. 96–104, 2013.
-  E. A. Cappo, A. Desai, M. Collins, and N. Michael, “Online planning for human–multi-robot interactive theatrical performance,” Autonomous Robots, pp. 1–16, 2018.
-  M. Turpin, K. Mohta, N. Michael, and V. Kumar, “Goal assignment and trajectory planning for large teams of interchangeable robots,” Autonomous Robots, vol. 37, no. 4, pp. 401–415, 2014.
-  A. Desai, E. A. Cappo, and N. Michael, “Dynamically feasible and safe shape transitions for teams of aerial robots,” in IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 2016, pp. 5489–5494.
-  A. Schöllig, F. Augugliaro, S. Lupashin, and R. D’Andrea, “Synchronizing the motion of a quadrocopter to music,” in IEEE International Conference on Robotics and Automation (ICRA), 2010, pp. 3355–3360.
-  H. Georgi, The physics of waves. Prentice Hall Englewood Cliffs, NJ, 1993.
-  H. W. Kuhn, “The hungarian method for the assignment problem,” Naval research logistics quarterly, vol. 2, no. 1-2, pp. 83–97, 1955.
-  D. Mellinger and V. Kumar, “Minimum snap trajectory generation and control for quadrotors,” in IEEE International Conference on Robotics and Automation (ICRA), 2011, pp. 2520–2525.
-  C. Richter, A. Bry, and N. Roy, “Polynomial trajectory planning for aggressive quadrotor flight in dense indoor environments,” in Robotics Research. Springer, 2016, pp. 649–666.
-  Y. Chen, M. Cutler, and J. P. How, “Decoupled multiagent path planning via incremental sequential convex programming,” in IEEE International Conference on Robotics and Automation (ICRA), 2015, pp. 5954–5961.
-  J. A. Preiss, W. Hönig, N. Ayanian, and G. S. Sukhatme, “Downwash-aware trajectory planning for large quadrotor teams,” in IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), 2017, pp. 250–257.
-  J. A. Preiss, W. Honig, G. S. Sukhatme, and N. Ayanian, “Crazyswarm: A large nano-quadcopter swarm,” in IEEE International Conference on Robotics and Automation (ICRA), 2017, pp. 3299–3304.
-  A. Wächter and L. T. Biegler, “On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming,” Mathematical programming, vol. 106, no. 1, pp. 25–57, 2006.