Robust and Efficient Trajectory Planning for
Formation Flight in Dense Environments

Abstract

Formation flight has a vast potential for aerial robot swarms in various applications. However, existing methods lack the capability to achieve fully autonomous large-scale formation flight in dense environments. To bridge the gap, we present a complete formation flight system that effectively integrates real-world constraints into aerial formation navigation. This paper proposes a differentiable graph-based metric to quantify the overall similarity error between formations. This metric is invariant to rotation, translation, and scaling, providing more freedom for formation coordination. We design a distributed trajectory optimization framework that considers formation similarity, obstacle avoidance, and dynamic feasibility. The optimization is decoupled to make large-scale formation flights computationally feasible. To improve the elasticity of formation navigation in highly constrained scenes, we present a swarm reorganization method that adaptively adjusts the formation parameters and task assignments by generating local navigation goals. A novel swarm agreement strategy called global-remap-local-replan and a formation-level path planner is proposed in this work to coordinate the global planning and local trajectory optimizations. To validate the proposed method, we design comprehensive benchmarks and simulations with other cutting-edge works in terms of adaptability, predictability, elasticity, resilience, and efficiency. Finally, integrated with palm-sized swarm platforms with onboard computers and sensors, the proposed method demonstrates its efficiency and robustness by achieving the largest scale formation flight in dense outdoor environments.

I Introduction

Formation flight has become a fundamental capability for autonomous swarms to achieve coordinated aerial maneuvers. In cluttered wilds and complex urban areas, formation navigation has a wide potential in search and rescue[1], collaborative mapping[2], package delivery[3], and so on. However, effectively integrating real-world constraints into an aerial formation remains an unresolved problem. This article aims to empower aerial swarms to maintain cooperative formation behaviors in dense environments by proposing a complete formation flight system.

Inspired by natural swarm systems like bird flocks and fish schools, an ideal formation flight system should possess the capability to flexibly adapt and deform in dense environments. By striving to maintain the swarm in a “critical state”, swarm formation can dynamically balance the conflicts between maintaining formation and avoiding obstacles.

Figure 1: Large-scale formation flight in the dense outdoor environment. (a) Snapshot of the moment the swarm robots prepare to fly through the woods. (b) Rviz diagram of the executed trajectories. The grid map is merged using log data offline. Please watch our attached videos for more information about the experiments at https://www.youtube.com/watch?v=uEMyvPxYqmA.

While extensive research works focus on navigation in formation, few achieve robust formation flights in obstacle-rich areas. Three core challenges limit practical formation applications: (a) The inherent conflict between formation maintenance and obstacle avoidance is inevitable and difficult to mitigate. (b) Predefined formations lack elastic adaptability in response to constrained environments. (c) The swarm system cannot rapidly recover from disordered states caused by unknown obstacles or sudden changes in the desired formation shape.

Based on the above challenges, we conclude that an ideal formation flight system should have the ability to maintain formation while avoiding obstacles, adjust swarm formation distributions according to constrained environments, and reorganize the formation quickly after emergencies. These characteristics are summarized as the PAPER criteria:

• Portability: Aerial robotic swarms should comprise lightweight platforms with scalable systems and distributed architecture. A scalable system means the main components, such as estimation, decision, planning, and control modules, are all the same on each robot. A distributed architecture is inherently robust against individual hardware failures. These are the basis for large-scale formation flight.

• Adaptability: When facing obstacles, robots should locally adapt their trajectories to avoid collision in a way that does the least damage to the overall formation performance. This ability mitigates the conflict between formation maintenance and obstacle avoidance.

• Predictability: Reactive local feedback methods are short-sighted and can not consider the constraints in advance. Robots should optimize the motions over a prediction horizon so that the formation can respond smoothly to the future environmental changes in its vicinity, which is necessary for dense areas.

• Elasticity: A feasible and safe trajectory for a fixed formation shape may not exist in constrained environments, such as narrow corridors or holes. Therefore, swarm robots need to have elastic and flexible deformation capabilities by adjusting formation distributions (such as the scale of shape or task assignments) while keeping the full maneuverability of the formation.

• Resilience: Formation flight could encounter many unfavorable situations caused by unknown obstacles or sudden changes of desired formation shape. The navigation system should be able to resiliently reorganize and guide the whole formation so that the flight can recover from disordered states timely.

A complete formation flight system should meet the above PAPER criteria and ensure that the conditions for each criterion are compatible with the others.

Our previous work[4] only partially met the first three terms of PAPER criteria. We tackled formation flight as a coupled collaborative trajectory optimization problem, suitable primarily for small-scale formation scenarios. However, resolving cooperative constraints of the formation using the graph-based similarity metric was computationally heavy, resulting in increased overhead during each optimization iteration. Moreover, integrating dynamic inter-robot relationships within the coupled trajectory optimization problem considerably affected the efficiency of the optimization process, rendering it less appropriate for larger formations or more complex scenarios.

In this paper, we present a complete formation flight system that satisfies all PAPER criteria. To address the challenges in [4], we introduce a decoupled formation optimization method to significantly improve computational efficiency. This method consists of two components. Firstly, an optimal formation position sequence is pre-computed, avoiding repetitive metric calculation during the optimization process. Secondly, a fixed time interval sampling method is used to convert dynamic inter-robot relationships into static constraints, greatly reducing the complexity of the optimization problem. These improvements make our method suitable for large-scale swarms. Besides, the previous method lacks the capability of reorganizing the swarm formation, which may lead to disordered formation flights under adverse conditions, especially when the initial positions or task assignments are inappropriate. To address this, we propose a swarm reorganization method that can elastically adjust formation distributions by optimizing formation parameters and task assignments in response to external constraints. Subsequently, we develop a swarm agreement strategy called global-remap-local-replan, which enables rapid implementation of the swarm reorganization results to achieve consensus among swarm agents. Additionally, a formation-level global path-finding method, which treats the swarm formation as one entity, is also designed to guide the swarm out of the obstacle deadlocks. Finally, we integrate the estimation, mapping, decision, planning, and control modules into palm-sized swarm platforms [5] with onboard computers and sensors, enabling large-scale formation flight in dense environments. Detailed contributions are as follows.

1. We introduce an optimal formation position sequence, pre-computed using the differentiable graph-based metric [4]. This sequence represents the optimal position with the lowest similarity error, reducing the need for repetitive computation during the optimization process.

2. We design a decoupled spatial-temporal trajectory optimization framework that effectively handles dynamic inter-robot relationships, obstacle avoidance, and dynamic feasibility. Compared to our prior work [4], we achieve higher computational efficiency for large-scale swarms.

3. We present a swarm reorganization method to achieve elastic deformation of swarm distributions, which simultaneously solves optimal formation alignment and task assignment problems (ALAS for short). This method improves the elasticity of swarm formation against constrained environments. It relieves the dependence on the appropriate formation alignments and task assignments.

4. We design a global-remap-local-replan strategy (GRLR for short) that leverages the advantages of centralized formation parameter remapping and decentralized local trajectory replanning. With this strategy, the distributed asynchronous swarm is able to quickly recover from disordered states and return to formation flight quickly.

5. We integrate all these modules into a hierarchical formation flight system. Extensive benchmarks and simulations are conducted to validate the PAPER criteria of our method. A series of real-world experiments are designed to demonstrate the outstanding performance of the proposed distributed autonomous formation flight system.

II Related works

II-A Distributed Swarm Trajectory Planning

Extensive works exist for trajectory planning of distributed swarms. The concept of velocity obstacle (VO) is leveraged and generalized by Van Den Berg et al.[6, 7, 8] to accomplish reciprocal collision avoidance for multiple robots. However, the smoothness of the resulting trajectories cannot be guaranteed by VO-based approaches, which significantly impairs the usability of the actual robot systems.

In order to produce high-quality collision-free trajectories, optimization-based methods are widely introduced in the literature on distributed multicopter swarms[9, 10, 11]. Zhou et al.[12] incorporate Voronoi cell tessellation into a receding horizon QP scheme to prevent collision among the robots while planning. In [13], Chen et al. employ SCP to address the multiagent planning problem in non-convex space by incrementally tightening the collision constraints. Baca et al.[14] combine MPC with a conflict resolution strategy to ensure mutual collision avoidance for outdoor swarm operations. Nevertheless, the computational load of the above optimization-based methods is large, which could hamper the applicability of the planners in highly dense scenarios.

Recently, Zhou et al.[5] present a distributed autonomous quadrotor swarm system using spatial-temporal trajectory optimization, which generates collision-free motions in dense environments merely in milliseconds. Our distributed formation trajectory optimization is based on this work.

II-B Formation Flight in Free Space

Various techniques have been proposed to achieve multi-robot navigation in formation, which include virtual structures[15], leader-follower[16], navigation functions[17], reactive behaviors[18], consensus-based local control laws[19], and barycentric-coordinate-based control [20]. However, most of the existing methods only consider obstacle-free cases.

Weinstein et al.[21] present a VIO-swarm system that performs all modules onboard and can execute formation flight without inter-robot collisions in free space. Parker et al.[22] present a distributed formation control method and relax the dependency of the common reference frame.

As the scale of swarms increases, researchers begin to notice that it is difficult to maintain the formation only by trajectory planning, especially when there are deadlocks between robots. Turpin et al.[23] consider the problem of concurrent assignment and collision-free trajectory generation. Turpin gives centralized and decentralized solutions to this problem, allowing flight formation on a large scale. Morgan et al.[24] also use model predictive control to solve task assignment and trajectory generation simultaneously when given the desired formation shape. In addition to considering task assignments, Agarwal and Akella[25] consider formation alignment problems to optimize the formation parameters such as scale and location. This method reduces the cost of forming formation and speeds up convergence. However, these methods ignore the influence of constrained environments, in which formation should elastically deform to navigate.

II-C Formation Flight in constrained Environments

In constrained environments, where various obstacles and limitations exist, formation flight can be a challenging task that requires constant adjustments to maintain the swarm structure. An immediate solution is to design composite control laws that combine formation flight and collision avoidance by using multiple layered potential fields[26], which are prone to deadlock. A better solution is to allow the formation shape to deform while maintaining the overall swarm structure. Han et al. [27] propose a complex-valued graph Laplacian-based formation controller that regulates the scaling of formation shape during swarm maneuvering like passing through corridors. In [28], Zhao proposes a leader-follower control law enabling the affine transformation of formation in response to environmental changes. And the bearing-based local controller [29, 30] exhibits translational, scaling, and rotational invariance of formation flight. However, these methods rely on leaders or predefined trajectories and struggle with complex obstacles or sudden potential collisions.

Compared to the local feedback methods, predictive optimization-based methods proactively plan the future motion of swarm robots, striking a balance between formation flight and obstacle avoidance. Alonso-Mora et al.[31] control swarm robots by optimally rearranging the desired formation and planning local trajectories for each drone. However, since there is no inter-vehicle coordination in the distributed planners, formation maintenance is not conducted during local planning. Peng et al.[32] propose a method to improve flight safety by enabling the affine transformation of formation shape and treating it as a soft constraint during B-spline optimization. However, this approach requires optimizing the trajectories of all robots simultaneously and cannot be applied to large-scale swarms. To tackle formation preservation, Parys et al.[33] propose a distributed model predictive formation controller. This framework imposes relative position constraints on the swarm and coordinates the agents to break passively once obstacles violate positional constraints. Overall, these approaches offer unique solutions for trajectory planning in swarm robotics, but they each have limitations when dealing with different scenarios and scales of robotic systems.

To address these drawbacks, we formulate the overall formation requirement with a differentiable metric in trajectory optimization. This allows us to fully utilize the collaboration ability of the swarm, effectively avoid deadlock, and foresee obstacle avoidance. Besides, we adopt a distributed and decoupled optimization method to ensure dynamic real-time performance. This approach can be applied to large-scale swarms while still maintaining efficient trajectory planning.

Figure 2: Illustration of the system architecture. In order to facilitate understanding, we divide the various modules of the formation flight system according to the PAPER criteria mentioned in Sec.I. The main challenge in integrating PAPER criteria into a swarm formation flight system is to ensure that the conditions for each criterion are compatible with the others. Therefore we build a hierarchical formation flight system. We use different colors to represent different levels in the system. Moreover, the information broadcast through the network is from the same color module.

III System Overview

This paper aims to optimize the autonomy of swarm robots in real-world environments by coordinating their movements to form a desired formation shape. To accomplish this, we adopt a distributed swarm aerial robot system and propose a spatial-temporal trajectory optimization for formation flight. To enhance the system’s robustness, we also address the case of swarm disorder by incorporating an adaptive swarm reorganization method and an efficient swarm agreement strategy.

III-A Swarm Aerial Robot System

The swarm aerial robot system is composed of palm-sized quadrotor platforms[5] with depth stereo camera for imagery and depth sensing, as shown in Fig. 2. The software modules, including state estimation, environment perception, decision-making, trajectory planning, and flight control, run in real-time on an onboard computer. This lightweight and scalable platform suit dense environments.

We use visual-inertial odometry (VIO)[34] to estimate each robot’s pose with respect to its start frame, and we recover the transformations related to the start frames with only anonymous bearing measurements in our previous work[35]. For simplicity, the transformations are known in our experiments by requiring robots to take off from pre-defined locations. To correct the localization drift between swarm robots, we use a drift correction method[5] with onboard ultra-wideband (UWB).

The distributed system architecture enables each robot to fully utilize its computing resources to process more information, relieving the pressure of network communication. The robots share only important information, such as trajectories, for high-fidelity wireless communication, and there is no ground station to send control inputs.

III-B Distributed Local Formation Trajectory Optimization

The local formation trajectory optimization is distributed and asynchronous, which enables each robot to generate its trajectory only depending on local information and does not require the same start timestamp or same time duration of trajectory. Each robot evaluates the formation state by calculating the formation similarity error and generates its optimal formation position sequence (Sec.IV) to maintain the desired formation shape. Then a subsequent trajectory optimization module generates spatial-temporal formation trajectories (Sec.V) for real-time navigation. The above process cycles periodically within the receding horizon. The distributed local formation trajectory optimization can always maintain the overall formation when navigating in a complex environment.

III-C Swarm Reorganization and Agreement Methods

Only relying on the distributed local trajectory optimization may lead to poor formation flight quality when encountering narrow corridors or instant transformation of formation shapes. Therefore, we propose swarm reorganization (Sec.VI) and agreement methods (Sec.VII) to achieve the swarm consensus quickly. Inspired by the flock flying behavior of birds, we design a global-remap-local-replan (GRLR for short) strategy, which only centrally remaps crucial parameters of swarm formation and makes the swarm formation converge quickly through distributed replanning local trajectory. Firstly, when the stable state of the swarm formation is destroyed or about to be destroyed, the swarm robots designate a leader (drone 0 in this paper). The crucial parameters of swarm formation are calculated separately by formation alignment and task assignment optimization (ALAS for short) and formation-level global pathfinding method. From the perspective of swarm reorganization and agreement, the conflict between swarm formation and obstacle is alleviated by optimizing formation alignment, and the speed of formation convergence is greatly accelerated by optimizing task assignment. GRLR strategy is simple but very effective by combining distributed methods’ efficiency and centralized ones’ optimality.

IV adaptive description of swarm formation

IV-A Graph-based Formation Definition

In this paper, a swarm formation of $N$ robots is modeled by an undirected graph $𝒢=(𝒱,ℰ)$, where $𝒱:=\{1,2,\mathrm{\dots },N\}$ is the set of vertices, and $ℰ\subset 𝒱\times 𝒱$ is the set of edges. In graph $𝒢$, the vertex $i$ represents the $i^{th}$ robot with position vector ${𝐩}_i=[x_i,y_i,z_i]\in {ℝ}^3$ . An edge $e_{ij}\in ℰ$ that connects vertex $i\in 𝒱$ and vertex $j\in 𝒱$ means that robot $i$ and $j$ can measure the geometric distance between each other. In our work, each robot can obtain the positions of all robots $\{{𝐩}_1,\mathrm{\dots },{𝐩}_i,\mathrm{\dots },{𝐩}_N\}$, thus the graph $𝒢$ is complete. Then we determine the adjacency matrix $𝐀\in {ℝ}^{N\times N}$ and degree matrix $𝐃\in {ℝ}^{N\times N}$ of the formation graph $𝒢$ by

$$A_{ij}=w_{ij}={\Vert {𝐩}_i-{𝐩}_j\Vert }^2,$$ (1)

(2)

where the non-negative edge weight $w_{ij}$ is the squared distance between the $i^{th}$ and $j^{th}$ robots, and $\parallel \cdot \parallel$ denotes the Euclidean norm. Thus, the corresponding Laplacian matrix is

$$𝐋=𝐃-𝐀.$$ (3)

With the above matrices, the symmetric normalized Laplacian matrix of graph $𝒢$ is defined as

$$\widehat{𝐋}={𝐃}^{-1/2}{\mathrm{𝐋𝐃}}^{-1/2}=𝐈-{𝐃}^{-1/2}{\mathrm{𝐀𝐃}}^{-1/2},$$ (4)

where $𝐈\in {ℝ}^{N\times N}$ is the identity matrix. $\widehat{𝐋}$ contains the information that is invariant to scale, translation, and rotation.

Finally, we use graph theory to describe various desired formation shapes, such as squares, hexagons, and pyramids. By specifying the positions ${𝐩}_i^d=[x_i^d,y_i^d,z_i^d]\in {ℝ}^3,i=1,\mathrm{\dots },N$, computing ${\widehat{𝐋}}_{des}$ is simple. It’s important to note that the desired formation shape is independent of the coordinate system as long as the relative positions are provided.

IV-B Differentiable Formation Similarity Error Metric

Figure 3: Illustration of optimal formation position sequence using a 2D formation. (a) The surface shows the profile of the similarity metric when one UAV moves in the plane and the other three remain still. The minimum suggests the optimal formation position to form the desired shape. (b) The sequence of optimal formation positions corresponds to the timestamps.

To assess the deviation from the desired formation, we propose a differentiable formation similarity error metric as

$f_s$ $=f_s({𝐩}_1,\mathrm{\dots },{𝐩}_i,\mathrm{\dots },{𝐩}_N)=f_s(𝐀,𝐃)=f_s(\widehat{𝐋},{\widehat{𝐋}}_{des})$ (5)

$={\Vert \widehat{𝐋}-{\widehat{𝐋}}_{des}\Vert }_F^2=tr\{{(\widehat{𝐋}-{\widehat{𝐋}}_{des})}^T(\widehat{𝐋}-{\widehat{𝐋}}_{des})\},$

where $tr\{\cdot \}$ denotes the trace of a matrix, $\widehat{𝐋}$ is the symmetric normalized Laplacian of the current swarm formation, ${\widehat{𝐋}}_{des}$ is the counterpart of the desired formation. Frobenius norm $\parallel \cdot {\parallel }_F$ is used in our distance metric. As a graph representation matrix, $\widehat{𝐋}$ contains information about the graph structure[36]. This allows $f_s$ to consider only the geometric shape of the formation, and not be influenced by scaling, translation, or rotation. Additionally, $f_s$ is a dimensionless value that solely reflects the error in formation shape similarity.

In particular, under the distributed framework, each robot can only change its positions to reduce the overall formation similarity error. Therefore, the only variable for robot $i$ in (5) is ${𝐩}_i$, and $f_s({𝐩}_1,\mathrm{\dots },{𝐩}_i,\mathrm{\dots },{𝐩}_N)$ can be simplified as $f_s({𝐩}_i)$.

Our metric is analytically differentiable with respect to the position of each robot. For robot $i$, we use the weights of its $N$ adjacent edges $\{e_{i1},e_{i2},\mathrm{\dots },,e_{iN}\}$ to form a weight vector ${𝐰}_i={[w_{i1},w_{i2},\mathrm{\dots },,w_{iN}]}^T$. By the chain rule, the gradient of $f_s$ with respect to ${𝐩}_i$ can be written as

$$\frac{\partial f_s}{\partial {𝐩}_i}=\frac{\partial f_s}{\partial {𝐰}_i^T}\frac{\partial {𝐰}_i}{\partial {𝐩}_i}.$$ (6)

According to our metric (5), the gradient of $f_s$ with respect to each weight $w_{ij}$ can be computed as follow

${\displaystyle \frac{\partial f_s}{\partial w_{ij}}}$ $=tr\{{({\displaystyle \frac{\partial f_s}{\partial \widehat{𝐋}}})}^T({\displaystyle \frac{\partial \widehat{𝐋}}{\partial w_{ij}}})\},$ (7)

${\displaystyle \frac{\partial f_s}{\partial \widehat{𝐋}}}$ $={\displaystyle \frac{\partial {\Vert \widehat{𝐋}-{\widehat{𝐋}}_{des}\Vert }_F^2}{\partial \widehat{𝐋}}}=2(\widehat{𝐋}-{\widehat{𝐋}}_{des}),$

${\displaystyle \frac{\partial \widehat{𝐋}}{\partial w_{ij}}}$ $=-{\displaystyle \frac{\partial ({𝐃}^{-1/2}{\mathrm{𝐀𝐃}}^{-1/2})}{\partial w_{ij}}}.$

Then the gradient $\partial f_s/\partial {𝐰}_i$ can be written as

$$\partial f_s/\partial {𝐰}_i={[\partial f_s/\partial w_{i1},\partial f_s/\partial w_{i2},\mathrm{\dots },\partial f_s/\partial w_{iN}]}^T.$$ (8)

As for $\partial {𝐰}_i/\partial {𝐩}_i$, the Jacobian can be easily derived since the weight function (1) is a differentiable quadratic form.

IV-C Optimal Formation Position Sequence

In our previous work [4], we incorporated $f_s$ directly into the trajectory optimization, making formation flight a coupled trajectory optimization problem. While this method is suitable for small-scale formation flight, it becomes computationally inefficient as the number $N$ of robots increases. Considering the simplified equation for coupled trajectory optimization

$$\underset{{𝐩}_{i,0},\mathrm{\dots },{𝐩}_{i,M_c}}{\min }\sum _{j=0}^{M_c}{f}_s({𝐩}_{i,j})+J_{other},$$ (9)

where ${𝐩}_{i,j}$ represent the $j^{th}$ sample point of $i^{th}$ robot trajectory in (19) for convenience. $J_{other}$ represents all other cost functions, and $M_c$ is the number of sample points with corresponding timestamps. The primary purpose of calculating $f_s$ is to supply gradient information for minimizing formation similarity error. However, since the graph $𝒢$ is a complete graph, computing $f_s$ has a complexity of $O(N^2)$. Consequently, the coupled trajectory optimization (9) also exhibits high complexity of $O(N^2)$ in each iteration, limiting its applicability to large-scale swarm operations.

To address this issue, we must identify an equivalent approach with reduced computational complexity to replace the function of $f_s$ in (9). We introduce the concept of optimal formation position ${𝐩}_{i,j}^{\ast }$ for robot $i$ at timestamp $j$, which is the position that minimizes the formation similarity error $f_s$. Fig. 3(a) illustrates this concept using a 2D formation as an example. It is evident from the figure that there exists an optimal formation position ${𝐩}_{i,j}^{\ast }$ that results in a minimal formation similarity error, and the partial derivative is $\partial f_s/\partial {𝐩}_{i,j}^{\ast }=0$. In the future period with a sequence of timestamps $\{0,\mathrm{\dots },j,\mathrm{\dots },M_c\}$, we represent the expected positions of robot $i$ with the optimal formation position sequence ${𝐩}_i^{\ast }=\{{𝐩}_{i,0}^{\ast },\mathrm{\cdots },{𝐩}_{i,j}^{\ast },\mathrm{\dots },{𝐩}_{i,M_c}^{\ast }\}$, as shown in Fig. 3(b). By precomputing ${𝐩}_i^{\ast }$, we can utilize its quadratic distance to replace the gradient information offered by $f_s$ in (11), thus decreasing the computational requirements as follows

$$f_s({𝐩}_{i.j})\Rightarrow {\Vert {𝐩}_{i,j}-{𝐩}_{i.j}^{\ast }\Vert }^2.$$ (10)

Since the optimal solutions of $f_s$ and quadratic distance cost are equivalent, the trajectory approaches the positions with minimal formation similarity error, maintaining the desired formation. Thus, we can effectively solve the coupled trajectory optimization with a two-step procedure

$\text{①}{𝐩}_i^{\ast }=\arg \min {\displaystyle \sum _{j=0}^{M_c}}{f}_s({𝐩}_{i.j}),$ (11)

$\stackrel{{𝐩}_i^{\ast }}{\Rightarrow }$ $\text{②}\underset{{𝐩}_{i,0},\mathrm{\dots },{𝐩}_{i,M_c}}{\min }{\Vert {𝐩}_{i.j}-{𝐩}_{i.j}^{\ast }\Vert }^2+J_{other}.$

As a result, the previously required calculation of $f_s$ in each trajectory optimization process is replaced by the computation of the quadratic distance, simplifying the optimization problem. This significantly reduces computational demands and enables large-scale swarm formation.

Formula (11) indicates that trajectory optimization in Sec.V is performed on discretized points. Non-uniform discretized points may lead to poor trajectories and sub-optimal performance. Therefore it is crucial to ensure a uniform distribution of these points to maintain the effectiveness of the optimization process. In engineering practice, since graphs $𝒢$ are constructed from a series of discretized timestamps as depicted in Fig. 3(b), each ${𝐩}_{i,j}^{\ast }$ is independent.

To ensure a smoother trajectory, we introduce the uniform optimal formation position sequence ${\widehat{𝐩}}_i^{\ast }$, which is generated by considering the formation similarity error $J_s$ and the uniform distribution cost $J_u$

$${\widehat{𝐩}}_i^{\ast }=\arg \min {\lambda }_s{J}_s+{\lambda }_u{J}_u,$$ (12)

$J_s$ $={\displaystyle \sum _{j=0}^{M_c}}{f}_s({\widehat{𝐩}}_{i,j}^{\ast }),$ (13)

$J_u$ $=𝔼({𝐔}^2)-𝔼{(𝐔)}^2={\displaystyle \frac{{\Vert 𝐔\Vert }_2^2}{M_c}}-{\displaystyle \frac{{\Vert 𝐔\Vert }_1^2}{{(M_c)}^2}},$

where ${\lambda }_s$ and ${\lambda }_u$ are the relative weights. $𝔼(\cdot )$ is mathematic expectation and the squared distance vector $𝐔\in {ℝ}^{M_c}$ is

$$𝐔=({\Vert {\widehat{𝐩}}_{i,1}^{\ast }-{\widehat{𝐩}}_{i,0}^{\ast }\Vert }_2^2,\mathrm{\cdots },{\Vert {\widehat{𝐩}}_{i,M_c}^{\ast }-{\widehat{𝐩}}_{i,M_c-1}^{\ast }\Vert }_2^2).$$ (14)

We use the quasi-Newton method [37] to solve this unconstrained optimization problem (12) and generate uniform ${\widehat{𝐩}}_i^{\ast }$ for the later trajectory optimization (18). By doing so, the trajectory resulting from these discretized points in Sec.V can be smoother and avoid sudden spatial changes.

V Spatial-temporal trajectory optimization for formation flight

V-A Trajectory Representation

The differential flatness of multicopters [38] benefits trajectory generation without integrating differential equations. Moreover, the motion planning of multicopters can be performed on low-dimensional smooth trajectories. In this paper, we adopt a state-of-the-art trajectory representation named MINCO [39] to achieve minimum control effort spatial-temporal trajectory planning for swarm aerial robots in three-dimensional environments. MINCO conducts spatial-temporal deformation of the flat-output $M$-piece trajectory $p(t)$ by decoupling the space and time parameters with a linear-complexity mapping $ℳ$

$$p(t)={ℳ}_{𝐪,𝐓}(t),\forall t\in [t_0,t_M],$$ (15)

where $𝐪={({𝐪}_1,\mathrm{\cdots },{𝐪}_{M-1})}^T\in {ℝ}^{3\times (M-1)}$ are the adjacent intermediate points between each pair of connected pieces and $𝐓={(T_1,\mathrm{\cdots },T_M)}^T\in {ℝ}_{>0}^M$ the time duration of each piece.

A $m$-dimensional $M$-piece trajectory $p(t)$ is represented by piecewise polynomials. And $i^{th}$ piece $p_i(t)$ is defined as a multi-degree polynomial ($Q=5$ in this paper)

$$p_i(t)={𝐜}_i^T𝜷(t),\forall t\in [0,T_i],$$ (16)

where ${𝐜}_i\in {ℝ}^{(Q+1)\times m}$ is the coefficient matrix and $𝜷(t)={[t^0,t^1,\mathrm{\cdots },t^Q]}^T$ is the natural basis.

For an s-integrator ($s=3$ in this paper) chain dynamics system, a $M$-piece $2s-1$ degree trajectory $p(t)$ is defined by constant boundaries and minimum control effort $\{𝐪,𝐓\}$. Furthermore, MINCO is advanced in convert $\{𝐪,𝐓\}$ to $\{𝐜,𝐓\}$ using a linear-time and space parameter mapping $𝐜=ℳ(𝐪,𝐓)$, where $𝐜={({𝐜}_1^T,\mathrm{\cdots },{𝐜}_M^T)}^T$ is polynomial coefficients.

V-B Problem Formulation

After determining the desired formation shape in Sec.IV, we expect a cluster of trajectories for swarm robots, which are smooth, collision-free, and formation maintained. In practice, navigating swarm robots in an unknown dense environment with FOV-limited sensors and onboard computer requires an efficient real-time planner focusing on local information. Besides, centralized optimization is limited by the scale of the swarm. Therefore, we choose a distributed local trajectory optimization for formation flight as follows

$\underset{𝐪,𝐓}{\min }$ ${\displaystyle {\int }_{t_0}^{t_M}}{\Vert p^{(s)}(t)\Vert }^2𝑑t+\rho \cdot T_{\mathrm{\Sigma }},$ (17a)

$s.t.$ $p(t)={ℳ}_{𝐪,𝐓}(t),\forall t\in [t_0,t_M],$ (17b)

${𝐩}^{[s-1]}(0)={\overline{𝐩}}_0,$ (17c)

${𝐩}^{[s-1]}(t_M)={\overline{𝐩}}_f,$ (17d)

$\mathscr{H}(p(t),\mathrm{\dots },p^{(s)}(t))⪯\mathrm{𝟎},\forall t\in [t_0,t_M].$ (17e)

We define costs (17a) for smoothness and aggressiveness to achieve smooth and efficient flight. $\rho$ is time regularization parameter, $T_{\mathrm{\Sigma }}={\sum }_{i=1}^M{T}_i$. The state of robot $p(t)$ (17b) is parameterized by the optimization variables $\{𝐪,𝐓\}$. ${𝐩}^{[s-1]}(t)={(p{(t)}^T,\dot{p}{(t)}^T,\mathrm{\dots },p^{(s-1)}{(t)}^T)}^T\in {ℝ}^{ms}$ represents the higher-order derivatives of a chain dynamic system with $s$-integrator. Boundary conditions involve initial state ${\overline{𝐩}}_0\in {ℝ}^{ms}$ (17c) and terminal state ${\overline{𝐩}}_f\in {ℝ}^{ms}$ (17d). Continuous-time constraints $\mathscr{H}$ (17e) include swarm formation similarity, dynamic feasibility, obstacle avoidance, and swarm reciprocal avoidance.

V-C Constraints Transcription

To solve the continuous constrained optimization problem (17) in real-time, we use the optimization variable of MINCO (15) to eliminate all kinds of equality constraints (17b)-(17d). And penalty function method [40] is used to deal with the inequality constraints (17e). Then, every integral is evaluated by a finite sum of sample points. Finally, the continuous constrained optimization problem is converted to a discrete unconstrained optimization problem

$\underset{𝐪,𝐓}{\min }{\displaystyle \sum _x}{\lambda }_{\star }{\tilde{J}}_{\star }(𝐪,𝐓,\delta ),$ (18)

where ${\tilde{J}}_{\star }$ are various terms of cost function or penalties, and ${\lambda }_{\star }$ are relative weights. Subscripts $\star =\{f,e,t,o,r,d\}$ ($f$) swarm formation similarity, ($e$) denote control effort, ($t$) total time, ($o$) obstacle avoidance, ($r$) swarm reciprocal avoidance, ($d$) dynamic feasibility. $\delta$ is the sampling time interval.

In our previous work [4], we used the fixed number sampling points ${\widehat{𝐩}}_{i,j}=p_i((j/{\kappa }_i)\cdot T_i)$ to transform the optimization problem, where $p_i(t)$ is the $i^{th}$ piece trajectory and ${\kappa }_i$ is the fixed sample number on this piece. However, considering that the total time $T_{\mathrm{\Sigma }}$ changes during the optimization process, the fixed number sampling points ${\widehat{𝐩}}_{i,j}$ are difficult to space on the whole trajectory equally. Therefore, we take fixed time interval sampling points for the whole trajectory to ensure the accuracy of the penalty function sampling transformation

${\widetilde{𝐩}}_j(t)=p_i(j\delta -{\displaystyle \sum _{l=1}^{i-1}}{T}_l),$ (19)

$j\in \{0,\mathrm{\cdots },\kappa \},\kappa =\lfloor {\displaystyle \frac{T_{\mathrm{\Sigma }}}{\delta }}\rfloor ,$

where $\kappa$ is the sample number and $T_l$ is the preceding time for any .

For the trajectory planning of swarm robots, the fixed time interval sampling points ${\widetilde{𝐩}}_j(t)$ can simplify the optimization problem. Compared with ${\widehat{𝐩}}_{i,j}$, the timestamp corresponding to ${\widetilde{𝐩}}_j(t)$ is fixed, so the states of other robots at this timestamp are also constant during the optimization process. Therefore, it is feasible to calculate the states of other robots w.r.t ${\widetilde{𝐩}}_j(t)$ according to the broadcast trajectories before optimization. Then we can solve the uniform formation position sequence optimization (12) in advance and use ${\widehat{𝐩}}_i^{\ast }$ to replace the formation similarity metric $f_s$ in trajectory optimization (17a) of $i^{th}$ robot. This decoupled formation trajectory optimization results in higher computational efficiency, making our method suitable for large-scale swarm robots.

Despite the optimization problem is not differentiable when sampling number $\kappa$ changes, the cost function remains continuous w.r.t. time duration $𝐓$. In this paper, we use the quasi-Newton method proposed in [37] to solve the non-smooth discrete unconstrained optimization problem (18).

V-D Cost Functions and Gradients

Given the fixed sampling time interval $\delta$, we can evaluate the cost functions and gradients of the whole trajectory by a finite sum of sampling points ${\widetilde{𝐩}}_j(t)$. The cost of various general purpose penalties at $j^{th}$ sampling points is

$${𝒫}_{\star }(𝐜,𝐓,j\delta )={𝒫}_{\star }({\widetilde{𝐩}}_j(t)),$$ (20)

then the cost function ${\tilde{J}}_{\star }$ in (18) is calculated as follows

${\tilde{J}}_{\star }(𝐪,𝐓,\delta )$ $=J_{\star }(𝐜,𝐓,\delta ),$

$=\delta {\displaystyle \sum _{j=0}^{\kappa }}{\overline{\omega }}_j{𝒫}_{\star }(𝐜,𝐓,j\delta )+$ (21)

${\displaystyle \frac{1}{2}}(T_{\mathrm{\Sigma }}-\kappa \delta )\left[{𝒫}_{\star }(𝐜,𝐓,\kappa \delta )+{𝒫}_{\star }(𝐜,𝐓,T_{\mathrm{\Sigma }})\right],$

where $({\overline{\omega }}_0,{\overline{\omega }}_1,\mathrm{\cdots },{\overline{\omega }}_{\kappa -1},{\overline{\omega }}_{\kappa })=(1/2,1,\mathrm{\cdots },1,1/2)$ are the orthogonal coefficients following the trapezoidal rule [41]. And MINCO allows any second-order continuous cost function ${\tilde{J}}_{\star }(𝐪,𝐓)$ to be represented by $J_{\star }(𝐜,𝐓)$. Hence, $\partial {\tilde{J}}_{\star }/\partial 𝐪$ and $\partial {\tilde{J}}_{\star }/\partial 𝐓$ can be efficiently obtained from $\partial J_{\star }/\partial 𝐜$ and $\partial J_{\star }/\partial 𝐓$ respectively, which is benefit to the construction and solution of the optimization problem. In (19), the sampling time $t=j\delta -{\sum }_{l=1}^{i-1}{T}_l$ is related to the preceding time $T_l$, so the gradient of $J_{\star }$ w.r.t ${𝐜}_i$ and $T_l$ are computed as

$$\frac{\partial J_{\star }}{\partial {𝐜}_i}=\frac{\partial J_{\star }}{\partial {𝒫}_{\star }}\frac{\partial {𝒫}_{\star }}{\partial {\widetilde{𝐩}}_j(t)}\frac{\partial {\widetilde{𝐩}}_j(t)}{\partial {𝐜}_i},$$ (22)

$$\frac{\partial J_{\star }}{\partial T_l}=\frac{\partial J_{\star }}{\partial {𝒫}_{\star }}\frac{\partial {𝒫}_{\star }}{\partial {\widetilde{𝐩}}_j(t)}\frac{\partial {\widetilde{𝐩}}_j(t)}{\partial t}\frac{\partial t}{\partial T_l},$$ (23)

(24)

where the calculation of $\partial J_{\star }/\partial {𝒫}_{\star }$ is simple and the details of ${𝒫}_{\star }({\widetilde{𝐩}}_j(t))$ for various general purpose are given as follow.

V-D1 Cost of Swarm Formation Similarity ${𝒫}_f$

In Sec.IV-C, we decouple the formation similarity error metric from trajectory optimization by constructing an unconstrained optimization problem to calculate the uniform optimal formation position sequence ${\widehat{𝐩}}_i^{\ast }$ for each sampling point. This improvement avoids multiple calculations of formation similarity metric $f_s$. Then, we use the quadratic form to calculate the cost of swarm formation similarity

$${𝒫}_f({\widetilde{𝐩}}_j(t))=\max {\{\parallel {\widetilde{𝐩}}_j(t)-{\widehat{𝐩}}_{i,j}^{\ast }{\parallel }^2,0\}}^3.$$ (25)

V-D2 Control Effort $J_e$

The $s^{th}$ ($s=3$ in this paper) control input for the trajectory and its gradients are written as

$$J_e=\sum _{i=1}^M{\int }_0^{T_i}{\Vert p_i^{(s)}(t)\Vert }^2𝑑t,$$ (26)

$$\frac{\partial J_e}{\partial {𝐜}_i}=2\left({\int }_0^{T_i}{𝜷}^{(s)}(t){𝜷}^{(s)}{(t)}^T𝑑t\right){𝐜}_i,$$ (27)

$$\frac{\partial J_e}{\partial T_i}={𝐜}_i^T{𝜷}^{(s)}(T_i){𝜷}^{(s)}{(T_i)}^T{𝐜}_i.$$ (28)

V-D3 Total Time $J_t$

In order to ensure the aggressiveness of the trajectory, we minimize the total time $J_t={\sum }_{i=1}^M{T}_i$. The gradients are given by $\partial J_t/\partial 𝐜=\mathrm{𝟎}$ and $\partial J_t/\partial 𝐓=\mathrm{𝟏}$.

V-D4 Cost of Obstacle Avoidance ${𝒫}_o$

Inspired by[42], obstacle avoidance penalty $J_o$ is computed using Euclidean Signed Distance Field (ESDF). We penalize the sampling points which are too close to the obstacles

$${𝒫}_o({\widetilde{𝐩}}_j(t))=\max {\{{\psi }_o({\widetilde{𝐩}}_j(t)),0\}}^3,$$ (29)

$${\psi }_o({\widetilde{𝐩}}_j(t))=d_o-d_o({\widetilde{𝐩}}_j(t)),$$ (30)

where $d_o$ is the safety threshold set according to the actual situation and $d_o({\widetilde{𝐩}}_j(t))$ is the distance between ${\widetilde{𝐩}}_j(t)$ and the closest obstacle around it. The gradient of ${𝒫}_o$ w.r.t ${\widetilde{𝐩}}_j(t)$ is

$$\frac{\partial {𝒫}_o}{\partial {\widetilde{𝐩}}_j(t)}=-\nabla d^T,$$ (31)

where the $\nabla d$ is the gradient of ESDF in ${\widetilde{𝐩}}_j(t)$.

V-D5 Cost of Swarm Reciprocal Avoidance ${𝒫}_r$

We penalize ${\widetilde{𝐩}}_j(t)$ when it is too close to the trajectories $p_{\varphi }(t),\varphi \in \mathrm{\Phi }$ at the fixed timestamp $t=j\delta$, where $\mathrm{\Phi }$ represents the all other robots in the swarm. Compared to our previous work[4], the state of other robots with fixed timestamp $p_{\varphi }(j\delta )$ are constant during the optimization process and do not produce a gradient w.r.t $𝐓$ for the cost function $J_r$. So the optimization problem and the gradients are simplified.

The cost of swarm reciprocal avoidance is defined as

$${𝒫}_r({\widetilde{𝐩}}_j(t))=\sum _{\mathrm{\Phi }}\max {\{{\psi }_r({\widetilde{𝐩}}_j(t),p_{\varphi }(j\delta )),0\}}^3,$$ (32)

$${\psi }_r({\widetilde{𝐩}}_j(t),p_{\varphi }(j\delta ))=d_r^2-{\Vert {\widetilde{𝐩}}_j(t)-p_{\varphi }(j\delta )\Vert }^2,$$ (33)

where $d_r$ is the safe clearance between each robot. And the gradient of ${𝒫}_r$ w.r.t ${\widetilde{𝐩}}_j(t)$ is

$$\frac{\partial {𝒫}_r}{\partial {\widetilde{𝐩}}_j(t)}=-2{({\widetilde{𝐩}}_j(t)-p_{\varphi }(j\delta ))}^T.$$ (34)

V-D6 Cost of Dynamic feasibility ${𝒫}_d$

We limit the maximum value of velocity and acceleration to guarantee that the robots can execute the trajectory.

${𝒫}_d({\widetilde{𝐩}}_j(t))$ $={𝒫}_{d,v}({\widetilde{𝐩}}_j(t))+{𝒫}_{d,a}({\widetilde{𝐩}}_j(t)),$ (35)

${𝒫}_{d,v}({\widetilde{𝐩}}_j(t))$ $=\max {\{\parallel {\dot{\widetilde{𝐩}}}_j(t){\parallel }^2-v_m^2,0\}}^3,$

${𝒫}_{d,a}({\widetilde{𝐩}}_j(t))$ $=\max {\{\parallel {\ddot{\widetilde{𝐩}}}_j(t){\parallel }^2-a_m^2,0\}}^3,$

where $v_m$ and $a_m$ are the maximum velocity and acceleration.

V-E Discussion on solution quality of trajectory optimization

The proposed trajectory optimization process (17) aims to solve a challenging multi-stage Linear Quadratic Minimum Time (LQMT) problem, which is inherently non-convex and non-linear. Additionally, incorporating ESDF for obstacle avoidance introduces further non-convex constraints. As a result, guaranteeing the global optimal solution with the quasi-Newton method is not always possible. To address concerns regarding local minima and infeasible solutions, we have implemented measures that prioritize safety and dynamic feasibility while maintaining high-performance formation flight.

Firstly, we utilize hybrid-A* searching algorithm[42] to generate initial trajectories that are collision-free and dynamically feasible, ensuring a valid final solution trajectory. During optimization, we give greater weight to obstacle avoidance and dynamic constraints to prioritize safety and feasibility. Additionally, we conduct collision checks on trajectories to enhance safety. Moreover, our distributed swarm optimization framework effectively mitigates the impact of local minima on overall formation performance. Implementing these measures, our method reliably achieves robust formation flight while maintaining computational efficiency.

VI Swarm reorganization method

During the formation navigation, the swarm could encounter many unfavorable conditions, such as highly constrained space, inappropriate assignment of tasks, and sudden formation switching commands. To recover from these situations, we present a swarm reorganization method. The method aims to generate high-quality local goals which satisfy the desired formation distribution and respect the current states of each robot. With these local goals, the swarm can reform the desired shape quickly, even in highly constrained environments such as narrow corridors or holes.

Unlike high-frequency distributed formation trajectory optimization, the swarm reorganization method only runs at a low frequency when the stable state of formation flight is destroyed or about to be destroyed. The method first calculates the formation constraint awareness and then solves an optimal formation ALignment and task ASsignment problem (ALAS). The former awareness distributedly quantifies the conflict between formation maintenance and obstacle avoidance of each robot, while the latter solves the ALAS problem centrally.

VI-A Formation Constraint Awareness

First, we need to derive weights to indicate how severely the environment constrains the robots. These weights are called formation constraint awareness, which should be determined by the current formation status and obstacle information.

Inspired by our previous work [43], we hope to describe the conflict degree based on the relationship between different gradient information. Firstly, we retrieve the ESDF distance $d_o(p)$ and the corresponding obstacle gradient $\nabla d(p)$. Meanwhile, we calculate the current gradient $\nabla f_s(p)$ of the formation similarity term. Secondly, we calculate the cosine $\beta$ of the angle between gradients $\nabla d(p)$ and $\nabla f_s(p)$

$$\beta =\frac{\nabla d(p)\nabla f_s(p)}{\Vert \nabla d(p)\Vert \Vert \nabla f_s(p)\Vert }.$$ (36)

Then we utilize the sigmoid function to map the cosine of the angle to a conflict coefficient $\eta$

$$\eta (\beta )=\frac{1}{1+e^{(\alpha \beta +\gamma )}},$$ (37)

where $\alpha$ regulates how fast this awareness rises as the cosine value $\beta$ increases, $\gamma$ regulates the dead zone and the activation zone of this angle-based awareness. The conflict coefficient $\eta$ is maximum when the directions of $\nabla d(p)$ and $\nabla f_s(p)$ are opposite, which indicates the most conflicting case. And $\eta$ reaches a minimum when the two gradients have the same direction, which means no conflict.

The formation constraint awareness should also consider the influence of the gradient magnitude and the distance of the current closest obstacle. Hence, we design the formation constraint awareness $g_i$ of $i^{th}$ robot as

$$g_i=\lambda \cdot \eta (\beta )\cdot \frac{\Vert \mathrm{▽}{J}_f(p)\Vert }{d_o(p)}.$$ (38)

We apply the calculation to each robot in the swarm and thus get a formation constraint awareness vector $𝐠=\{g_1,\mathrm{\dots },g_N\}$ of the whole swarm. To distinguish the most constrained ones, we use softmax function to amplify the variance of awareness vector $𝐠$ and normalize the vector

$$𝐰=\text{softmax}(𝐠).$$ (39)

$\lambda$ in (38) is used to adjust the variance of elements in $𝐰$. $𝐰$ is the final awareness vector describing the degree of formation-obstacle conflict of the robots in the swarm.

Figure 4: Illustration of formation alignment and task assignment. (a) Before solving ALAS problem, the initial formation goals suffer from a large transition distance to robots and disordered assignments that may lead to deadlock. (b) With robots at the same positions, after solving ALAS, the formation goals enjoy low distance costs and better assignments.

Figure 5: Illustration of the weighted formation alignment problem in a narrow corridor. When a hexagon formation enters the corridor, the original formation distribution(pink points) cannot be maintained anymore. The upper and lower robots(red points) are severely pressed by the obstacles, and hence have the largest constraint awareness. After solving the awareness-weighted formation alignment, the swarm obtains a new desired formation distribution(blue points) that best matches up with the constrained robots.

VI-B Formation Alignment and Task Assignment Optimization

For formation flights in constrained scenes, e.g. in narrow corridors, the unconstrained robots with lower constraint awareness possess larger space to freely coordinate with other robots, since the obstacles don’t hinder the formation requirement. On the contrary, the constrained robots with larger awareness always fall into the conflict between formation maintenance and obstacle avoidance. Hence, refining the positions of unconstrained robots to match up with the constrained ones is more reasonable when generating local goals for formation reorganization. In this work, we use $𝐰$ in Sec.VI-A to weigh the robots when adjusting the formation distributions and place more weights on the constrained robots.

In this section, we only elaborate on the ALAS problem for local goal generation. Afterward, the global-remap-local-replan strategy uses the generated local goals to reorganize the formation, which is detailed in Sec.VII-A.

Let ${\mathrm{𝐥𝐠}}_i={[l{g}_{ix},l{g}_{iy},l{g}_{iz}]}^T,i=1,\mathrm{\dots },N$ represents the current positions of robots. The desired formation shape template is given by $N$ positions ${𝐪}_j={[q_{jx},q_{jy},q_{jz}]}^T,j=1,\mathrm{\dots },N$. Then, the aligned formation positions ${𝐪}_j^{\prime }$ can be written as

$${𝐪}_j^{\prime }=s\cdot {𝐪}_j+𝐝,$$ (40)

where $s\in ℝ$ is the scaling factor, $𝐝\in {ℝ}^3$ denotes the translation factor. In this work, the formation alignment is determined by a scaling factor and a translation factor.

ALAS is composed of formation alignment and task assignment as shown in Fig. 4. The former aims to find the optimal alignment of the desired formation based on a weighted Euclidean distance cost. And the latter solves the optimal assignment that matches the agents with the local goals.

The task assignment problem is formulated as

$$\underset{\sigma }{\min }\sum _i^n{\Vert {\mathrm{𝐥𝐠}}_i-(s^{\ast }\cdot {𝐪}_{\sigma (i)}+{𝐝}^{\ast })\Vert }^2,$$ (41)

where $\sigma \in S_N$ is the assignment map of the formation task and $S_N$ is the symmetric group of all permutations from the set $\{1,\mathrm{\dots },N\}$ to itself. Problem (41) solves the assignment that minimizes the overall transition distance between current robots and the aligned local goals.

The formation alignment problem is formulated as

$$\underset{s,𝐝}{\min }\sum _i^n{w}_i\cdot {\Vert {\mathrm{𝐥𝐠}}_i-(s\cdot {𝐪}_{{\sigma }^{\ast }(i)}+𝐝)\Vert }^2,$$ (42)

where ${\sigma }^{\ast }$ represents the optimal assignment, $w_i$ is the awareness weight from Sec.VI-A. Problem (42) generates a standard formation that best fits into the current robot positions according to a distance cost weighted by the constraint awareness. Fig. 5 illustrates how the alignment adjusts the formation distribution when the swarm is traversing a corridor.

Problems (41) and (42) are coupled. The whole ALAS problem has three decision variables: scaling factor $s$, translation $𝐝$, and assignment $\sigma$. The goal of ALAS is to find an optimal set of decision variables that minimize both (42) and (41). Note that there is no awareness weight $w_i$ multiplied in formulation (41). Because in the assignment optimization, we only care about the total Euclidean distance cost, which is irrelevant to the degree of the constraint of any agent.

However, for formation alignment using only scaling factor $s$ and translation $𝐝$, [25] proves that the corresponding assignment $\sigma$ can be optimized in a decoupled manner, rather than alternating the two optimization phases iteratively. In [25], the optimal assignment solution ${\sigma }^{\ast }$ is shown invariant w.r.t the changes in formation scaling factor $s$ and translation $𝐝$. And the solution of (41) can be directly optimized by solving the following integer programming with new pseudo costs ${\kappa }_{ij}$

$\underset{\sigma =(x_{ij})}{\min }{\displaystyle \sum _{i=1}^n}{\displaystyle \sum _{j=1}^n}{\kappa }_{ij}{x}_{ij},$ (43)

$\text{where}{\kappa }_{ij}=-{\mathrm{𝐥𝐠}}_i^T{𝐪}_j.$

The formulation (LABEL:equ:task_cost) is independent of the scale parameter $s$ and translation $𝐝$. Hence, (LABEL:equ:task_cost) can be first solved prior to the formation alignment phase. Then we determine the best alignment using the optimized assignment ${\sigma }^{\ast }$. The formation alignment problem with awareness weights is still convex and the closed-form solution to (41) is given In Appendix.-A.

Given the solution of ALAS, the position of generated local goal ${\mathrm{𝐥𝐠}}_i^{\prime }$ for the $i^{th}$ robot is calculated by

$${\mathrm{𝐥𝐠}}_i^{\prime }=s^{\ast }\cdot {𝐪}_{{\sigma }^{\ast }(i)}+{𝐝}^{\ast }.$$ (44)

After the ALAS optimization, the distribution of generated local goals is in the desired formation shape, and respects the formation-obstacle conflict of each robot.

VII Swarm agreement method

VII-A Global-remap-local-replan Strategy

Figure 6: Illustration of GRLR strategy. (a) The local-replan strategy generates initial local goals within the planning horizon. Due to the conflict reference trajectories, the swarm robots are expected to deadlock with the current formation behavior. Therefore, the global-remap strategy calls ALAS method to realign the shape and reassign the task of initial local goals. (b) After solving the refined local goals, the global-remap strategy generates conflict-free reference trajectories. In this way, the swarm formation will converge quickly with the new formation behavior.

1:Global reference trajectory ${𝒯}_{ref}$,   Local trajectory ${𝒫}_i$,   Initial local goals ${\mathrm{𝐥𝐠}}_i$,   Global goals ${𝐆}_i$,   Assignments $\sigma$,   Scale $s$,   Translation $𝐝$ ;

2:Initialize: $CallGlobalRemap\leftarrow False$,

3:for each robot i do

4: $g_i\leftarrow$ ConstrainedAwareness(${𝒫}_i$); $▷$ detailed in (38)

5: if $g_i>g_d$ then $▷$ $g_d$ is threshold of $g_i$

6: $CallGlobalRemap\leftarrow True$;

7: $𝐰\leftarrow \text{softmax}(𝐠)$;

8: end if

9:end for

10:if SimilarityError($\cdot$)$>e_{sim,d}$ then $▷$ detailed in (5)

11: $CallGlobalRemap\leftarrow True$;

12: $𝐰\leftarrow \text{softmax}(𝐈)$;

13:end if

14:if $CallGlobalRemap$ then

15: ${\sigma }^{\ast }\leftarrow$ Assignment($\sigma ,s,𝐝$);$▷$ detailed in (41)

16: $s^{\ast },{𝐝}^{\ast }\leftarrow$Alignment($𝐰,{\sigma }^{\ast }$);$▷$ detailed in (42)

17: ${\mathrm{𝐥𝐠}}_i^{\ast }\leftarrow$RemapLocalGoals(${\sigma }^{\ast },s^{\ast },{𝐝}^{\ast }$);

18: ${𝐆}_i^{\ast }\leftarrow$RemapGlobalGoals(${𝐆}_i,{\sigma }^{\ast }$);

19: ${𝒯}_{ref}^{\ast }\leftarrow$GlobalTrajectoryReplan(${\mathrm{𝐥𝐠}}_i^{\ast },{𝐆}_i^{\ast }$);

20: Return ${𝒯}_{ref}^{\ast }$;

21:end if

Algorithm 1 Global-remap strategy

The swarm system cannot quickly recover from disordered states caused by unknown obstacles or sudden changes in the desired formation shape. To address this challenge, we propose a novel approach that utilizes a global-remap-local-remap (GRLR) strategy for trajectory generation. This approach enables us to efficiently navigate complex environments while maintaining the swarm in a “critical state” of the desired formation, effectively balancing coordination and adaptability.

For distributed framework, communication helps the robot to obtain information about others and generate better coordination behavior. Especially in the case of formation transformation, collaborative decision-making can make the swarm formation converge quickly. However, circular dependencies sometimes occur due to communication delays, making it difficult to guarantee the consistency of decisions. Therefore, We design a Global-Remap-Local-Replan (GRLR) strategy for the swarm formation system, which only centrally remaps crucial parameters of formation and maintains the formation coordination through distributed replanning local trajectory.

GRLR strategy comprises the local-replan for a single robot and the global-remap for a formation-level system. The local-replan is a receding horizon incremental planning strategy[44], which allows each robot plans a trajectory within its limited sensing range. The local goals are selected on the global reference trajectories within planning horizon ${\mathrm{\Psi }}_p$, as shown in Fig. 6 (a). The global-remap is an efficient centralized strategy that only remaps the local goals by solving ALAS method and refines the global reference trajectory, as shown in Fig. 6 (b). GRLR strategy is very suitable for distributed asynchronous systems, and there is no deadlock in swarm systems even in the presence of network delays.

The main workflow of the proposed global remap strategy is described in Algorithm.1. Before generating the new formation behavior, the global-remap strategy checks if there are any emergence events (Line 1-12), such as the stable state of the swarm formation being destroyed (Line 9) or about to be destroyed (Line 4). Unlike the local-replan strategy is triggered at a fixed frequency, the global-remap strategy is started by emergent events (Line 13). Then the ALAS method is called to solve the optimal assignment ${\sigma }^{\ast }$ and alignment $s^{\ast },{𝐝}^{\ast }$ (Line 14-15). Global-remap strategy remaps the local goals ${\mathrm{𝐥𝐠}}_i^{\ast }$ and global goals ${𝐆}_i^{\ast }$ and generates a new global trajectory ${𝒯}_{ref}^{\ast }$ for each robot (Line 16-19). Finally, robots form the new formation by executing the local-replan strategy.

In this work, we utilize this semi-distributed GRLR strategy to make the swarm formation adaptable to unknown obstacles or sudden changes in the desired formation shape by replanning local trajectories at 1 Hz and checking emergent events for triggering the global-remap strategy at 20 Hz.

VII-B Formation-level Global Path Finding

Algorithm 2 Formation-level Global Path Finding

1:Tree $𝒯$,  State $𝐳$,  Path cost $c$,  Path $𝐏$;

2:Initialize: ${𝒯}_a\leftarrow \mathrm{\varnothing }\cup \{{\text{z}}_{start}\}$,  ${𝒯}_b\leftarrow \mathrm{\varnothing }\cup \{{\text{z}}_{goal}\}$,   $c_{best}\leftarrow \mathrm{\infty }$,  $FoundSolution\leftarrow False$;

3:for $i=1$ to $N$ do

4: ${\text{z}}_{random}\leftarrow$ Sample(${\text{z}}_{start},{\text{z}}_{goal},c_{best}$);

5: if not $FoundSolution$ then

6: ${\text{z}}_{new}\leftarrow$ GreedyExtendTree(${𝒯}_a,{\text{z}}_{random}$);

7: ${\text{z}}_{conn}\leftarrow$ NearestVertice(${\text{z}}_{new},{𝒯}_b$);

8: $c_{new}\leftarrow$ Connect(${\text{z}}_{new},{\text{z}}_{conn},{𝒯}_a,{𝒯}_b$);

9: if then

10: $c_{best}\leftarrow c_{new}$;

11: $FoundSolution\leftarrow True$;

12: end if

13: else

14: ${\text{z}}_{new}\leftarrow$ ExtendTree(${𝒯}_a,{\text{z}}_{random}$);

15: ${𝒵}_{near}\leftarrow$ NearVertex(${\text{z}}_{new},{𝒯}_a$);

16: ${𝒯}_a\leftarrow$ Rewire(${\text{z}}_{new}$,${𝒵}_{near}$);

17: ${\text{z}}_{conn}\leftarrow$ NearestVertice(${\text{z}}_{new},{𝒯}_b$);

18: $c_{new}\leftarrow$ Connect(${\text{z}}_{new},{\text{z}}_{conn},{𝒯}_a,{𝒯}_b$);

19: if then

20: $c_{best}\leftarrow c_{new}$;

21: end if

22: end if

23: SwapTrees(${𝒯}_a,{𝒯}_b$);

24:end for

25:$𝐏\leftarrow$ RetrievePath(${𝒯}_a,{𝒯}_b$);

26:Return $𝐏$;

We propose a method for formation-level global path finding. Given a start and goal configuration, the planner generates a feasible path connecting them with collision-free intermediate formations. A bidirectional RRT approach is employed to address this path-finding problem.

Many navigation tasks expect the formation to maneuver with a desired scale. In practice, an oversized formation could reduce the vehicle’s communication quality, while an overly small formation scale could increase the risk of inter-vehicle collisions. Unlike the method in [45] which only samples position $𝐩\in {ℝ}^3$ of the formation center, our method adds the formation scale $s$ into the sampling space and makes the formation configuration $𝐳=\{𝐩,s\}\in {ℝ}^3\times {ℝ}^{+}$. In this way, the objective of maintaining desired scale, i.e., minimizing the changes in scale along the path, can be handled by minimizing the $L_2$-norm distance of path in the configuration space $𝐳$.

Navigation in dense environments requires the robots to maintain a formation while letting the obstacles pass through the formation. The method in [46] samples the center position $𝐩$, and then the scale factor $s$ is solved by optimizing the formation placement in obstacle-free convex regions. However, this approach does not allow any obstacle to intersect with the convex hull of the formation and hence wastes many solutions. In contrast, our method directly samples the whole states of the formation configuration $𝐳$ to fully explore the solution space. For each edge of our RRT algorithm, a collision check is conducted on each robot rather than the formation’s whole convex hull to allow obstacles to pass.

The main workflow of our bidirectional RRT planner is described in Algorithm.2, where two trees ${𝒯}_a$ and ${𝒯}_b$ grow towards each other from the initial state ${\text{z}}_{start}$ and the goal state ${\text{z}}_{goal}$ respectively. Before the first solution is found, the bidirectional planner extends the trees in an RRT-Connect[47] manner (Line 5-7). In GreedyExtendTree() and Connect(), the greedy heuristic[47] is adopted to aggressively explore the environment and make tree-connection attempts. After a feasible solution is found, i.e. a finite path cost $c_{new}$ is returned by Connect(), the function Sample() computes an informed sampling set with the new cost $c_{new}$ as depicted in [48]. Then the standard Bidirectional-RRT*[49] procedures are conducted in each loop to update the trees (Line 13-17). Since the path cost is $L_2$-norm distance in the configuration space $𝐳$, informed sampling[48] and Bidirectional-RRT*[49] can guarantee the asymptotic optimality of the path solution.

This formation-level path planner is deployed to render a global path when the global environment information is available. Then global trajectories connecting the waypoints of the global path are generated using MINCO[39], and the framework in Sec.V is employed to optimize the local motions.

VIII Benchmark

In the benchmark, it is important to assess the distortion degree of the current formation ${ℱ}^c$ fairly relative to the desired one ${ℱ}^d$ during flight. Inspired by [22], we solve the following nonlinear optimization problem to find the best similarity transformation ($Sim(3)$ transformation) that aligns ${ℱ}^c$ with ${ℱ}^d$. Then the average formation distance degree ${\overline{e}}_{dist}$ is calculated at the normalized formation scale

$${\overline{e}}_{dist}=\frac{1}{s_o\cdot L}{\int }_{ℒ}\underset{𝐑,𝐭,s}{\min }\sum _{i=1}^n{\Vert {𝐩}_i^d-(s𝐑{𝐩}_i^{𝐜}+𝐭)\Vert }^{2}dl,$$ (45)

where ${𝐩}_i^d$ and ${𝐩}_i^c$ represent the position of $i^{th}$ robot in formation ${ℱ}^d$ and ${ℱ}^c$, respectively. The $Sim(3)$ transformation is composed of a rotation $𝐑\in SO(3)$, a translation $𝐭\in {ℝ}^3$ and a scale expansion $s\in {ℝ}_{+}$. Moreover, $s_o$ is the initial formation scale, and $L$ is the length of formation trajectory $ℒ$. Optimizing the transformation in (45) and applying it to formations, the influence of scaling and rotation is squeezed out so that all the formations can be equitably rated by measuring the position error w.r.t the desired formation. A larger ${\overline{e}}_{dist}$ represents a larger distortion from the desired formation ${ℱ}^d$. Besides, we also calculate the average formation similarity degree ${\overline{e}}_{sim}$

$${\overline{e}}_{sim}=\frac{1}{s_o\cdot L}{\int }_{ℒ}{\Vert \widehat{𝐋}-{\widehat{𝐋}}_{des}\Vert }_F^2𝑑l,$$ (46)

where $\widehat{𝐋}$ and ${\widehat{𝐋}}_{des}$ are detailed in (5) and the formation similarity error ${\Vert \widehat{𝐋}-{\widehat{𝐋}}_{des}\Vert }_F^2$ is proposed in Sec.IV. We show important parameters in Table I used in the following benchmarks, simulations, and real-world experiments. All benchmarks are run on a desktop with an Intel i7-12700 CPU.

TABLE I: Formation parameters of the proposed method

Parameter	Symbol	Value
Similarity error threshold	$e_{sim,d}$	0.05
Constraint awareness threshold	$g_d$	$2/N$
Parameter for regulation in (37)	$\alpha$	5
Parameter for regulation in (37)	$\lambda$	25
Parameter for regulation in (37)	$\gamma$	-1
Sampling time interval ($s$)	$\delta$	0.5
Planing Horizon ($m$)	${\mathrm{\Psi }}_p$	7.5
Max velocity ($m/s$)	$v_m$	1.0
Max acceleration ($m/s^2$)	$a_m$	6.0
Weight for control effort	${\lambda }_e$	10000.0
Weight for total time	${\lambda }_t$	80.0
Weight for swarm reciprocal avoidance	${\lambda }_r$	10000.0
Weight for obstacle avoidance	${\lambda }_o$	10000.0
Weight for swarm formation similarity	${\lambda }_f$	10000.0
Weight for dynamic feasibility	${\lambda }_d$	10000.0

Figure 7: Random forest map and formation types of benchmarks. (a) Random forest map. (b) Regular hexagon shape. (c) Irregular shape. (d) Triangular prism shape. (e) Octahedron shape.

TABLE II: Performance Comparison between Formation Definition Methods

	Formation type	Regular hexagon		Irregular shape		Triangular prism		Octahedron
Scenario		${\overline{e}}_{dist}(\%)$	${\overline{e}}_{sim}(\%)$	${\overline{e}}_{dist}(\%)$	${\overline{e}}_{sim}(\%)$	${\overline{e}}_{dist}(\%)$	${\overline{e}}_{sim}(\%)$	${\overline{e}}_{dist}(\%)$	${\overline{e}}_{sim}$(%)
Same to same	Position	39.120	0.384	35.649	0.384	34.506	0.374	21.534	0.341
	Displacement	16.231	0.172	15.023	0.159	14.952	0.125	11.645	0.153
	Distance	15.489	0.164	14.295	0.131	14.009	0.118	10.285	0.113
	Ours	15.443	0.161	14.287	0.128	14.012	0.119	10.281	0.112
Rotation change	Position	57.456	0.945	51.298	0.732	49.821	0.612	34.124	0.542
	Displacement	39.456	0.412	31.012	0.439	29.546	0.345	23.125	0.353
	Distance	27.513	0.312	22.312	0.234	21.031	0.201	14.173	0.159
	Ours	19.234	0.218	17.032	0.171	15.013	0.151	12.146	0.138
Scale change	Position	59.654	1.098	58.416	0.784	53.246	0.741	37.845	0.555
	Displacement	42.516	0.629	40.021	0.624	39.412	0.398	29.845	0.395
	Distance	59.542	1.030	59.105	0.799	54.126	0.632	38.451	0.578
	Ours	18.332	0.192	18.196	0.185	16.023	0.179	12.264	0.164
Scale & rotation change	Position	62.584	1.304	60.124	0.796	56.213	0.832	41.856	0.635
	Displacement	45.627	0.755	40.194	0.631	40.168	0.423	31.288	0.504
	Distance	62.154	1.250	61.059	0.804	54.317	0.684	42.138	0.684
	Ours	20.231	0.243	18.345	0.204	16.851	0.183	12.357	0.175

TABLE III: Performance Comparison between Formation Navigation Methods

	Formation type	Regular hexagon
Scenario		$successrate(\%)$	$length(m)$	${\overline{e}}_{dist}(\%)$	${\overline{e}}_{sim}(\%)$
Sparse	VRB[26]	75	22.978	57.962	0.984
	Spatial-only	100	21.923	15.023	0.152
	Spatial-temporal	100	21.756	11.240	0.138
Medium	VRB[26]	25	-	-	-
	Spatial-only	100	22.130	14.927	0.158
	Spatial-temporal	100	21.932	13.274	0.153
Dense	VRB[26]	0	-	-	-
	Spatial-only	100	22.283	17.630	0.185
	Spatial-temporal	100	22.133	15.443	0.161

VIII-A Adaptability of Graph-based Formation Definition

To demonstrate the adaptability of graph-based formation definition in Sec.IV, we conduct numerous benchmarks compared to the mainstream formation definition methods concluded in[50], which are categorized based on the controlled variables, namely position-based[51], distance-based[52] and displacement-based methods[53].

We implement these methods in our framework and adapt them to the dense environments by replacing the original cost $J_s$ in (12) to generate uniform optimal formation position sequence ${\widehat{𝐩}}_i^{\ast }$ for each robot $i$. For the position-based method, we set drone_0 as the leader and predefined the absolute relative positions for all other robots to specify the desired formation. So its cost is $J_{s,1}=0$. The distance-based method optimizes the error of desired inter-agent distances

$$J_{s,2}=\sum _{j\in N}{\left(\Vert {𝐩}_i-{𝐩}_j\Vert -\Vert {𝐩}_i^d-{𝐩}_j^d\Vert \right)}^2,$$ (47)

where $N$ is the number of robots, and ${𝐩}_i^d$ is the desired position vector for the $i^{th}$ robot. The displacement-based method optimizes the error of desired relative displacements

$$J_{s,3}=\sum _{j\in N}{\Vert ({𝐩}_i-{𝐩}_j)-({𝐩}_i^d-{𝐩}_j^d)\Vert }^2.$$ (48)

Then we simulate four different geometric formation types in a high-density environment of $40\times 15m$ size with randomly generated obstacles, as shown in Fig. 7 (a). 2D and 3D formations with irregular and regular geometries are considered, namely formation types in a regular hexagon, irregular geometry, triangular prism, and octahedron, as shown in Fig. 7 (b)-(e). To fully compare the adaptability of these methods, we design four different scenarios considering both scaling and rotational variation of formation shape. The formation’s initial and final positions may differ in scale and rotation. Then the scenarios are corresponding categorized as ’Same to same’, ’Rotation change’, ’Scale change’, and ’Scale & rotation change’. We test each method 20 times for each scenario and formation type. The corresponding results over ${\overline{e}}_{dist}$ (45) and ${\overline{e}}_{sim}$ (46) are summarized in Table II.

Unlike our graph-based formation definition, in other methods, changing the scale and rotation of the formation is not permitted during the flight. As shown in Table II, in the same formation type and same scenario, the data states that our method achieves promising results with almost the lowest ${\overline{e}}_{sim}$ and ${\overline{e}}_{dist}$. Moreover, our method shows the lowest error growth rate when the scenario becomes more complicated. In the same scenario, the distortion degrees of all methods decrease with the change of formation type from 2D to 3D centrosymmetric structure, which shows that the formation maintenance is also related to the structural stability of the formation itself. In addition, the distance-based method is invariant to the rotation and achieves relatively acceptable performance in the ’Rotation change’ scenario. Nevertheless, it can not handle size-variant cases. Similarly, other methods are sensitive to rotation or scaling, leading to significant performance degradation in such scenarios. Generally speaking, our graph-based formation definition method achieves scaling and rotational invariance. The invariance improves the formation flight’s adaptability and outperforms the mainstream methods in complicated scenarios.

Figure 8: Comparison of formation flight under improper initial conditions with and without ALAS. (a) In the case without ALAS, the executed trajectories (orange lines) are winding, and the formation shape (red lines) converges slowly due to the crossed global trajectories (blue lines). (b) In the case of ALAS, swarm reorganization makes the formation flight process orderly.

Figure 9: Comparison of formation flight in a constrained environment with and without ALAS. The light blue area represents a wall with a hole in which swarm robots should shrink to pass through. (a) In the case without ALAS, the formation shape (red lines) is deformed to force through the area. (b) In the case with ALAS, the formation shape actively gets smaller to pass the area. (c) Formation flight with ALAS can decrease each robot’s constraint awareness (red dotted lines). This improvement maintains the formation shape with lower similarity error (red lines).

VIII-B Predictability of Spatial-Temporal Trajectory Optimization

To prove the predictability of formation trajectory optimization in Sec.V, we compare our work with the virtual rigid body (VRB) method [26], a SOTA formation control framework that avoids obstacles using potential fields. Moreover, we also compare the performance between the spatial-only and spatial-temporal optimization to illustrate the importance of the time domain for formation flight. We simulate seven drones flying in a regular hexagon from one side to another with a velocity limit of $0.5m/s$. The cluttered area is of $30\times 15m$ size, and three obstacle densities are tested for comparison. Parameters are finely tuned for the best performance of each method.

The results are summarized in Table III, which indicates that the VRB method [26] has an unsatisfactory success rate when dealing with medium and dense obstacles. This is mainly due to the short-term obstacle avoidance generated by multiple interacting potential fields, which often leads to local minima near the corridors, causing robots to become trapped and fail. However, optimization methods consider the future movement of formation, so they can balance the formation maintenance and obstacle avoidance but not break the formation shape. Therefore, optimization methods achieve better performance and maintain the success rate.

We can also conclude that the spatial-temporal method is much more effort-efficient, robust, and flexible when considering temporal optimization. The spatial-only method cannot adjust the trajectory in the time domain, which leads to excessive spatial deformations of the trajectory. So the trajectory length and the formation error ${\overline{e}}_{sim}$ and ${\overline{e}}_{dist}$ are larger in the spatial-only method.

VIII-C Elasticity of Swarm Reorganization Method

To validate the swarm reorganization methods in Sec.VI, we design two benchmarks to illustrate the necessity of task assignment and the adaptability of formation alignment.

Firstly we design a comparison of regular hexagon formation flight under improper initial conditions with and without ALAS to validate the necessity of formation task assignment, as shown in Fig. 8. The blue lines represent each robot’s global trajectory and its assigned tasks in the formation. In Fig. 8(a), the global trajectories are partially crossed due to inappropriate task assignment, leading to trajectory optimization conflicts. So the executed trajectories shown by orange lines look very disordered, and the formation shape shown by red lines converges slowly. In Fig. 8(b), the above problems are effectively resolved by considering ALAS. After one calculation of ALAS, the swarm robots reassign formation tasks and quickly reach a swarm consensus. Then the swarm formation smoothly converges to the desired shape and navigates to the destination in an energy-efficient way. The results of this test validate the necessity of task assignment.

Then, we compare formation flight with and without ALAS when passing through a constrained hole to display the adaptability of ALAS. The results in Fig. 9(a) and Fig. 9(b) show that the formation shape may be deformed when passing through the corridor without ALAS. Otherwise, the case with ALAS can adaptively adjust the formation shape to the constrained environments. From the quantitative analysis results in Fig. 9(c), the case with ALAS can quickly adjust the formation scale and make the swarm reach a consensus so that the formation similarity error and constraint awareness of each robot decline rapidly. However, in the case without ALAS, a higher similarity error and constraint awareness are maintained until the swarm formation leaves the hole, which means the swarm formation is always within the limitations of the environment so that the formation shape cannot converge. This benchmark proves the adaptability of formation alignment.

TABLE IV: Comparison between Global path finding methods

	Alonso-Mora’s method[46]	Ours
Sampling time $(s)$	2.0	2.0
Desired scale $(m)$	3.0	3.0
Path length $(m)$	50.77	24.10
Min scale along the path $(m)$	1.88	2.98

Figure 10: Comparison of Alonso-mora’s method and our formation-level path finding method in the same constrained map. (a) This method samples convex regions (purple polyhedra) in free space and connects them if the intersections are traversable in formation. Because the convex regions must be generated in the safe space, this method is too conservative to generate a longer path with a smaller scale. (b) Our method directly samples the whole states of the 3D-scale formation configuration to fully explore the solution space to allow the formation to pass through tiny obstacles.

Figure 11: Comparison of time consumption between previous work and proposed work [4]. We compare the performance of two methods in the 15-robots scenario. The detailed comparison results are shown in Table V.

VIII-D Resilience of Swarm Agreement Method

To highlight the resilience of our swarm agreement method, we compare Alonso-Mora’s global planning method[46] and our formation-level path finding method in a constrained map, as shown in Fig. 10. This map comprises several blocks and some tiny obstacles, in which swarm robots need to find the path that allows the formation to pass safely. We test the planners 20 times, and Table IV shows the averaged resultant data. The results state that Alonso-Mora’s method[46] is unsatisfactory in dealing with this scenario. Method [46] has no penalty for scale changes of formation, which could choose the corridor route that leads to sudden changes in formation scale. Moreover, it can not handle tiny obstacles and thus yield to a longer path with smaller scales, as shown in Fig. 10(a). Unlike Alonso-Mora’s method, our formation-level path-finding method directly samples in the augmented 3D-scale space and can better maintain the desired formation scale. As shown in Fig. 10(b), our method generates a much shorter path and only sacrifices a small quantity of formation scale. Therefore, our method can handle the map with blocks and tiny obstacles and find safe guidance for swarm formation, which is more suitable for dense environments.

VIII-E Efficiency of Decoupled Formation Optimization

TABLE V: Performance comparison in 15-drones scenario

	Previous method[4]	Proposed method
Time consumption $(ms)$	141.7	38.2
success rate $(\%)$	95.0	100.0
length $(m)$	47.387	45.282
${\overline{e}}_{dist}$ $(\%)$	12.439	11.724
${\overline{e}}_{sim}$ $(\%)$	0.147	0.139

We compare our proposed decoupled formation optimization with the previously coupled formation optimization[4] which directly calls formation similarity distance metric (5) multiple times in the optimization process. To ensure a fair comparison, we exclude the ALAS problem during this benchmark. Both methods’ results are shown in Fig. 11. The previous method[4] only supports small-scale swarm formation since the computation time grows exponentially. Thanks to the decoupled approach, the time consumption of our proposed method for a swarm of 42 robots is not more than 150ms, which can support the real-time application for large-scale swarms.

We select the experimental data from the 15-drones scenario, as presented in Table V. Our method not only achieves significantly shorter computation times than the previous method, but also demonstrates better performance in terms of formation maintenance, success rate, and trajectory length. This validates the effective decoupling of our method, leading to comprehensive performance improvements.

IX Real world experiments and simulation

IX-A Real world Experiments

Our method is integrated with an autonomous distributed aerial swarm system stated in Sec.III. The swarm shares some information, such as trajectories, through a broadcast network, which is the only connection among all robots. As shown in Fig. 12, we use a palm-sized quadrotor platform [5] with local sensors and an onboard computer. Software modules such as estimation, perception, planning, and control are all running onboard in real-time. The maximum number of swarm robots during real-world experiments is 16. Three different real-world experiments are designed to verify the proposed formation flight system’s characteristics fully.

Figure 12: Illustration of palm-sized swarm aerial robots.

Figure 13: Composite snapshots of a regular tetrahedron formation passing through a corridor. The swarm rotates and compresses the formation shape to fly through the narrow space from right to left. The blue line shows the outline of the formation shape.

Figure 14: Illustration of 3-D formation transformation experiment. The blue circle represents quadrotors assigned to the upper position, and the yellow circle represents a lower position. The white line represents the outline of the formation shape. (a) Formation flight starts from unconverged initial positions and improper initial task assignments. So swarm robots call the GRLR strategy to reorganize the formation parameters. (b) After 3 seconds, swarm robots converge to the desired cube shape. Then swarm receives a formation transformation command and quickly calls the GRLR strategy. (c) After 3 seconds, swarm robots converge to the desired double-arrow shape. (e) The executed trajectories (light blue lines) indicate that swarm formation is convergent most time. (d) A more accurate numerical analysis states that similarity error decreases quickly after calling the GRLR strategy.

In the first experiment, as shown in Fig. 13, four quadrotors in a 3-D regular tetrahedron formation manage to pass through a narrow corridor safely. During the flight, the swarm adaptively rotates and compresses the formation shape in response to environmental changes. This test proves that the scaling and rotational invariance provides more flexibility for formation flights in constrained spaces.

Figure 15: Formation multiple transformation simulation. The different colors of the trajectories and robots correspond to the different timestamps. We choose four specific timestamps, and the corresponding formation shapes constitute “FAST” (http://zju-fast.com/).

Figure 16: Formation flight in a maze map. The formation-level global path (blue line) avoids walls while ignoring small obstacles. Then swarm formation follows the global path and generates local trajectories to avoid all obstacles and maintains the formation shape.

Then we design a 3-D formation shape transformation experiment to testify the reorganization ability of our method, as shown in Fig. 14. In Fig. 14(a), the desired formation is cube shape, but the swarm robots navigate from the unconverged initial positions and improper initial task assignments. In the beginning, the swarm robots quickly call the GRLR strategy, then make the formation shape quickly converge to the desired square shape, as shown in Fig. 14(b). Then, the swarm robots receive a formation transformation command from the station laptop and converge to a double-arrow shape, as shown in Fig. 14(c). The top view of the navigation process is shown in Fig. 14(e). From the light blue executed trajectories, we can see that the flight behavior of swarm robots tends to be consistent during time $[4,10]$ and time $[14,20]$. Moreover, during time $[0,4]$ and time $[10,14]$, swarm robots try to reach a swarm consensus and frequently adjust flight behavior to form the desired formation shape. A more accurate numerical analysis can be seen from Fig. 14(d). When the formation system is far from the convergence state or meets a formation transformation command, swarm robots adjust formation alignment and task assignment by calling the GRLR strategy at the time corresponding to the dotted line. According to the similarity error represented by the red line, we can see that except for the non-convergence state at the initial moment and sudden formation transformation, the swarm formation can maintain the desired shape while avoiding obstacles.

Finally, we conduct a 16-drone swarm formation flight experiment outdoors. To the best of our knowledge, this is the largest fully autonomous formation flight experiment in a complex outdoor environment. As shown in Fig. 1, 16 drones form a triangular queue shape and successfully traverse an obstacle-rich area without collision. This area has many street trees, stakes, and street lamps. This experiment proves the robustness and large-scale ability of our proposed method. For more details, please view the experimental video.

Figure 17: Large-scale formation flight from far to near. It can be seen from the executed trajectories that the swarm robots still maintain the formation while avoiding obstacles.

IX-B Simulation Experiments

In order to comprehensively show the characteristics of the proposed formation flight method, we also conduct several simulation experiments to supplement the real-world experiments. All simulation experiments are run on a desktop with an Intel i7-12700 CPU in real-time.

Firstly, we design a formation multiple transformation simulation to testify swarm reorganization ability when coping with emergency changes in dense environments, as shown in Fig. 15. The different colored trajectories represent the time flow of the formation flight from left to right. The robots of the same color correspond to the formation shape at that timestamp. The command of formation transformation is given in real-time instead of being set in advance. The swarm robots must overcome the instantaneous change of desired formation shape and quickly converge to the new formation state. Thanks to the excellent ability of the proposed method to rapidly transform formation in complex environments, we finally generated the acronym ”FAST” for our laboratory.

Then, we set up a special maze simulation environment consisting of walls and many small obstacles such as posts and rings. In this simulation, we aim to verify swarm agreement ability. The formation-level global path finding method first runs in the global map, which is a centralized method proposed in Sec.VII-B. Then, it generates a global path that considers the scale of formation shape, as shown by the blue lines in Fig. 16. The blue cube shapes represent sampling points. After that, 8 robots form a cube formation, navigate following the global path and generate local formation trajectories using distributed methods in Sec.V. This simulation demonstrates that our formation flight system can better accommodate the obstacle constraints and provide safer guidance for the formation flight.

Finally, to test our method’s effectiveness with large-scale irregular formations, we design a double-arrow formation consisting of 30 drones. As depicted in Fig. 17, the swarm successfully avoids the obstacles, and the desired formation is well preserved during the flight.

X Conclusion and future work

In this paper, we analyze the core dilemmas to achieve formation flight in dense environments in detail and accurately summarize PAPER criteria to solve the above problems. Then we propose a hierarchical formation flight architecture composed of graph-based formation definition, distributed formation trajectory optimization, swarm reorganization method, and swarm agreement methods. The proposed complete formation flight system satisfies all PAPER criteria and achieves excellent performance in maintaining cooperative formation flight in dense environments. We design comprehensive benchmarks in terms of adaptability, predictability, elasticity, resilience, and efficiency to verify the outstanding performance of our proposed method. Finally, we conduct abundant real-world experiments and simulations to prove that we solve the problem of autonomous formation flight in dense environments within large-scale swarms.

In the future, we intend to further improve the efficiency of distributed formation flight method through local information propagation of sub-graphs. Furthermore, while our work currently requires an operator to manually determine the formation shape, the research on optimizing formation shape is a promising area. We believe it has the potential to showcase more intelligent and cooperative swarm behavior, ultimately leading to an enhanced task capacity.

-A The Closed-form Solution of Alignment Problem

Now we derive the closed-form solution to the formation alignment problem. Assume after the assignment phase, we have the current robot position ${\mathrm{𝐥𝐠}}_i$ and desired formation ${𝐪}_i$ with optimal matches. We use $c_i$ to denote the $i^{th}$ term of the alignment objective in (42). Then we have

$c_i$ $=w_i\cdot {\Vert {\mathrm{𝐥𝐠}}_i-(s{𝐪}_i+𝐝)\Vert }^2,$ (49)

$=w_i{({\mathrm{𝐥𝐠}}_i-s{𝐪}_i-𝐝)}^T({\mathrm{𝐥𝐠}}_i-s{𝐪}_k-𝐝),$

$=w_i{\mathrm{𝐥𝐠}}_i^T{\mathrm{𝐥𝐠}}_i+s^2{w}_i{𝐪}_i^T{𝐪}_i-2s{w}_i{\mathrm{𝐥𝐠}}_i^T{𝐪}_i+$

$\mathrm{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}2}s{w}_i{𝐪}_i^T𝐝-2w_i{\mathrm{𝐥𝐠}}_i^T𝐝+w_i{𝐝}^T𝐝.$

We use $F$ to denote the alignment objective in (42). Then the objective $F$ can be written as

$F(s,𝐝)$ $={\displaystyle \sum _i^n}{c}_k,$ (50)

$=b_{lg}+s^2{b}_q-2s{b}_{lg,q}$

$+2s{\widehat{𝐪}}^T𝐝-2{\widehat{\mathrm{𝐥𝐠}}}^T𝐝+{𝐝}^T𝐝,$

where

$b_{lg}$ $={\displaystyle \sum _i^n}{w}_i{\mathrm{𝐥𝐠}}_i^T{\mathrm{𝐥𝐠}}_i,$ (51)

$b_q$ $={\displaystyle \sum _i^n}{w}_i{𝐪}_i^T{𝐪}_i,$

$b_{lg,q}$ $={\displaystyle \sum _i^n}{w}_i{\mathrm{𝐥𝐠}}_i^T{𝐪}_i,$

$$\widehat{𝐪}=\sum _i^n{w}_i{𝐪}_i,\widehat{\mathrm{𝐥𝐠}}=\sum _i^n{w}_i{𝐩}_i.$$ (52)

Since the awareness weights $w_i$ are the outputs of a softmax function, we have the property ${\sum }_i^n{w}_i=1$. Note that we use this property to simplify the last term in (50). Now we need to prove that the objective $F$ in (50) is convex w.r.t. scale parameter $s$ and translation $𝐝$. The Hessian matrix of the objective is given by

(53)

The eigenvalues $\lambda$ of the Hessian are as follows

$$\lambda =\{2,\mathrm{\hspace{0.17em}\hspace{0.17em}2},(b_q+1)\pm \sqrt{{(b_q+1)}^2-4b_w}\},$$ (54)

where

$b_w$ $=b_q-{\widehat{𝐪}}^T\widehat{𝐪},$ (55)

$={\displaystyle \sum _i^n}{w}_i{𝐪}_i^T{𝐪}_i-{({\displaystyle \sum _i^n}{w}_i{𝐪}_i)}^T({\displaystyle \sum _i^n}{w}_i{𝐪}_i).$

We need to prove that (55) is non-negative, then all the eigenvalues in (54) will be non-negative, thus the Hessian is semi-positive definite, the objective in (50) will be convex w.r.t. variables $s$ and $𝐝$.

For the non-negative weights $w_i$, we have the property ${\sum }_i{w}_i=1$ from the softmax function, so we construct a convex function $h(𝐱)={𝐱}^T𝐱$, and apply weighted Jensen’s Inequality. We get

$$w_{1}h({𝐪}_1)+\mathrm{\cdots }+w_{n}h({𝐪}_n)\ge h(w_1{𝐪}_1+\mathrm{\cdots }+w_n{𝐪}_n),$$ (56)

$$\Rightarrow \sum _i^n{w}_i{𝐪}_i^T{𝐪}_i\ge {(\sum _i^n{w}_i{𝐪}_i)}^T(\sum _i^n{w}_i{𝐪}_i).$$ (57)

Thus (55) is nonnegative, the objective $F$ is convex. Then we can obtain the closed-form solution of this alignment problem by solving

(58)

The close-form solution of problem (42) is given by

$s^{\ast }$ $={\displaystyle \frac{b_{lg,q}-{\widehat{\mathrm{𝐥𝐠}}}^T\widehat{𝐪}}{b_q-{\widehat{𝐪}}^T\widehat{𝐪}}},$ (59)

${𝐝}^{\ast }$ $=\widehat{\mathrm{𝐥𝐠}}-s^{\ast }\widehat{𝐪}.$

Lun Quan received the B.Eng. degree in control science and engineering in 2019 from Zhejiang University, Hangzhou, China, where he is currently working toward a Ph.D. degree in control science and engineering. His research interests include motion planning for multi-robot systems, swarm intelligence, and autonomous navigation for unmanned vehicles.

Longji Yin received the B.Eng. degree in automation from Zhejiang University, Zhejiang, China, in 2019 and the M.Sc. degree in robotics from Johns Hopkins University, Maryland, U.S. in 2021. In 2022, he worked as a research assistant at Zhejiang University, Zhejiang, China. He is currently working toward a Ph.D. degree in mechanical engineering at the University of Hong Kong, Hong Kong. His research interests include motion planning and mapping for autonomous aerial robots.

Tingrui Zhang received the B.Eng. degree in automotive engineering from Beijing Institute of Technology, Beijing, China, in 2022. He is currently pursuing an M.Phil. degree in control engineering at Zhejiang University, Hangzhou, China. His research interests include motion planning for autonomous robots and numerical optimization.

Mingyang Wang received the B.Eng. degree in control science and engineering from Zhejiang University, Hangzhou, China, in 2022. He is currently working toward an M.Phil. degree in control science and engineering from Zhejiang University, Hangzhou, China. His research interests include motion planning and control for aerobatic flight.

Ruilin Wang received the B.Eng. degree in control science and engineering from Zhejiang University, Hangzhou, China, in 2023. He is currently working toward an M.Phil. degree in artificial intelligence at Sun Yat-sen University, Zhuhai, China. His research interests include aerial robots and motion planning.

Sheng Zhong received the B.Eng. degree in control science and engineering from Zhejiang University, China, in 2023. He is currently working toward a Ph.D. degree in control science and engineering at Hunan University, Changsha, China. His research interests include motion estimation based on event cameras.

Xin Zhou received the B.Eng. degree in electrical engineering and automation from China University of Mining and Technology, Xuzhou, China, in 2019. He is currently working toward the Ph.D. degree in control engineering from Zhejiang University, Hangzhou, China. His research interests include motion planning and mapping for aerial swarm robotics.

Yanjun Cao received his Ph.D. degree in computer and software engineering from the University of Montreal, Polytechnique Montreal, Canada, in 2020. He is currently an associate researcher at the Huzhou Institute of Zhejiang University, as a PI in the Center of Swarm Navigation. He leads the Field Intelligent Robotics Engineering (FIRE) group of the Field Autonomous System and Computing Lab (FAST Lab). His research focuses on key challenges in multi-robot systems, such as collaborative localization, autonomous navigation, perception, and communication.

Chao Xu received his Ph.D. in Mechanical Engineering from Lehigh University in 2010. He is currently Associate Dean and Professor at the College of Control Science and Engineering, Zhejiang University. He serves as the inaugural Dean of ZJU Huzhou Institute, as well as plays the role of the Managing Editor for IET Cyber-Systems & Robotics. His research expertise is Flying Robotics, Control-theoretic Learning. Prof. Xu has published over 100 papers in international journals, including Science Robotics (Cover Paper), Nature Machine Intelligence (Cover Paper), etc. Prof. Xu will join the organization committee of the IROS-2025 in Hangzhou.

Fei Gao received the Ph.D. degree in electronic and computer engineering from the Hong Kong University of Science and Technology, Hong Kong, in 2019. He is currently a tenured associate professor at the Department of Control Science and Engineering, Zhejiang University, where he leads the Flying Autonomous Robotics (FAR) group affiliated with the Field Autonomous System and Computing (FAST) Laboratory. His research interests include aerial robots, autonomous navigation, motion planning, optimization, and localization and mapping.