A Mean-Field Game Control for Large-Scale Swarm Formation Flight in Dense Environments

Abstract

As an important part of cyberphysical systems (CPSs), multiple aerial drone systems are widely used in various scenarios, and research scenarios are becoming increasingly complex. However, planning strategies for the formation flying of aerial swarms in dense environments typically lack the capability of large-scale breakthrough because the amount of communication and computation required for swarm control grows exponentially with scale. To address this deficiency, we present a mean-field game (MFG) control-based method that ensures collision-free trajectory generation for the formation flight of a large-scale swarm. In this paper, two types of differentiable mean-field terms are proposed to quantify the overall similarity distance between large-scale 3-D formations and the potential energy value of dense 3-D obstacles, respectively. We then formulate these two terms into a mean-field game control framework, which minimizes energy cost, formation similarity error, and collision penalty under the dynamical constraints, so as to achieve spatiotemporal planning for the desired trajectory. The classical task of a distributed large-scale aerial swarm system is simulated by numerical examples, and the feasibility and effectiveness of our method are verified and analyzed. The comparison with baseline methods shows the advanced nature of our method.

1. Introduction

As an important part of cyberphysical systems (CPSs), multiple aerial drone systems have shown great advantages and value in performing tasks such as environmental reconnaissance, information transmission & reception, and target tracking in the fields of disaster rescue, geological exploration, and smart cities [1]. Compared with common multi unmanned aerial vehicle (UAV) systems, hundreds or thousands of larger-scale UAV swarms have better maneuverability performance and mission realization efficiency [2]. When the drone swarm performs actual missions, it will fly in formation as needed. In the process of approaching the target, it not only needs to avoid obstacles in time but also needs to consider the trajectory collision avoidance between the drones in the swarm. As a result, the drone swarm requires a large amount of communication and computing resources to guarantee real-time communication information interaction and perception information calculation, which will grow exponentially with the scale.

Research work on the flight control of aerial swarms is extensive. A large number of research studies focus on the formation flight, trajectory planning, and obstacle avoidance of the swarm, and their combination forms are more common. Specifically, Zhang et al. [3] proposed a formation control method based on the communication consensus mechanism and leader–follower strategy, which realized the efficient flight of multiple UAVs in formation around a single static obstacle. In the same year, Zhang et al. [4] proposed an obstacle avoidance control algorithm based on virtual structure and 3-D spatial leader–follower strategy, which realized a triangular formation of multiple UAVs flying through multiple static obstacles with high efficiency. Shao et al. [5] proposed an improved particle swarm optimization algorithm for formation-obstacle avoidance-trajectory planning in a 3-D environmen and completed an experiment in which multiple UAVs formed a formation after avoiding the static cone and ellipsoid obstacles. Combined with an improved artificial potential field method (APF), Wu et al. [6] performed a swarm control for multiple UAVs based on a virtual core structure, solved the local minima problem in APF, and successfully achieved the swarm control of multiple UAVs and the flexible obstacle avoidance of adaptive formation flight in an urban environment. In 2021, the Fei Gao team [7] proposed a distributed swarm trajectory optimization method for formation flight in a dense environment, which realizes multiple UAVs maintain formation flight through unknown obstacle-rich scenarios. Peng et al. [8] put forward a perception sharing and swarm trajectory global optimal algorithm, which experimentally showed that the success rate of obstacle avoidance in a dense obstacle environment could reach at least 80%. In the latest research so far, the formation-obstacle avoidance scenarios of multi-UAV systems have become more and more complex, and the continuously developed optimal control algorithms provide decision-making guarantees for mission flight in the corresponding scenarios. However, the current situation shows that the research schemes in the above literature are all limited by the large-scale of UAVs because the traditional optimization and control technologies deal with the dynamic interaction between individuals, respectively. With the increase of the number of UAVs, the difficulty of solving the optimal strategy will increase significantly.

For the large-scale UAVs flight control problem, some researchers have contributed interesting work. For example, [9] used the mean-field game (MFG) model to solve the real-time flight control problem of large-scale UAVs on the 2-D plane, and the proposed algorithm effectively realizes interaircraft collision avoidance and reduces the energy consumption of UAVs. Shiri et al. [10] combined a federal learning method that can share the parameters of neural network models on UAVs and completed a simulation experiment based on the MFG model for large-scale-UAV swarms flying to the destination under 2-D conditions. Chen et al. [11] proposed a mean-field trust-region policy optimization control method based on multi-agent reinforcement learning, which solved the problem of limited communication range for large-scale UAV flight control. Xu et al. [12] proposed a dual-fields approach, in which the improved APF method overcomes the failure of traditional APF in collision avoidance; the MFG method greatly reduces the communication interference of the UAV group, couples these two methods to adjust the flight trajectory and power consumption, and realizes the flight obstacle avoidance of large-scale UAVs in multiple static obstacles with height differences. Gao et al. [13] proposed an energy-efficient velocity control algorithm for a large number of UAVs based on MFG, which expressed the speed control of large-scale UAVs as a differential game and used the original double mixed gradient method to solve the problem. Their experiments completed that large-scale UAVs bypassed static obstacles in the case of 2-D while minimizing energy consumption. The above documents all use the mean-field method to convert the traditional control 1-vs.-N game into a 1-vs.-1 game, which greatly reduces the interaction frequency of agent systems. Thus, they can deal with large-scale mission scenarios to effectively reduce the communication resources and flight energy consumption. However, it should be pointed out that the state quantities in their numerical solution of the mean-field are all low-dimensional, and the relevant numerical solution problems can be solved using the grid-based method. In terms of solving high-dimensional MFGs, [14] perfectly avoids spatial grids of high-dimensional MFGs based on the machine learning framework, but their work is limited to the deterministic setting ($\sigma$ = 0). Ref. [15] is the first document to solve high-dimensional MFGs under a random setting ($\sigma$ > 0), which exploits the natural connection between MFGs and generative adversarial neural networks (GANs) [16]. Ref. [17] further proposed a numerical solution for solving multi-population high-dimensional stochastic MFGs, a coupled alternating neural network (CA-Net). The numerical experiments completed the coordinated flight of large-scale multipopulation quadcopters to the destination and realize the intergroup interaction and intragroup collision avoidance. The MF method, MFG model, and its numerical solution methods gathered from these documents provide a preliminary solution idea for the flight control problem of large-scale aerial drones in 3-D scenes. However, the current situation shows that these large-scale control research schemes are limited by complex scenarios because most of the state quantities in MF terms are directly used for real-time error measurement. With the complexity of application scenarios, a large number of UAV states may change frequently, resulting in a significant decline in mission success rate.

To summarize, a large-scale agent trajectory planning method that can effectively manage both formation and obstacle avoidance in dense environments is lacking in the literature. On the one hand, obstacle-dense avoidance control for aerial swarm navigation in formation needs to break through the large-scale restrictions; on the other hand, large-scale UAV flight control methods cannot be directly applied to complex scenarios. In practice, a single agent will stay away from obstacles for safety. However, when formation imposes agents tracking targets, it may oppose obstacle avoidance, that is, sometimes formation and obstacle avoidance are contradictory; large-scale itself has the difficulty of communication and calculation, which is exacerbated by frequent formation–neighbor communication and obstacle avoidance environment interactions. How to systematically trade off the conflicting requirements of large-scale, formation, and obstacle avoidance is the key point to accomplishing large-scale noncolliding formation flights.

To bridge the gap, we propose a large-scale agent trajectory optimization method capable of navigating large-scale swarms in formations while avoiding dense obstacles. First, extending the formation method in [7], we model the formation using undirected graphs based on probability densities and define a mean-field-based differentiable metric that assesses the difference between large-scale formation shapes in 3-D workspaces. Our formation similarity metric also retains the translation, rotation, and scale invariance of the formation methods in [7], enabling quantitative evaluation of the overall performance of large-scale formation maintenance. To solve the communication and computing difficulties caused by complex scenarios of large-scale multiagent systems and formulate the trade-offs between formation and obstacle avoidance in large-scale contexts, we design a mean-field game control framework that simultaneously satisfies the dynamic constraints and minimizes energy cost, formation similarity error, and collision penalty. The latter two of them serve as mean-field terms directly related to complex environments, and the functional form we define provides greater flexibility for agent maneuvering. Finally, to verify the effectiveness and practicability of the method, we conduct extensive simulation experiments of distributed aerial large-scale UAVs.

Our main contributions are summarized as:

• A differentiable graph-theory-based mean-field term that quantifies the similarity distance between large-scale three-dimensional formations; a differentiable ellipsoid-based mean-field term that inscribes the potential energy value of dense three-dimensional obstacles.
• A general control framework for complex scenarios of large-scale multiagent systems—mean-field game control, which jointly takes the amount of communication and computation, operating energy consumption, formation similarity, and obstacle avoidance into account.
• A series of simulations with a distributed large-scale aerial swarm system validates the efficiency and robustness of our method. The comparison with baseline methods shows the advanced nature of our method.

2. Two Types of Mean-Field Terms

An aerial swarm formed by a large number of UAVs is considered in this paper, which is traversing an obstacle-rich area in an expected formation. Herein, two types of differentiable mean-field terms for formation similarity metric and obstacle potential energy value engraving are constructed in the following two subsections.

2.1. Formation Mean-Field Term

The large-scale system $\mathcal{N}$ with N agents have a probability distribution of the spatial state $\mathbf{x}$ varying with time t, that is, the mean field $\rho (\mathbf{x},t)$, where $\mathbf{x}=[x_1,x_2,x_3,...,x_n]\in {\mathbb{R}}^n\sim \rho$. Usually the position vector $\mathbf{p}=[x,y,z]=[x_1,x_2,x_3]$. A formation of N agents is modeled by an undirected graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where $\mathcal{V}:=\{1,2,\dots ,N\}$ is the set of vertices, and $\mathcal{E}\subset \mathcal{V}\times \mathcal{V}$ is the set of edges [7]. In graph $\mathcal{G}$, the vertex i represents the $i^{\mathrm{th}}$ agent with position vector ${\mathbf{p}}_i=\left[x_i,y_i,z_i\right]\in {\mathbb{R}}^3$. An edge $e_{ij}\in \mathcal{E}$ that connects vertex $i\in \mathcal{V}$ and vertex $j\in \mathcal{V}$ means the agent i and j can measure the geometric distance between each other. In this work, each agent communicates with all other agents, thus, the formation graph $\mathcal{G}$ is complete. Each edge of the graph $\mathcal{G}$ is associated with a non-negative number as a weight. The weight of edge $e_{ij}$ is given by

$$w_{ij}={∥{\mathbf{p}}_i-{\mathbf{p}}_j∥}^2,(i,j)\in \mathcal{E},$$

(1)

where $\parallel ·\parallel$ denotes the Euclidean norm. Now, the adjacency matrix $\mathbf{A}=\left(A_{ij}\right)\in {\mathbb{R}}^{N\times N}$ and degree matrix $\mathbf{D}=\left(D_{ij}\right)\in {\mathbb{R}}^{N\times N}$ of the formation $\mathcal{G}$ are determined,

$$A_{ij}=w_{ij},$$

(2)

$$D_{ij}=\left\{\begin{array}{cc}\sum _j{A}_{ij}\hfill & \mathrm{if}i=j\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(3)

Thus, the corresponding Laplacian matrix is given by

$$\mathbf{L}=\mathbf{D}-\mathbf{A}.$$

(4)

With the above matrices, the symmetric normalized Laplacian matrix of the graph $\mathcal{G}$ is defined as

$$\widehat{\mathbf{L}}={\mathbf{D}}^{-1/2}{\mathbf{LD}}^{-1/2}=\mathbf{I}-{\mathbf{D}}^{-1/2}{\mathbf{AD}}^{-1/2},$$

(5)

where $\mathbf{I}\in {\mathbb{R}}^{N\times N}$ is the identity matrix.

As a graph representation matrix, Laplacian contains information about the graph structure [18]. To achieve the expected formation for the large-scale swarm, we propose a formation similarity metric as

$${\mathcal{F}}^1\left(\rho (\mathbf{x},t)\right)={\int }_{\Omega }{\int }_{\Omega }{f}^1(\mathbf{x},\widehat{\mathbf{x}})d\rho (\mathbf{x},t)d\rho (\widehat{\mathbf{x}},t),$$

(6)

$$f^1(\mathbf{x},\widehat{\mathbf{x}})={∥\widehat{\mathbf{L}}\left(\mathbf{x}\right)-{\widehat{\mathbf{L}}}_e\left(\widehat{\mathbf{x}}\right)∥}_F^2=tr\left\{{\left(\widehat{\mathbf{L}}-{\widehat{\mathbf{L}}}_e\right)}^T\left(\widehat{\mathbf{L}}-{\widehat{\mathbf{L}}}_e\right)\right\},$$

(7)

where $\Omega$ denotes state space of $\mathbf{x}$, $tr\{·\}$ is the trace of a matrix, $\widehat{\mathbf{L}}$ is the symmetric normalized Laplacian of the current swarm formation, and ${\widehat{\mathbf{L}}}_e$ is the counterpart of the expected formation. Frobenius norm ${\parallel ·\parallel }_F$ is used in our distance metric.

Remark 1.

The function ${\mathcal{F}}^1$ is invariant to the translation, rotation, and scale of the formation because ${\mathcal{F}}^1$ is an integral of $f^1$. $f^1$ is natively invariant to translation and rotation of the formation since the corresponding graph is weighted by the absolute distance between agent positions; scaling invariance is achieved by normalizing graph Laplacian with the degree matrix in (5) [7].

2.2. Obstacle Mean-Field Term

The large-scale system $\mathcal{N}$ with N UAVs needs to traverse an obstacle-rich area. There are K static obstacles. The position of UAV is available. We assume that obstacles can be detected by the UAV vision sensor. In a 3-D urban environment, without loss of generality, obstacles can be considered rectangular solids ${\Omega }_{obs,k}=\left\{(x,y,z)\right||x-x_{0,k}|\le v_{1,k},|y-y_{0,k}|\le v_{2,k},|z-z_{0,k}|\le v_{3,k}\}$, $k=1,...,K$ where $(x_{0,k},y_{0,k},z_{0,k})$ is the center of corresponding obstacle, $(\pm v_{1,k},\pm v_{2,k},\pm v_{3,k})$ are its vertices that are parallel to the $x,y,z$-axis, respectively. For training, we formulate obstacles as

$${\Omega }_{obs,trn,k}:Q_k(x,y,z)=\frac{1}{3v_{1,k}^2}{(x-x_{0,k})}^2+\frac{1}{3v_{2,k}^2}{(y-y_{0,k})}^2+\frac{1}{3v_{3,k}^2}{(z-z_{0,k})}^2\le 1.1k=1,...,K.$$

(8)

We encode ellipsoidal repulsion of obstacles for the generic UAV as

$$f^2\left(\mathbf{x}\right)=\left\{\begin{array}{cc}\sum _{k=1}^K{\gamma }_{obs,k}\frac{1}{Q_k}\hfill & \mathrm{if}\mathbf{x}\in {\Omega }_{obs,trn}\hfill \\ 0\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(9)

where ${\gamma }_{obs,k}$, k=1, ⋯, K, is repulsive force gain coefficient between the generic UAV and other obstacles, ${\Omega }_{obs,trn}={\bigcup }_{k=1}^K{\Omega }_{obs,trn,k}$.

Remark 2.

The cuboid obstacle repulsion is encoded as the unit circumscribed ellipsoid repulsion making this function differentiable [19]. ${\Omega }_{obs,trn}$ contains an obstacle radial bound ten percent more than in circumscribed ellipsoid of ${\Omega }_{obs}$ because we found this additional training buffer alleviates collisions during validation. By training with ellipsoidal repulsion, which has gradient information within the obstacles, we incentivize the model to learn trajectories avoiding the obstacles [20].

To achieve all obstacles avoidance for the large-scale swarm, we propose an obstacle potential energy function as

$${\mathcal{F}}^2\left(\rho (\mathbf{x},t)\right)={\int }_{\Omega }{f}^2\left(\mathbf{x}\right)\rho (\mathbf{x},t)d\mathbf{x}.$$

(10)

Remark 3.

The function ${\mathcal{F}}^2$ can avoid dense obstacles for large-scale swarms, which benefits from the proper construction of 3-D obstacles on the one hand, and the organic combination of obstacle potential and mean-field method on the other hand.

3. A Mean-Field Game Control Framework for Complex Scenarios

In this section, we find the optimal state control for large-scale UAVs which can systematically weigh the energy consumption, formation similarity, and obstacle avoidance of each UAV when it travels complex scenarios. Specifically, we first formulate the optimal control problem for a single UAV. The interaction relationship between drones forms an N-agent noncooperative differential game. As the number of UAVs grows large, the complexity of solving the differential game in the optimal control problem can increase significantly as traditional methods need to separately deal with the increasing interactions. Therefore, we reformulate the large-scale UAVs optimal state control problem in complex scenarios as a mean-field game control problem. Subsequently, we propose a generative-adversarial-network (GAN)-based algorithm to solve the mean-field game.

3.1. Single-UAV Optimal Control Problem

We consider that the continuous-time state $\mathbf{x}$ of the ith UAV has the dynamic–kinematic equation of the following form

$$d{\mathbf{x}}_i=h_i({\mathbf{x}}_i,{\mathbf{x}}_i^{-},{\mathbf{u}}_i)dt+\sigma d{\omega }_i\left(t\right),$$

(11)

where ${\mathbf{x}}_i^{-}$ denotes the states of all other UAVs, ${\mathbf{u}}_i\in {\mathcal{U}}_i$ is the control input (strategy). In the stochastic case, we add a noise term to the dynamics, $\sigma$ is a volatility term—a fixed coefficient matrix, and ${\omega }_i\left(t\right)$ denotes a Wiener process (standard Brownian motion), which is identical and independent among all UAVs. The interpretation here is that we are modeling the situation when the quadcopter suffers from noisy measurements [15]. For example, under wind perturbations, $\sigma$ can be the covariance matrix of the wind velocity [10,13].

Without loss of generality, denoting gravity as g, the acceleration of a UAV with mass m can be written as

$$\left\{\begin{array}{c}\ddot{x_i}=\frac{u_i}{m}(sin\left({\psi }_i\right)sin\left({\phi }_i\right)+cos\left({\psi }_i\right)sin\left({\theta }_i\right)cos\left({\phi }_i\right))\hfill \\ \ddot{y_i}=\frac{u_i}{m}(-cos\left({\psi }_i\right)sin\left({\phi }_i\right)+sin\left({\psi }_i\right)sin\left({\theta }_i\right)cos\left({\phi }_i\right))\hfill \\ \ddot{z_i}=\frac{u_i}{m}cos\left({\theta }_i\right)cos\left({\phi }_i\right)-g\hfill \\ \ddot{{\psi }_i}={{\tau }_{\psi }}_i\hfill \\ \ddot{{\theta }_i}={{\tau }_{\theta }}_i\hfill \\ \ddot{{\phi }_i}={{\tau }_{\phi }}_i\hfill \end{array}\right.,$$

(12)

where $(x,y,z)$ is the spatial position of the UAV, $(\psi ,\theta ,\phi )$ is the angular orientation with corresponding torques ${\tau }_{\psi },{\tau }_{\theta },{\tau }_{\phi }$, and u is the main thrust directed out of the bottom of the aircraft [21]. To fit a control framework, the above second-order system can be turned into a first-order system

$$\dot{{\mathbf{x}}_i}=h_i({\mathbf{x}}_i,{\mathbf{x}}_i^{-},{\mathbf{u}}_i)⟹\left\{\begin{array}{c}\dot{x_i}={v_x}_i\hfill \\ \dot{y_i}={v_y}_i\hfill \\ \dot{z_i}={v_z}_i\hfill \\ \dot{{\psi }_i}={v_{\psi }}_i\hfill \\ \dot{{\theta }_i}={v_{\theta }}_i\hfill \\ \dot{{\phi }_i}={v_{\phi }}_i\hfill \\ {\dot{v_x}}_i=\frac{u_i}{m}(sin\left({\psi }_i\right)sin\left({\phi }_i\right)+cos\left({\psi }_i\right)sin\left({\theta }_i\right)cos\left({\phi }_i\right))\hfill \\ {\dot{v_y}}_i=\frac{u_i}{m}(-cos\left({\psi }_i\right)sin\left({\phi }_i\right)+sin\left({\psi }_i\right)sin\left({\theta }_i\right)cos\left({\phi }_i\right))\hfill \\ {\dot{v_z}}_i=\frac{u_i}{m}cos\left({\theta }_i\right)cos\left({\phi }_i\right)-g\hfill \\ {\dot{v_{\psi }}}_i={{\tau }_{\psi }}_i\hfill \\ {\dot{v_{\theta }}}_i={{\tau }_{\theta }}_i\hfill \\ {\dot{v_{\phi }}}_i={{\tau }_{\phi }}_i\hfill \end{array}\right.,$$

(13)

where $\mathbf{x}={[x,y,z,\psi ,\theta ,\phi ,v_x,v_y,v_z,v_{\psi },v_{\theta },v_{\phi }]}^T\in {\mathbb{R}}^{12}$ is the state with velocities $\mathbf{v}$, and $\mathbf{u}={[u,{\tau }_{\psi },{\tau }_{\theta },{\tau }_{\phi }]}^T\in {\mathbb{R}}^4$ is the control.

For the ith UAV, we consider the cost function of the following form

$${\mathcal{J}}_i\left({\mathbf{x}}_i,{\mathbf{x}}_i^{-},{\mathbf{u}}_i\right)={\int }_0^T{L}_i\left({\mathbf{x}}_i,{\mathbf{u}}_i\right)+{\mathcal{F}}_i\left({\mathbf{x}}_i,{\mathbf{x}}_i^{-}\right)dt+{\mathcal{G}}_i\left({\mathbf{x}}_i^T\right),$$

(14)

where $L_i$ is a running cost incurred by the ith UAV based solely on its actions, ${\mathcal{F}}_i$ is a running cost incurred by the ith UAV based on its interactions with the rest of the swarm, and ${\mathcal{G}}_i$ is a terminal cost incurred by the ith UAV based on its final state.

The running cost $L_i$ is given by

$$L_i\left({\mathbf{x}}_i,{\mathbf{u}}_i\right)=c_1{\parallel {\mathbf{u}}_i({\mathbf{x}}_i,t)\parallel }^2,$$

(15)

where $c_1$ is a constant. The running cost $L_i$ denotes the control effort implemented by the ith UAV. The terminal cost ${\mathcal{G}}_i$ is given by

$${\mathcal{G}}_i\left({\mathbf{x}}_i^T\right)=c_2{\parallel {\mathbf{x}}_i^T-{\mathbf{x}}_0\parallel }_2,$$

(16)

where $c_2$ is a constant, ${\mathbf{x}}_i^T$ is the final state of the ith UAV, and ${\mathbf{x}}_0$ is the target state on which we want the UAVs to reach the objective. The terminal cost ${\mathcal{G}}_i$ denotes the distance between UAV i’s target state and the desired state. The interaction cost ${\mathcal{F}}_i$ is given by

$${\mathcal{F}}_i\left({\mathbf{x}}_i,{\mathbf{x}}_i^{-}\right)=l_1{\mathcal{F}}_i^1\left({\mathbf{x}}_i,{\mathbf{x}}_i^{-},t\right)+l_2{\mathcal{F}}_i^2\left({\mathbf{x}}_i,{\mathbf{x}}_i^{-},t\right)+l_3{\mathcal{F}}_i^3\left({\mathbf{x}}_i,{\mathbf{x}}_i^{-},t\right),$$

(17)

where $l_1$, $l_2$, $l_3$ are constants, ${\mathcal{F}}_i^1\left({\mathbf{x}}_i,{\mathbf{x}}_i^{-},t\right)$ is UAV i’s formation cost,

$${\mathcal{F}}_i^1\left({\mathbf{x}}_i,{\mathbf{x}}_i^{-},t\right)=\frac{1}{N}\sum _{i=1}^N\frac{1}{N}\sum _{j=1}^N{∥\widehat{\mathbf{L}}({\mathbf{x}}_i,{\mathbf{x}}_i^{-},t)-{\widehat{\mathbf{L}}}_e({\widehat{\mathbf{x}}}_j,{\widehat{\mathbf{x}}}_j^{-},t)∥}_F^2,$$

(18)

${\mathcal{F}}_i^2\left({\mathbf{x}}_i,{\mathbf{x}}_i^{-},t\right)$ is UAV i’s obstacle avoidance cost,

$${\mathcal{F}}_i^2\left({\mathbf{x}}_i,{\mathbf{x}}_i^{-},t\right)=\frac{1}{N}\sum _{i=1}^N\left\{\begin{array}{cc}\sum _{k=1}^K{\gamma }_{obs,k}\frac{1}{Q_k\left({\mathbf{x}}_i,t\right)}\hfill & \mathrm{if}{\mathbf{x}}_i\in {\Omega }_{obs,trn}\hfill \\ 0\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(19)

${\mathcal{F}}_i^3\left({\mathbf{x}}_i,{\mathbf{x}}_i^{-},t\right)$ is UAV i’s collision avoidance cost,

$${\mathcal{F}}_i^3\left({\mathbf{x}}_i,{\mathbf{x}}_i^{-},t\right)=\frac{1}{N}\sum _{j\ne i}^N{\mathbb{1}}_{\parallel E^{\frac{1}{2}}\left({\mathbf{x}}_i\left(t\right)-{\mathbf{x}}_i^{-}\left(t\right)\right)\parallel \le e_0},$$

(20)

$${\mathbb{1}}_{\mathbb{A}}\left(\xi \right)=\left\{\begin{array}{cc}1\hfill & \mathrm{if}\xi \in \mathbb{A}\hfill \\ 0\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

where $\mathbb{1}$ means the indicator function, $e_0$ is the safe distance between UAVs, $E=diag(1,1,1/c)$ transforms Euclidean distance into ellipsoidal distance. The interaction cost ${\mathcal{F}}_i$ denotes the sum of formation maintenance and trajectory collision loss about UAV i.

In summary, the optimal control problem faced by UAV i is given by

$$\begin{array}{cc}\hfill & \underset{{\mathbf{u}}_i\in {\mathcal{U}}_i}{inf}{\mathcal{J}}_i\left({\mathbf{x}}_i,{\mathbf{x}}_i^{-},{\mathbf{u}}_i\right)\hfill \\ \mathrm{s}.\mathrm{t}.\hfill & d{\mathbf{x}}_i=h_i({\mathbf{x}}_i,{\mathbf{x}}_i^{-},{\mathbf{u}}_i)dt+\sigma d{\omega }_i\left(t\right).\hfill \end{array}$$

(21)

To this end, the cooperative control of a large number of UAVs in complex environments is formulated as a noncooperative differential game. Therefore, this is an N-player noncooperative game whose well-known solution is the Nash equilibrium (NE), i.e., the control decisions under which no UAV can unilaterally decrease its cost [22].

3.2. Mean-Field Game Control Formulation

As the UAV-swarm number N increases, the complexity of solving the differential game in (21) will increase significantly. Mean-field games (MFGs) are a class of problems that encode large populations of interacting agents into systems of coupled partial differential equations, which overcomes the difficulties of large-scale and information sharing. For the current UAV, the neighborhood interaction of formation and collision avoidance, and the environment interaction of obstacle avoidance will cause the UAV to need a great amount of communication resources and computing resources. Traditional optimization and control methods need to separately deal with the increasing interactions, which leads to the problem of dimensionality explosion. Therefore, we reformulate the large-scale UAVs optimal state control problem in complex scenarios as a mean-field game control (MFGC) problem. Under the framework of MFGC, a generic agent only react to the collective behaviors (mean-field) of all agents instead of the behavior of each agent, which greatly reduces the amount of communication and computation. Now, we can drop the index i since agents are indistinguishable in MFG. Let $\rho (\mathbf{x},t)$ denote the probability of state $\mathbf{x}$ at time t; then, the cost functional (14) is now transformed into

$$\mathcal{J}\left(\mathbf{x},\rho ,\mathbf{u}\right)={\int }_0^T\left\{{\int }_{\Omega }L\left(\mathbf{x},\mathbf{u}\right)\rho (\mathbf{x},t)d\mathbf{x}+\mathcal{F}\left(\mathbf{x},\rho (\mathbf{x},t)\right)\right\}dt+\mathcal{G}\left(\rho (·,T)\right).$$

(22)

Meanwhile, the state dynamics in (11) are transformed into

$$d\mathbf{x}=h\left(\mathbf{x}\right(t),\rho (\mathbf{x},t),\mathbf{u}(\mathbf{x},t\left)\right)dt+\sigma d\omega \left(t\right).$$

(23)

According to Ito’s lemma [23], (23) can be represented in terms of the mean field $\rho (\mathbf{x},t)$ and then will be equivalent to the Fokker–Planck equation given by

$${\partial }_t\rho -\frac{{\sigma }^2}{2}\Delta \rho +\nabla ·\left(\rho h\right)=0.$$

(24)

with cost functional (22) and the Fokker–Planck equation in (24), the high-dimensional MFG which describes the complex-scenario cooperative control of a large number of UAVs is summarized as

$$\begin{array}{cc}\hfill & \underset{\rho ,\mathbf{u}}{inf}\mathcal{J}\left(\mathbf{x}\left(t\right),\rho (\mathbf{x},t),\mathbf{u}(\mathbf{x},t)\right)\hfill \\ \mathrm{s}.\mathrm{t}.\hfill & {\partial }_t\rho -\frac{{\sigma }^2}{2}\Delta \rho +\nabla ·\left(\rho h\right)=0\hfill \\ \hfill & \rho (\mathbf{x},0)={\rho }_0\left(\mathbf{x}\right),\hfill \end{array}$$

(25)

where ${\rho }_0\left(\mathbf{x}\right)$ is the initial probability distribution of the UAV swarm’s states. To this end, the state control for a large number of UAVs is formulated as a high-dimensional MFG. The agents forecast a distribution of the population ${\left\{\rho (·,t)\right\}}_{t=0}^T$ and aim at minimizing their cost. Therefore, at a Nash equilibrium, we have that, for every $\mathbf{x}\in {\mathbb{R}}^n$,

$$\mathcal{J}\left(\mathbf{x},\rho ,\widehat{\mathbf{u}}\right)\le \mathcal{J}\left(\mathbf{x},\rho ,\mathbf{u}\right),\forall \mathbf{u}\in \mathcal{U},$$

(26)

where $\widehat{\mathbf{u}}$ is the equilibrium strategy of an agent at state $\mathbf{x}$. Here, we assume that agents are small, and their unilateral actions do not alter the density $\rho$ [14]. With finite UAVs, it yields an MF approximation that achieves the $ϵ$-NE [10,22].

3.3. GAN-Based Algorithm for Complex-Scenario MFGs

In this section, we propose a generative-adversarial-network (GAN)-based approach to solving the high-dimensional MFG in (25). Inspired by Wasserstein GANs [24], APAC-Net [15,25], we use the variational primal–dual structure of MFG and phrase (25) as a convex–concave saddle-point problem. Then, we propose a GAN-based algorithm to solve the complex-scenario MFG in (25).

Now, we show the underlying primal–dual structure of the complex-scenario MFG and derive the convex–concave saddle-point problem equivalent to (25). Denote $\varphi$ as the Lagrange multiplier, we can put the differential constraint in (24) (Fokker-Planck equation) into the cost function in (22) to obtain the following extended cost function

$$\underset{\varphi }{sup}\underset{\rho ,\mathbf{u}}{inf}{\int }_0^T\left\{{\int }_{\Omega }\rho (\mathbf{x},t)L\left(\mathbf{x},\mathbf{u}(\mathbf{x},t)\right)d\mathbf{x}+\mathcal{F}\left(\rho (·,t)\right)\right\}dt+\mathcal{G}\left(\rho (·,T)\right)-{\int }_0^T{\int }_{\Omega }\varphi (\mathbf{x},t)\left({\partial }_t\rho -\frac{{\sigma }^2}{2}\Delta \rho +\nabla ·\left(\rho (\mathbf{x},t)h(\mathbf{x},t)\right)\right)d\mathbf{x}dt.$$

(27)

Minimizing concerning $\mathbf{u}$ to get the Hamiltonian via $H(\mathbf{x},p)={inf}_u\{ph+L(\mathbf{x},\mathbf{u})\}$ and integrating by parts, we can rewrite (27) as

$$\underset{\rho }{inf}\underset{\varphi }{sup}{\int }_0^T\left\{{\int }_{\Omega }\left({\partial }_t\varphi +\frac{{\sigma }^2}{2}\Delta \varphi +H(\mathbf{x},\nabla \varphi )\right)\rho (\mathbf{x},t)d\mathbf{x}+\mathcal{F}\left(\rho (·,t)\right)\right\}dt+{\int }_{\Omega }\varphi (\mathbf{x},0){\rho }_0\left(\mathbf{x}\right)d\mathbf{x}-{\int }_{\Omega }\varphi (\mathbf{x},T)\rho (\mathbf{x},T)d\mathbf{x}+\mathcal{G}\left(\rho (·,T)\right).$$

(28)

The Formular (28) is the cornerstone of our algorithm.

We solve (28) by training a GAN-like neural network. Although the idea is inspired by Wasserstein GANs [24] and APAC-Net [15,25], the loss function of our algorithm considers the formation similarity ${\mathcal{F}}^1$ and dense obstacle avoidance ${\mathcal{F}}^2$, which provides greater flexibility for large-scale agents maneuvering; details are shown in Algorithm 1. We show the structure and training process of our GAN-based neural network in Figure 1.

Algorithm 1 GAN-based algorithm for complex-scenario MFGs

Require:$\sigma$ diffusion parameter, H Hamiltonian, g terminal cost, f interaction term. Require: Initialize neural networks $N_{\omega }$ and $N_{\theta }$, batch size B. Require: Set ${\varphi }_{\omega }$ and $G_{\theta }$ as in (29). while not converged do          train ${\varphi }_{\omega }$:          Sample batch ${\left\{\left({\mathbf{z}}_b,t_b\right)\right\}}_{b=1}^B$ where ${\mathbf{z}}_b\sim {\rho }_0$ and $t_b\sim \mathrm{Unif}(0,T)$          Obtaining generated data ${\left\{\left({\mathbf{x}}_b,t_b\right)\right\}}_{b=1}^B,{\mathbf{x}}_b\leftarrow G_{\theta }\left({\mathbf{z}}_b,t_b\right)$.          Update discriminator parameters $\omega$ to minimize ${\ell }_{\mathrm{total}}={\ell }_0+{\ell }_t+{\ell }_{\mathrm{HJB}}$                   $\omega \leftarrow \omega -{\eta }_1\nabla {\ell }_{\mathrm{total}}\left(\omega \right)$          train $G_{\theta }:$          Sample batch ${\left\{\left({\mathbf{z}}_b,t_b\right)\right\}}_{b=1}^B$ where ${\mathbf{z}}_b\sim {\rho }_0$ and $t_b\sim \mathrm{Unif}(0,T)$          Update generator parameters $\theta$ to minimize ${\zeta }_t$                   $\theta \leftarrow \theta -{\eta }_2\nabla {\zeta }_t\left(\theta \right)$ end while

First, we initialize the neural networks $N_{\omega }(\mathbf{x},t)$ and $N_{\theta }(\mathbf{z},t)$, then set

$${\varphi }_{\omega }(\mathbf{x},t)=(1-t)N_{\omega }(\mathbf{x},t)+tg\left(\mathbf{x}\right),G_{\theta }(\mathbf{z},t)=(1-t)\mathbf{z}+t{N}_{\theta }(\mathbf{z},t),$$

(29)

where $\mathbf{z}\sim {\rho }_0$ are samples drawn from the initial distribution. The formulation of ${\varphi }_{\omega }$ and $G_{\theta }$ in (29) automatically satisfy the terminal and initial condition, respectively. Moreover, our algorithm encodes the underlying structure of complex-scenario MFGs via (28) and (29), which absolves the neural networks learning for the entire solution of the MFG from scratch.

Our approach for training neural networks includes parallel alternately training $G_{\theta }$ (the state density distribution) and ${\varphi }_{\omega }$ (the value function). In order to obtain the equilibrium of the complex-scenario MFG, we train ${\varphi }_{\omega }$ by first sampling a batch ${\left\{{\mathbf{z}}_b\right\}}_{b=1}^B$ from ${\rho }_0$, and ${\left\{t_b\right\}}_{b=1}^B$ uniformly from $[0,T]$. Then, we compute the push-forward states ${\mathbf{x}}_b=G_{\theta }\left({\mathbf{z}}_b,t_b\right)$ for $b=1,\cdots ,B$. The main loss item for training the discriminator ${\varphi }_{\omega }$ is

$${\mathrm{loss}}_{\varphi }=\underset{{\ell }_0}{\underbrace{\frac{1}{B}\sum _{b=1}^B{\varphi }_{\omega }\left({\mathbf{x}}_b,0\right)}}+\underset{{\ell }_t}{\underbrace{\frac{1}{B}\sum _{b=1}^B{\partial }_t{\varphi }_{\omega }\left({\mathbf{x}}_b,t_b\right)+\frac{{\sigma }^2}{2}\Delta {\varphi }_{\omega }\left({\mathbf{x}}_b,t_b\right)+H\left({\nabla }_{\mathbf{x}}{\varphi }_{\omega }\left({\mathbf{x}}_b,t_b\right),{\mathbf{x}}_b\right)}},$$

(30)

where we can optionally add a regularization term

$${\ell }_{\mathrm{HJB}}=\lambda \frac{1}{B}\sum _{b=1}^B\parallel {\partial }_t{\varphi }_{\omega }\left({\mathbf{x}}_b,t_b\right)+\frac{{\sigma }^2}{2}\Delta {\varphi }_{\omega }\left({\mathbf{x}}_b,t_b\right)+H\left({\nabla }_{\mathbf{x}}{\varphi }_{\omega }\left({\mathbf{x}}_b,t_b\right),{\mathbf{x}}_b\right)+f\left({\mathbf{x}}_b,t_b\right)\parallel$$

(31)

to penalize deviations from the HJB equations. This extra regularization term has also been found effective in, e.g., Wasserstein GANs. Finally, we backpropagate the total loss to update the weights of the discriminator ${\varphi }_{\omega }$. To train the generator, we again sample ${\left\{{\mathbf{z}}_b\right\}}_{b=1}^B$ and ${\left\{t_b\right\}}_{b=1}^B$ as before. The loss of the generator is given by

$${\mathrm{loss}}_G=\underset{{\zeta }_t}{\underbrace{\frac{1}{B}\sum _{b=1}^B{\partial }_t{\varphi }_{\omega }\left(G_{\theta }\left({\mathbf{z}}_b,t_b\right),t_b\right)+\frac{{\sigma }^2}{2}\Delta {\varphi }_{\omega }\left(G_{\theta }\left({\mathbf{z}}_b,t_b\right),t_b\right)+H\left({\nabla }_{\mathbf{x}}{\varphi }_{\omega }\left(G_{\theta }\left({\mathbf{z}}_b,t_b\right),t_b\right),{\mathbf{x}}_b\right)+f\left(G_{\theta }\left({\mathbf{z}}_b,t_b\right),t_b\right)}}.$$

(32)

At last, we backpropagate this loss to update the weights of $G_{\theta }$. In conclusion, in each time slot $t\in [0,T]$, the neural networks will be trained. The generator $G_{\theta }$ will generate the state distribution at time t and the discriminator ${\varphi }_{\omega }$ will get the result of the value function at time t.

4. Simulation Results

In this section, to verify the effectiveness and robustness of our method, we conduct a series of simulation experiments of distributed aerial large-scale UAVs. To demonstrate the competitiveness of our method, we compare it against baseline methods regarding their performances at last.

4.1. Simulation Parameters

For the model hyperparameters, we set $c_1=0.5$ (the weight of the control effort in (15)), $c_2=5$ (the weight of the terminal cost in (16)), and $l_1=5$, $l_2=10$, $l_3=5$ (the weights of the interaction cost in (17)). For the neural networks, there are three linear hidden layers with 100 hidden units per layer in both networks. Residual neural networks (ResNet) are used for both networks, with a skip connection weight of $0.5$. Tanh activation function is used in ${\varphi }_{\omega }$ while ReLU activation function is used in $G_{\theta }$. We choose ADAM as the optimizer with $\beta =(0.5,0.9)$, ${\eta }_1=4\times {10}^{-4}$ (the learning rate for ${\varphi }_{\omega }$), ${\eta }_2=1\times {10}^{-4}$ (the learning rate for $G_{\theta }$), weight decay of ${10}^{-4}$ for both networks, batch size 25, and $\lambda =1$ (the weight of the HJB regularization term in (31)).

4.2. Formation Performance of Large-Scale UAVs

To validate the feasibility and robustness of our method for large-scale UAVs formation, we simulated twenty-five drones flying in desired formations from the bottom left side to the top right side in the 3-D scene, as shown in Figure 2. We experimented with two forms of formations—spiral formation and wave formation, which could confirm the generality of our approach. In Figure 2a, we set the initial points of the drones as $(2cos\left(\frac{2\pi i}{25}\right),2sin\left(\frac{2\pi i}{25}\right)+18,\frac{5i}{24}-10)$, $i=0,1,...,24$; the target points of the drones are $(2cos\left(\frac{2\pi i}{25}\right),2sin\left(\frac{2\pi i}{25}\right)-18,\frac{5i}{24}+5)$, $i=0,1,...,24$. In Figure 2b, the initial points of the drones are $(\frac{20i}{25}-10,sin(\frac{20i}{25}-10)+18,\frac{5i}{24}-10))$, $i=0,1,...,24$; the target points of the drones are $(\frac{20i}{25}-10,sin(\frac{20i}{25}-10)-18,\frac{5i}{24}+5))$, $i=0,1,...,24$.

The test results are shown in Figure 3, which demonstrates that our method can maintain the formation of large-scale swarms in 3-D environments. To certify the generality of our approach in large-scale formations, we design two forms of formations consisting of twenty-five quadrotors—spiral and wave formations. As depicted in Figure 3a,b, the $xoy$ -plane projections illustrate the shape contour of the formation maintenance, which means that the swarm maintains the desired formation well during the flight.

The convergence of our proposed MFG control algorithm and the effectiveness and stability of formation control under different forms of formations are shown in Figure 4. In Figure 4a, we plot the HJB residual errors, i.e., ${\ell }_{\mathrm{HJB}}$ in Algorithm 1, which shows the convergence of MFGC. Without an efficient strategy control, the HJB residuals under different formations are relatively high. HJB residuals drop fast after we apply a series of controls. After around $5\times {10}^4$ iterations, the error curves tend to be stable when we obtain an optimal control for UAVs. In Figure 4b, we plot the formation loss. Without a proper design of the formation cost term, there is no effective formation control, and the formation loss under different formations is relatively high. Formation loss drops fast after we apply our formation cost term. After around $1\times {10}^5$ iterations, the error curves tend to be stable when we obtain an optimal formation control for UAVs.

4.3. Effect of Volatility Term $\sigma$ on Formation

We investigate the effect of the volatility term $\sigma$ on the behavior of the MFGC solutions—the executed trajectories for the formation flight of large-scale UAVs. In Figure 5, we show the solutions for MFGCs using $\sigma =0$, $\sqrt{2}\times {10}^{-1}$ and $\sqrt{10}\times {10}^{-1}$. As $\sigma$ increases, the density of UAVs widens along the paths, and the desired formation maintenance of UAVs decreases due to the added diffusion term in the MFG control equations in (25). From

Table 1. Comparison with interaction costs.

Volatility Parameter $\mathit{\sigma }$	Formation Cost ${\mathcal{F}}^1$	Collision Avoidance Cost ${\mathcal{F}}^3$
0	$1.41\times {10}^{-3}$	$6.35\times {10}^{-2}$
$\sqrt{2}\times {10}^{-1}$	$6.34\times {10}^{-2}$	$6.34\times {10}^{-2}$
$\sqrt{10}\times {10}^{-1}$	$4.30\times {10}^{-1}$	$6.14\times {10}^{-2}$

, it can be seen that larger volatility $\sigma$ corresponds to greater formation cost ${\mathcal{F}}^1$, but the collision avoidance cost ${\mathcal{F}}^3$ is almost invariant. These results are consistent with those in Figure 5.

4.4. Performance of Dense Obstacle Avoidance for Large-Scale UAVs Formation Flying

To verify the effectiveness of our method with large-scale UAVs formation flying through an obstacle-rich area, we design a spiral formation consisting of twenty-five quadrotors. A packed area of nine cuboid obstacles is set up in the simulation. As depicted in Figure 6, the swarm successfully avoids the obstacles and the desired formation is well preserved during the flight. Under our MFG control, the formation similarity error is steadily close to 0 as shown in Figure 7a; the executed trajectories are zero-collision, which can be seen from Figure 7b,c.

4.5. Comparison with Baselines

In the last part of the simulation results, we verify the superiority of our method by comparing it to some typical baseline methods in

Table 2. Comparison with baseline methods.

Method	Scene	Scale of UAVs	Scene Complexity	Communication
[4]	Formation flight and obstacle avoidance	Small	0.67 1	$\mathcal{O}\left(N\right)$ 2
[9]	Cluster flight	Large	0.5	$\mathcal{O}\left(1\right)$
[13]	Cluster flight and obstacle avoidance	Large	0.67	$\mathcal{O}\left(1\right)$
[7]	Formation flight and dense obstacle avoidance	Small	0.83	$\mathcal{O}\left(N\right)$
Ours	Formation flight and dense obstacle avoidance	Large	1	$\mathcal{O}\left(1\right)$

¹ The measurement method of scene complexity is as follows: here it is a scoring system, [cluster flight, formation flight] = [1′, 2′]; [obstacle avoidance, dense obstacle avoidance] = [1′, 2′]; [small, large] = [1′, 2′]. We accumulate the scores for each literature experiment scene according to each item and finally normalize them. ² $\mathcal{O}$() is infinitesimal of the same order.

. From Table 2, it can be seen that our approach simplifies large-scale communication and handles the most complex application scenarios. For more details of these baseline methods, please refer to the related documents [4,7,9,13], etc.

5. Conclusions

In this paper, we propose a large-scale agent trajectory optimization method capable of navigating large-scale swarms in formations while avoiding dense obstacles. First, two types of differentiable mean-field terms are developed, where the overall similarity distance between large-scale 3-D formations and the potential energy value of dense 3-D obstacles are quantified respectively. Second, we embed these two terms into a designed mean-field game control framework that simultaneously satisfies the dynamic constraints and minimizes energy cost, formation similarity error, and collision penalty. This framework solves the communication and computing difficulties caused by complex scenarios of large-scale multiagent systems and formulates the trade-offs between formation and obstacle avoidance in large-scale contexts. Finally, the solid performance of our method in simulations of the formation and dense obstacle avoidance for large-scale UAVs validates its practicality and efficiency. The comparison with baseline methods shows the advanced nature of our method. In the future, we will consider the application and adaption of this mean-field game control framework in more complex scenarios of multiagent systems, such as complex communication conditions (time-varying, partial rejection) and dense mixed obstacles (static, dynamic).