# B-splines-based mechanism design in population games for distributed evolutionary dynamics formation control

Cong Khanh Dinh, Vincent Marguet, Ionela Prodan, Carlos

Ocampo-Martinez, Nicanor Quijano, Julian Barreiro-Gomez

► To cite this version:

Cong Khanh Dinh, Vincent Marguet, Ionela Prodan, Carlos Ocampo-Martinez, Nicanor Quijano, et al.. B-splines-based mechanism design in population games for distributed evolutionary dynamics formation control. Control Engineering Practice, 2026, 168, 10.1016/j.conengprac.2025.106676. hal-05497433

HAL Id: hal-05497433

https://hal.science/hal-05497433v1

Submitted on 6 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# B-splines-based mechanism design in population games for distributed evolutionary dynamics formation control

Cong Khanh Dinh ID a,*, Vincent Marguet ID a, Ionela Prodan ID a, Carlos Ocampo-Martinez ID b, Nicanor Quijano ID c, Julian Barreiro-Gomez ID d

<figure><img src="image_2.png" /></figure>

## ARTICLE INFO

### Keywords:

- B-spline parameterization

- Game theory

- Multi-copter

- Formation reconfiguration

## ABSTRACT

This paper presents a novel formation control framework for autonomous multi-agent systems based on distributed evolutionary dynamics and time-varying population games. We design the fitness functions using B-spline parameterization, enabling smooth, constraint-aware trajectory generation that incorporates both equilibrium targets and transient behavior. The proposed method is integrated with multiple distributed dynamics, namely replicator, projection, and Smith dynamics, and ensures convergence under connectivity and stability assumptions. Compared to conventional static fitness designs, the B-spline-based approach provides enhanced performance in formation reconfiguration and motion planning. The effectiveness of the proposed strategy is validated through extensive simulations and real-world experiments with multiple aerial nanodrones.

## 1. Introduction

Formation control in multi-agent systems poses significant challenges due to decentralization, communication constraints and dynamic environments. Factors such as agent heterogeneity, uncertainties and topology changes further increase complexity. Addressing these issues requires control strategies that are robust, scalable and adaptable to real-world constraints. Existing formation control strategies are commonly classified into two main categories: direct and indirect approaches. Direct methods formulate explicit cost functions and hard constraints, often implemented through techniques such as Model Predictive Control (MPC) (Sforni et al., 2024; Vargas et al., 2022). These methods enable precise enforcement of formation objectives and system limitations within the optimization routine, but typically result in increased computational complexity, which may hinder real-time implementation or scalability. Indirect approaches, on the other hand, encode formation objectives implicitly by introducing penalty or barrier functions into the cost formulation (Miele et al., 2024; Zhang et al., 2024). Rather than enforcing constraints explicitly, these methods steer the system toward the desired formation through soft constraint mechanisms, often reducing the computational burden and facilitating online applicability.

An example of an indirect approach is the use of population games, in which agents interact through local utility functions that depend on the aggregate behavior of the group. In this framework, formation emerges as the outcome of agents optimizing their strategies in response to shared incentives or penalties, without the need for explicitly imposed formation constraints (Quijano et al., 2017). This decentralized and distributed structure makes population games particularly well-suited for scalable and adaptive formation control. Agents iteratively adjust their strategies based on local interactions, following various revision protocols that give rise to different evolutionary dynamics. For example, the pairwise imitation protocol leads to replicator dynamics, the pairwise comparison protocol yields Smith dynamics, the logit-choice protocol generates logit dynamics, and a modified pairwise comparison protocol results in projection dynamics (Barreiro-Gomez & Tembine, 2018; Sandholm, 2010). These evolutionary dynamics serve as distributed optimization processes for resource allocation, where the resulting Nash equilibria coincide with the Karush-Kuhn-Tucker conditions under the full-potential game framework.

Studies have focused on designing revision protocols to generate evolutionary dynamics that maintain properties such as Nash stationarity and positive correlation, ensuring convergence to Nash

<footnote>a Univ. Grenoble Alpes, Grenoble INP†, LCIS, †Institute of Engineering Univ. Grenoble Alpes, Valence, 26000, France</footnote>

<footnote>b Department of Automatic Control, Universitat Politècnica de Catalunya - BarcelonaTECH, Barcelona, 08028, Spain</footnote>

<footnote>c Department of Electrical and Electronics Engineering, Universidad de los Andes, Bogota, 111711, Colombia</footnote>

<footnote>d Department of Computer and Information Engineering, KU Center for Autonomous Robotic Systems, Khalifa University, Abu Dhabi, 127788, UAE</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: cong-khanh.dinh@lcis.grenoble-inp.fr (C.K. Dinh), vincent.marguet@lcis.grenoble-inp.fr (V. Marguet), ionela.prodan@lcis.grenoble-inp.fr (I. Prodan), carlos.ocampo@upc.edu (C. Ocampo-Martinez), nquijano@uniandes.edu.co (N. Quijano), julian.barreirogomez@ku.ac.ae (J. Barreiro-Gomez).</footnote>
equilibria. The work in Barreiro-Gomez et al. (2016b) introduced re-
vision protocols with migration constraints, conditioning interactions
among decision makers by their strategy selection. These constraints are
modeled through a connected, potentially incomplete graph, resulting in
distributed evolutionary dynamics. This approach, combined with full-
potential games, establishes population games as a powerful framework
for distributed optimization. However, a limitation of both centralized
and distributed evolutionary dynamics is their restriction to positivity
and simplex constraints. To address this issue, works such as Barreiro-
Gomez and Tembine (2018), Martinez-Piazuelo et al. (2022a) have pro-
posed revision protocols incorporating both migration constraints and
carrying capacities, expanding the set of feasible constraints. The de-
sign of evolutionary dynamics largely depends on the choice of revision
protocol. Although different protocols may share the same equilibrium
point, they can exhibit different transient behaviors. The question of op-
timal protocol selection was recently explored in Barreiro-Gomez and
Park (2025), where a connection between evolutionary dynamics, pop-
ulation games, and mean-field games has been established.

In contrast to approaches focused on designing revision protocols,
other works have investigated the design of fitness functions and their
influence on the resulting evolutionary dynamics. For instance, the au-
thors of Park and Leonard (2025) have proposed a Kullback-Leibler
divergence-regularized learning scheme to mitigate the effects of de-
layed payoffs, which can otherwise induce undesired oscillations. Simi-
larly, the authors of Martinez-Piazuelo et al. (2022b) have introduced a
class of payoff dynamics tailored to ensure convergence to generalized
Nash equilibria.

With all the aforementioned advances in the field, many works have studied distributed control design, and even machine learning (Barreiro-Gomez & Poveda, 2023), based on distributed evolutionary dynamics. In particular, we focus on the work presented in Barreiro-Gomez et al. (2021), where distributed evolutionary dynamics was used for the multi-agent formation control. In this work, the population game is designed exclusively and only taking care of the equilibrium point, but there was no design on the transitory behavior, i.e., the fitness function depends only on the population states. Building on Barreiro-Gomez et al. (2021), we enhance the design of population games using B-splines (Lyche et al., 2018), allowing the fitness function to incorporate both the desired equilibrium and an optimal time-dependent trajectory, thus improving the distributed formation algorithm.

This aspect is where B-splines become highly useful. Unlike piecewise-linear or discrete approximations, B-splines ensure that small changes in the strategy space lead to smooth adjustments, reflecting gradual adaptation in real-world populations. One of the main benefits of B-splines is their ability to represent smooth functions using a small number of control points. B-splines are widely used in trajectory generation and optimization due to their desirable mathematical properties, such as local support, unity partition, smoothness, and the convex hull property (Lyche et al., 2018). These properties make B-splines particularly useful for representing smooth and flexible curves while maintaining computational efficiency. In trajectory planning, B-splines provide a structured way to parametrize motion while ensuring continuity and controllability. Their recursive definition allows for efficient evaluation, and their differentiability properties enable analytical computation of velocities and accelerations, which is crucial for constrained motion planning (Lyche et al., 2018; Marguet et al., 2024b).

In this paper, we enhance the formation control approach proposed in Barreiro-Gomez et al. (2016a) with a data-driven evolutionary game theory extension presented in Barreiro-Gomez et al. (2021). These works focus solely on steady-state design and do not address transitory events, meaning they cannot guarantee collision-free formations or account for other physical or operational constraints. Additionally, they rely on the assumption that the local controller can accurately track the expected formation of follower agents, which is often unrealistic in practice. Finally, the implementation presented in Barreiro-Gomez et al. (2021) was limited to two physical followers and one virtual leader, which restricted

the distributed communication topology to the path graph. Here, we present experimental results and validation that comprise four physical agents. Then, we may test different distributed graph topologies.

While these approaches show promise, in this paper, we address their limitations by integrating a trajectory planner into the population game framework that leverages B-splines properties. This integration refines the fitness function into a time-varying condition, making it more suitable for real-time applications. Simulations and experiments<fnref n="1" /> are conducted with multiple nanodrones for formation control and reconfiguration to validate the performance of the proposed method.

Therefore, the main contributions of this work are as follows:

* We develop an online formation control framework based on distributed population games, and establish its stability through theoretical analysis and experimental validation with multiple aerial nanodrones.

* We design a novel time-varying fitness function using B-spline parameterization, which enables smooth transitions between formations while inherently satisfying motion constraints.

* We conduct a comparative analysis of multiple revision protocols in distributed settings, quantifying their effects on implementation complexity, convergence speed, and overall control performance.

The remainder of this paper is organized as follows. Section 2 provides background on distributed population games and B-spline properties, which form the foundation of the proposed method. Section 3 presents the main formation control framework, including the time-varying fitness function design, stability analysis, and integration with tracking control. Section 4 validates the proposed approach through simulation and experimental scenarios involving multiple aerial nanodrones. Finally, Section 5 concludes the paper and discusses directions for future research.

Throughout this work, the following notation is employed. Bold letters signify vectors, whereas matrices are indicated by bold capital letters. The set of real, positive, and non-negative real numbers is denoted by $\mathbb{R}$, $\mathbb{R}_{>0}$, $\mathbb{R}_{+}$, respectively. Without loss of generality, $p_i^x$ refers to the state of strategy i-th in the x population. The notation $[\mathbf{P}]_j$ is used to denote the j-th point in the matrix of control points $\mathbf{P} = [\mathbf{P}_1, \mathbf{P}_2, \dots]$. The composition of functions is denoted by $\circ$. Finally, operator $[\cdot]_+ \triangleq \max(\cdot, 0)$. The Euclidean norm of a vector $\mathbf{x}$ is denoted by $\|\mathbf{x}\|$.

## 2. Preliminaries

This section introduces the graph characterization within the multi-agent setting, the evolutionary dynamics (Barreiro-Gomez et al., 2016b; Sandholm, 2010) and further the B-splines properties (Lyche et al., 2018) that will be instrumental in the fitness function design.

### 2.1. Settings and problem statement

Let $\mathcal{V} = \{1, 2, \dots, \nu\}$ be the set of $\nu > 2$ agents. One of them is defined as the leader $\ell \in \mathcal{V}$. The leader will determine the behaviors of the remaining $\nu - 1$ followers defined in $\mathcal{F} = \mathcal{V} \setminus \{\ell\}$. Each agent $i \in \mathcal{V}$ is located at position $\mathbf{z}_i = [x_i \ y_i \ z_i]^\top \in \mathbb{R}_+^3$, which is assumed to be measurable and non-negative<fnref n="2" /> We assume that each agent $i \in \mathcal{V}$ has a limited communication capability denoted by $\rho_i \in \mathbb{R}_+$. For simplicity, we consider that all agents have identical communication capabilities, i.e., $\rho_i = \rho_j$ for all $i, j \in \mathcal{V}$. As a result, the communication network evolves over time according to the spatial distribution of the agents. The connectivity of the system is represented by a weighted graph $\mathcal{G}(t) = (\mathcal{V}, \mathcal{E}(t), \mathbf{A}(t))$, where

$$ \mathcal{E}(t) = \{ (i, j) \mid i, j \in \mathcal{V}, d_{ij}(t) \le \rho_i \} \qquad (1) $$

<footnote n="1">¹ The link for the experimental video is GTB spline.</footnote>

<footnote n="2">² Note that this is extensible to consider more parameters for the agents in a straightforward manner, e.g., $\mathbf{z}_i = [x_i \ y_i \ z_i \ \phi_i]^\top$, including the yaw angle.</footnote>
is the set of edges that describes the communication topology among the nodes representing the agents where $d_{ij}(t) = \|z_i - z_j\|$, with $\| \cdot \|$ denoting the Euclidean norm. Besides, $A(t) = [a_{ij}]$ is a $v \times v$ adjacency matrix whose entries satisfy

$$a_{ij} = \begin{cases} 1, & \text{if } (i, j) \in \mathcal{E}, \\ 0, & \text{otherwise.} \end{cases} \qquad (2)$$

From this point, we define $\mathcal{N}_i = \{j | (i, j) \in \mathcal{E}\}$ as the collection of neighbors of a node $i \in \mathcal{V}$, meaning the set of agents that the $i$-th agent can communicate with. In this formation control problem, the leader follows a precomputed reference trajectory while each follower $i \in \mathcal{F}$ tracks the leader by satisfying a desired formation given by a reference distance on each axis, namely $r_i = [r_i^x \ r_i^y \ r_i^z]^T$. Given potential non-complete communication nature in the formation, the distributed population games discussed below will tackle this issue.

## 2.2. Population game

We begin by reviewing key concepts from distributed population games, as introduced in Barreiro-Gomez et al. (2021). Fundamentally, a population is characterized by a total mass $m \in \mathbb{R}_+$, representing a large but finite number of players engaged in strategic interactions. Each player has a set of $v$ available strategies given by $\mathcal{V} = \{1, 2, ..., v\}$. Note that this set corresponds to the set of agents in the formation problem, i.e., there is a strategy corresponding to each agent. A non-negative scalar $p_i \in \mathbb{R}_+$ represents the fraction of population mass that chooses the strategy $i \in \mathcal{V}$. The population state is denoted by vector $\mathbf{p} = [p_1 \ p_2 \ ... \ p_v]^T \in \mathbb{R}_+^v$, which represents the strategic distribution. Due to limitation on the resource $m$, the admissible set of population states<fnref n="3" /> is a simplex

$$\Delta = \left\{ \mathbf{p} \in \mathbb{R}_+^v \mid \sum_{i \in \mathcal{V}} p_i = m \right\}. \qquad (3)$$

**Definition 1** (Fitness function and Population games). Players playing the strategy $i \in \mathcal{V}$ will obtain a reward defined as fitness function $f_i : \Delta \to \mathbb{R}$. A population game thus can be fully characterized by the fitness vector $\mathbf{f} = [f_1 \ f_2 \ ... \ f_v]^T$.

**Definition 2** (Full potential game). Let $\mathbf{f} : \mathbb{R}_+^v \to \mathbb{R}^v$ be a population game. If there exists a continuously differentiable potential function $\Phi(\mathbf{p})$ such that $\nabla\Phi(\mathbf{p}) = \mathbf{f}(\mathbf{p})$, $\forall \mathbf{p} \in \mathbb{R}_+^v$, then $\mathbf{f}$ is a full potential game.

**Definition 3** (Nash equilibrium Sandholm, 2010). The point $p^* \in \Delta$ is a Nash equilibrium of the population game if $f_i(p^*) \ge f_j(p^*)$, $\forall i, j \in \mathcal{V}$. The set of Nash equilibria is thus defined as $\mathbb{N}\mathbb{E}(\mathbf{f}) = \{p \in \Delta | p_i > 0 \Rightarrow f_i(p) \ge f_j(p), \forall i, j \in \mathcal{V}\}$.

Decision makers reconsider their choices of strategies based on the current fitness function and aggregate the behavior of the overall system aiming for improvement of their benefits. Typically, it is not inherent to justify the assumption of Nash equilibrium knowledge in the context of large number of agents. This feature leads to the idea of modeling this evolution using a notion called revision protocol.

**Definition 4** (Revision protocol). The function $\eta : \mathbb{R}^v \times \Delta \to \mathbb{R}_+^{v \times v}$ is known as revision protocol describing the timing and decision of players in strategic interaction. The element $\eta_{ij}(f_i, f_j, p_i, p_j)$ denotes the conditional switching rate from strategy $i$ to $j$ given the fitness reward $f_i, f_j$ and the population states $p_i, p_j$. To simplify the notation, we denote $\eta_{ij}(f_i, f_j, p_i, p_j) = \eta_{ij}$.

Incorporating a communication-aware revision protocol and population games in the mean dynamics (Sandholm, 2010), which corresponds to the Kolmogorov-Forward equation as shown in Barreiro-Gomez and Park (2025), and adapting to a multi-agent setting with limited information exchange, we obtain the following distributed population dynamics:

$$\dot{p}_i = \sum_{j \in \mathcal{N}_i} p_j \eta_{ji} - p_i \sum_{j \in \mathcal{N}_i} \eta_{ji}, \quad \forall i \in \mathcal{V} \qquad (4)$$

In this work, various distributed population dynamics are recalled from Barreiro-Gomez et al. (2016a) and Barreiro-Gomez and Tembine (2018), which can then be integrated into the control scheme for multi-agent application presented in the next section.

**Definition 5** (Replicator dynamics, Barreiro-Gomez et al., 2016b). Choosing the well-known replicator revision protocol for the population dynamics model

$$\eta_{ij} = p_j [f_j - f_i]_+, \qquad (5)$$

we can obtain the Distributed Replicator Dynamics (DRDs)

$$\mathbf{p} = \mathbf{L}^{(\mathbf{p})} \mathbf{f}(\mathbf{p}), \qquad (6)$$

where $L^{(\mathbf{p})} = (l_{ij})^{(\mathbf{p})}$, and

$$l_{ij}^{(\mathbf{p})} = \begin{cases} -a_{ij}p_i p_j & \text{if } i \neq j, \\ \sum_{k \in \mathcal{V} \setminus \{i\}} a_{ik} p_i p_k & \text{if } i = j. \end{cases} \qquad (7)$$

**Definition 6** (Distributed Smith dynamics (DSDs)). Applying pairwise comparison protocol to the distributed mean dynamics

$$\eta_{ij} = [f_j - f_i]_+, \qquad (8)$$

the DSDs are defined as follows:

$$\bar{\mathbf{p}} = \tilde{\mathbf{L}}^{(\mathbf{p})} \mathbf{f}(\mathbf{p}), \qquad (9)$$

where $\tilde{\mathbf{L}}^{(\mathbf{p})} = \tilde{\eta}_{ij}^{(\mathbf{p})}$, and

$$\tilde{\eta}_{ij}^{(\mathbf{p})} = \begin{cases} -\frac{1}{2} a_{ij}(p_i(1 - \phi_{ij}) + p_j(1 + \phi_{ij})) & \text{if } i \neq j, \\ \frac{1}{2} \sum_{k \in \mathcal{V}} a_{ik}(p_i(1 - \phi_{ik}) + p_k(1 + \phi_{ik})) & \text{if } i = j, \end{cases} \qquad (10)$$

with $\phi_{ij} = \text{sign}(f_i - f_j)$, $\text{sign}(\cdot)$ denoting the sign function.

**Definition 7** (Projection dynamics). By introducing the modified pairwise comparison protocol in the mean dynamics

$$\eta_{ij} = \frac{[f_j - f_i]_+}{p_i}, \qquad (11)$$

the distributed Projection dynamics (DPDs) are given by

$$\mathbf{p} = \mathbf{L}\mathbf{f}(\mathbf{p}), \qquad (12)$$

where $\mathbf{L}$ is the Laplacian of the graph $\mathcal{G}$.

Hence, $\tilde{\mathbf{L}}^{(\mathbf{p})}$, $\mathbf{L}^{(\mathbf{p})}$, and $\mathbf{L}$ represent the Laplacian matrices associated with the distributed Smith, replicator, and projection dynamics, respectively. These dynamics are particularly relevant to our setting, as the proposed formation control framework is designed to operate online under a time-varying communication topology.

**Theorem 1** (Asymptotic stability of DSDs, DRDs and DPDs). Let $\mathbf{f}$ be a full potential game with strictly concave potential function $\Phi(\mathbf{p})$, that is, $\mathbf{f}(\mathbf{p}) = \nabla\Phi(\mathbf{p})$, and let $\mathbf{p}^* \in \mathbb{N}\mathbb{E}(\mathbf{f})$. If the graph is connected and $\mathbf{p}^* \in \text{int}\Delta$, then $\mathbf{p}^*$ is asymptotically stable. (cf. Barreiro-Gomez et al., 2016b)

**Corollary 1** (Time-varying graph). The asymptotic stability of $\mathbf{p}^*$ is also valid for the time-varying graph $\mathcal{G}(t) = (\mathcal{V}, \mathcal{E}(t), \mathbf{A}(t))$ as discussed in Barreiro-Gomez et al. (2016a).

With its versatility, the Corollary 1 can be used to a wide range of systems that share the same properties of the connected graph.

For further use in the evolutionary dynamics formation control, the splines properties are next introduced.

<footnote n="3">³ For introducing the theoretical concepts, in this section, we consider a single population. For the formation control in the next section, multiple populations will be employed.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Distributed formation control scheme for multi-drone system.</figcaption></figure>

## 2.3. B-splines

Concepts related to B-splines will also be beneficial for the design of constrained reference trajectories (Lyche et al., 2018).

A trajectory denoted by $z(t) \in \mathbb{R}^3$ is parameterized as a linear combination of control points, $\mathbf{P} \in \mathbb{R}^{3 \times n}$

$$ \mathbf{P} = [\mathbf{P}_0 \dots \mathbf{P}_{n-1}], \qquad (13) $$

and B-splines basis functions $\mathbf{B}_{d,\tau}(t) = [B_{0,d,\tau}(t) \dots B_{n-1,d,\tau}(t)]^T$

$$ \mathbf{z}(t) = \sum_{k=0}^{n-1} \mathbf{P}_k \mathbf{B}_{k,d,\tau}(t) = \mathbf{P} \mathbf{B}_{d,\tau}(t), \quad \forall t \in [0, T_f], \qquad (14) $$

with $d$ the degree and $\tau = \{\tau_0 \le \tau_1 \le \dots \le \tau_{n+d}\}$ a knot sequence starting at 0 and ending at $T_f$.

The basis functions are defined recursively and possess various mathematical and geometric properties, such as local support, partition of unity, smoothness, and the convex hull property, as discussed in Marguet et al. (2024a). Specifically, the smoothness property refers to the fact that the r-th order derivatives of a B-spline basis function can be expressed as linear combinations of B-spline basis functions of lower degree. By selecting a uniform clamped knot vector, continuity of derivatives up to order $d-1$ is ensured. This characteristic enables the representation of the derivative of a B-spline trajectory as a B-spline curve of reduced degree. Consequently, the control points of the derivative curve can also be deduced from the original control points via a conversion matrix, denoted by $\mathbf{M}_r$, i.e.,

$$ \mathbf{P}^{(r)} = \mathbf{P}\mathbf{M}_r, \qquad (15) $$

where $r$ represents the degree of derivation, which gives the following interpretation of the derivative curve:

$$ \mathbf{z}^{(r)}(t) = \mathbf{P}\mathbf{M}_r\mathbf{B}_{d-r,\tau}(t), \quad \forall t \in [0, T_f]. \qquad (16) $$

In trajectory planning, it is common to formulate and solve a constrained optimization problem, such as

$$ \min_{\mathbf{z}(t)} \int_{0}^{T_f} C(\mathbf{z}(t))dt \qquad (17a) $$

$$ \text{s.t.} \quad g_i(\mathbf{z}(t)) \le 0, \quad 0 \le t \le T_f, \qquad (17b) $$

$$ h_j(\mathbf{z}(t)) = 0, \quad 0 \le t \le T_f. \qquad (17c) $$

In (17c), we find a trajectory $\mathbf{z}(t)$ as in (14) that minimizes the cost $C(\cdot)$ from (17a) while simultaneously satisfying (17b) and (17c) over the time interval $[0, T_f]$:

* the cost $C(\cdot)$ can be taken as the energy consumed along the time horizon $[0, T_f]$ and should be minimized;

* the inequalities (17b) may consider bounds on the thrust magnitude and torques;

* the equalities in (17c) may represent constraints for passing through way-points.

Next, all these elements are merged to formulate the distributed formation control problem with off-line and on-line phases.

## 3. Distributed formation control

The proposed framework builds on population game dynamics, which naturally align with distributed consensus schemes and enable agents to coordinate locally to achieve a shared formation objective (Obando et al., 2013, 2022). By embedding the fitness functions in a B-spline representation, the approach produces smooth and constraint-aware trajectories with formal stability guarantees. Because the resulting dynamics reduce to first-order nonlinear differential equations (i.e., a zero type learning process in the sense of Vidal, 2001), the implementation remains lightweight and scalable, requiring only local information for real-time updates.

Population game dynamics provide the core mechanism for our formation control framework, enabling each UAV to adapt its behavior based on locally computed fitness values while B-spline parameterization ensures smooth constraint-aware trajectories, combining scalability, distributed coordination, and real-time implementability in a single approach.

The proposed scheme comprises both off-line and on-line stages as shown in Fig. 1. In the offline stage, the formation generation using B-spline considers system constraints and the desired formation shape to produce the leader reference, along with the time-varying gaps between followers and the leader $\mathbf{r}_i$ via B-spline curves. From the game theoretical perspective, this off-line procedure corresponds to a population game mechanism design able to incorporate physical and operational constraints over the system dynamics. On-line, the generated trajectories are used as the inputs for the distributed evolutionary dynamics to address the B-spline designed population game. As the dynamics approaches the Nash equilibrium, each agent acquires the necessary reference for the control loop. Utilizing the tracking controller equipped with a safety filter, developed according to Marguet et al. (2024b), the formation of many drones will be able to execute the intended configuration in real time.
This section introduces first the evolutionary dynamics which will be tested for the distributed formation control. The particularity resides in the fitness function design that is time-varying due to the use of B-splines. Stability proofs are based on Lyapunov and input-to-state stability (ISS) arguments.

### 3.1. Dynamical model and references generation

In this work, we consider the quadcopter model. The platform for testing is the nanodrone Crazyflie, which will be used for validation of the formation control framework. This section recalls the translational dynamics of a conventional quadcopter (Cao & Lynch, 2021)

$$ \ddot{z} = -ge_3 + \mathbf{R}_{\psi}\Omega(\mathbf{u}), \quad (18) $$

$$ \Omega(\mathbf{u}) = \begin{bmatrix} T \cos \phi \sin \theta \\ T \sin \phi \\ T \cos \phi \cos \theta \end{bmatrix}, \quad (19) $$

$$ \mathbf{R}_{\psi} = \begin{bmatrix} \cos \psi & \sin \psi & 0 \\ \sin \psi & -\cos \psi & 0 \\ 0 & 0 & 1 \end{bmatrix}, \quad (20) $$

where $\mathbf{z} \in \mathbb{R}^3$ is the position of the drone in the inertial frame, $\mathbf{u} = [T, \phi, \theta]^\top \in \mathbb{R}^3$ denotes the control input, which gathers the normalized thrust $T$, the roll angle, $\phi$ and the pitch angle $\theta$. Moreover, $\mathbf{R}_{\psi}$ is the rotation matrix that performs the conversion from the inertial reference frame to the quadcopter's body frame that has zero roll and pitch relative to the inertial frame given the measured yaw angle $\psi$, and $g = 9.81 \text{ m/s}^2$ is gravitational acceleration. The exact linearized model is obtained by using flat representation as described in Do et al. (2021) and Nguyen et al. (2021) given as follows:

$$ \dot{\xi} = H\xi + B\mu, \quad (21) $$

where $\xi = [\mathbf{z}^\top \quad \mathbf{v}^\top]^\top$ gathers the position $\mathbf{z}$ and velocity $\mathbf{v} = \dot{\mathbf{z}}$ of the drone. Besides, $H = \begin{bmatrix} \mathbf{0}_3 & \mathbf{I}_3 \\ \mathbf{0}_3 & \mathbf{0}_3 \end{bmatrix}$, $B = \begin{bmatrix} \mathbf{0}_3 \\ \mathbf{I}_3 \end{bmatrix}$, and $\mu = [\mu_x \quad \mu_y \quad \mu_z]^\top \equiv \mathbf{R}_{\psi}\Omega(\mathbf{u}) - ge_3$ denote the virtual input in the flat-output space.

The reference trajectory of the leader, namely $c_\ell(t)$, is supposed to be known a priori and can be obtained for example with the algorithm presented in Marguet et al. (2024a). This fact allows to impose constraints on the position, velocity, thrust, angles and angular velocity of the leader.

### 3.2. Fitness function design based on B-spline

In this work, B-spline curves are integrated into the fitness function of the population game to transform it into a time-varying condition, enabling real-time trajectory planning and reconfiguration. This enhancement allows the formation control approach to consider transient behaviors, ensuring collision avoidance and accommodating physical or operational constraints. Hence, to the best of our knowledge, it improves upon prior methods that focused solely on steady-state formations.

Consider a multi-population system involving 3 populations denoted by the set $Q = \{x, y, z\}$. For a population $q \in Q$, the value $p_i^q$ represents the reference coordinate on $q$-axis of the $i$-th strategy (agent in the formation control problem). Thus, vector $\mathbf{p}^q = [p_1^q \ p_2^q \ \dots \ p_\ell^q]^\top \in \mathbb{R}_+^\nu$ provides information about the strategic distribution of all the decision makers in the population $q$. Consequently, $\mathbf{p}_i = [p_i^x \ p_i^y \ p_i^z]^\top \in \mathbb{R}_+^3$ is the total desired position of the $i \in \mathcal{F}$ follower agent. Without losing generality and for simplicity, if the population is not mentioned, it can be understood as one population and the superscript $q$ is omitted in the notation.

**Proposition 1** (Leader-independent population mass Barreiro-Gomez et al., 2016a). Given one of distributed dynamics by (6), (9), (12), a fitness function can be designed so that the specified formation corresponds to Nash equilibrium. Such a fitness function would be

$$ f_i^q = \begin{cases} r_i^q - p_i^q, & \forall i \in \mathcal{F}, \\ -c_\ell^q, & i = \ell, \end{cases} \quad (22) $$

where

* $c_\ell^q$ is the current reference position of the leader for population $q \in Q = \{x, y, z\}$,

* $\mathbf{r}_i = [r_i^x \ r_i^y \ r_i^z]^\top$ represents the reference gap distance for agent $i \in \mathcal{F}$ w.r.t leader.

**Assumption 1.** Within the game-theoretic framework, and without loss of generality, $c_\ell^q(t)$ is assumed to be strictly positive for all $t > 0$. This assumption can be satisfied by applying a suitable offset to the original trajectory.

**Proof.** Given the fitness function designed in (22), this is a full potential game. Indeed, let us choose a concave potential function

$$ \Phi(\mathbf{p}) = -\frac{1}{2} \sum_{i \in \mathcal{F}} (r_i - p_i)^2 - c_\ell p_\ell, \quad (23) $$

where $\mathbf{p} \in \mathbb{R}_+^\nu$. To verify whether (23) is a potential function, we need to show that the change in the potential function corresponds to the fitness function for any strategy, i.e.,

$$ \nabla \Phi(\mathbf{p})_i = r_i - p_i, \quad \forall i \in \mathcal{F}, \quad (24) $$

$$ \nabla \Phi(\mathbf{p})_\ell = -c_\ell, \quad i = \ell. \quad (25) $$

Since this condition is satisfied, the game is a potential game with the potential function (23). Thus, the given fitness function defines a potential game, and this population dynamics satisfies the asymptotic stability for DRDs, DSDS, and DPDs according to Theorem (1). For more details, see (Barreiro-Gomez et al., 2016b), Section 3. □

Note that, according to (3), we have $\sum_{i \in \mathcal{F}} p_i + p_\ell = m$. Assume that $m$ is selected such that the maximum point in (23) is attained, i.e., $\mathbf{p}^* \in \text{int}\Delta$ (relative interior), and $\nabla\Phi(\mathbf{p})_i = \nabla\Phi(\mathbf{p})_\ell = 0$, for all $i \in \mathcal{F}$. Then, as $p_i \ge 0$, for all $i \in \mathcal{V}$, we have

$$ m = \sum_{i \in \mathcal{F}} r_i + (v - 1)c_\ell + p_\ell, \quad (26) $$

to guarantee that $\mathbf{p}^* \in \text{int}\Delta$,

$$ m > \sum_{i \in \mathcal{F}} r_i + (v - 1)c_\ell, \quad (27) $$

which is the minimal required mass.

The work in Barreiro-Gomez et al. (2016a) proposed to design the fitness functions depending only on the population states and a known measurable position for the leader, in addition to the time independence. In such a design, the transitory behavior of the population dynamics was not under any control, while only guarantees on the convergence point were considered in the design. Moreover, constraints related to the system were not considered, but only constraints over the population states could be satisfied. The work presented here overcomes these limitations by using the B-splines for the time-dependent fitness design.

In (22), the fitness function is designed with asymptotic stability for fixed formation. We extend this notion to a time-varying formation where the reference of the leader and the gaps relative to the leader are time-dependent. This aspect is quite beneficial for switching formation tasks.

First, let us define the $\mathcal{L}_\infty$ norm and some results of the Input-to-State Stability (ISS) framework.

**Definition 8** ($\mathcal{L}_\infty$ norm). Consider a function $u: [0, \infty) \to \mathbb{R}^m$. The so-called supremum norm of this function is defined as

$$ ||u||_\infty = \sup_{t_0 \le \tau \le t} ||u(\tau)||. \quad (28) $$

The set of all such functions satisfying (28) is denoted by $\mathcal{L}_\infty$, where $|| \cdot ||$ stands for the Euclidean norm.

Let us consider the nonlinear system

$$ \dot{x}(t) = f(x(t), u(t)), \quad t \ge 0, \quad (29) $$

where $x(t) \in \mathbb{R}^n$, $u(t) \in \mathbb{R}^m$ for some $n, m \in \mathbb{N}$.
**Definition 9** (Comparison functions Khalil & Grizzle, 2002). A function $\alpha: [0, \infty) \to [0, \infty)$ is said to belong to a class $\mathcal{K}$ if it is strictly increasing and $\alpha(0) = 0$. Furthermore, function $\alpha \in \mathcal{K}_\infty$ if it belongs to $\mathcal{K}$ and $\alpha(r) \to \infty$ as $r \to \infty$.

Similarly, a function $\beta: [0, \infty) \times [0, \infty) \to [0, \infty)$ is said to belong to class $\mathcal{KL}$ if it is continuous, $\beta(., t)$ strictly increasing, $\beta(r, .)$ strictly decreasing, and $\beta(r, t) \to 0$ as $t \to \infty$.

**Definition 10** (Input-to-state stability Sontag & Wang, 1995). The system (29) is said to be input-to-state stable (ISS) if there exists $\beta \in \mathcal{KL}$ and $\gamma \in \mathcal{K}_\infty$ such that $\forall x(t_0)$ and $u(t)$ is bounded,

$$||x(t)|| \le \beta(||x(t_0)||, t - t_0) + \gamma(||u||_\infty). \quad (30)$$

**Remark 1.** In Definition 10, the $\gamma$ is a gain for one of upper bound in of ISS, which is not unique.

**Theorem 2.** (Khalil & Grizzle, 2002; Sontag, 2008) The system (29) is ISS if and only if there exists a smooth function $V(x) : \mathbb{R}^n \to \mathbb{R}_+$ continuously differentiable such that

$$\alpha_1(||x||) \le V(x) \le \alpha_2(||x||), \quad (31)$$

$$\frac{\partial V}{\partial x} f(x, u) \le -W_3(x), \quad \forall ||x|| \ge \zeta(||u||) > 0, \quad (32)$$

$\forall x \in \mathbb{R}^n, u \in \mathbb{R}^m, \alpha_1, \alpha_2 \in \mathcal{K}_\infty, \zeta \in \mathcal{K}$ and $W_3(x)$ positive definite. Moreover, $\gamma$ is referred to as the Lyapunov gain function. Subsequently, the system (29) is input-to-state stable (ISS) with gain $\gamma = \alpha_1^{-1} \alpha_2 \zeta$.

**Proposition 2.** [Time-varying fitness function design based on B-spline] Represent the fitness in (22) under the B-spline-generated references given in the previous section, i.e.,

$$f_i^q = \begin{cases} r_i^q(t) - p_i^q, & \forall i \in \mathcal{F}, \\ -c_{\ell}^q(t), & i = \ell. \end{cases} \quad (33)$$

Here, $q \in \{x, y, z\}$ denotes the population for each axis. Furthermore, the leader trajectory is $c_\ell(t) = [c_\ell^x(t) \ c_\ell^y(t) \ c_\ell^z(t)]^\top$, which is further described as

$$c_{\ell}(t) = \mathbf{P}_{\ell} \mathbf{B}_{d,\ell}(t), \quad (34)$$

with $\mathbf{P}_\ell \in \mathbb{R}^{3\times n}$ as in (13) and the basis functions $\mathbf{B}_{d,\ell}(t)$ as in (14). This approach aims to manage follower trajectories in a distributed manner, allowing the formation to be reconfigured dynamically along the path. Each follower is assigned a time-varying gap profile defined $\mathbf{r}_i(t) = [r_i^x(t) \ r_i^y(t) \ r_i^z(t)]^\top$ which is further described by

$$\mathbf{r}_i(t) = \mathbf{P}_i \mathbf{B}_{d,\ell}(t), \quad (35)$$

representing its expected relative position, or offset, with respect to the leader. This gap remains constant when the formation is stable, but evolves during transitions. The B-spline curve $\mathbf{r}_i(t)$, which defines this transition for agent $i$ across all three axes, is determined by solving the following optimization problem:

$$\min_{\mathbf{r}_i(t)} \int_{t_0}^{T_f} ||\dot{\mathbf{r}}_i(t)||^2 dt, \quad (36a)$$

$$\text{s.t.} \quad \dot{\mathbf{r}}_i(t_0) = \dot{\mathbf{r}}_i(T_f) = \mathbf{0}_{3 \times 1}, \quad (36b)$$

$$\dot{\mathbf{r}}_i(t_0) = \dot{\mathbf{r}}_i(T_f) = \mathbf{0}_{3 \times 1}, \quad (36c)$$

$$\mathbf{r}_i(t_k) = \mathbf{p}_{i,k}^{wp}, \quad k = \{0, 1, \dots, n_{wp} - 1\}, \quad (36d)$$

$$\underline{\mathbf{r}} \le \mathbf{r}_i(t) \le \bar{\mathbf{r}}, \quad (36e)$$

where $T_f$ is the final time of the trajectory, $\mathbf{p}_{i,k}^{wp}$ represents the pass-through way-points, and $n_{wp}$ is the total number of way-points. $\underline{\mathbf{r}}$ and $\bar{\mathbf{r}}$ denote the lower and upper bounds for the relative distance between followers and leader.

The cost (36a) corresponds to the square of the velocity magnitude, (36b) and (36c) correspond, respectively, to the initial and final acceleration and velocity. Next, (36d) is the way-point constraint and (36e) ensures that the relative gaps are constrained within a box around the leader, which is beneficial in maintaining the connectivity of the overall

graph, hence allowing distributed online application of the population dynamics.

The optimization problem in (36) can be rewritten as a constrained problem in terms of control points as follows:

$$\min_{\mathbf{P}_i} \sum_{j,k} (\mathbf{P}_i M_1)_j^\top \Gamma(t) (\mathbf{P}_i M_1)_k \quad (37a)$$

$$\text{s.t.} \quad [\mathbf{P}_i M_2]_0 = [\mathbf{P}_i M_2]_{T_f} = \mathbf{0}_{3 \times 1}, \quad (37b)$$

$$[\mathbf{P}_i M_1]_0 = [\mathbf{P}_i M_1]_{T_f} = \mathbf{0}_{3 \times 1}, \quad (37c)$$

$$\mathbf{P}_i \mathbf{B}_d(t_k) = \mathbf{p}_k^{wp}, \quad k = \{0, 1, \dots, n_{wp} - 1\}, \quad (37d)$$

$$\underline{\mathbf{r}} \le \mathbf{P}_i \mathbf{B}_d(t) \le \bar{\mathbf{r}}, \quad (37e)$$

with $\Gamma(t) = \int_{t_0}^{t} n_{wp}^{-1} (B_{j,d-1}(t) B_{k,d-1}(t)) dt$.

To draw a connection between the population dynamics and a control-theoretic framework, we highlight the structural similarity with a non-autonomous system formulation. This fact leads to the following remark.

**Remark 2.** The population dynamics given in (4) can be identified similarly to the non-autonomous system given in (29), where the state $x$ corresponds to the population state $\mathbf{p}$ and the control input $u$ corresponds to the collection of generated trajectories $\mathbf{r}_i(t)$ and $c_\ell(t)$, as defined in the fitness function in (33), i.e.,

$$\dot{\mathbf{p}} = g(\mathbf{p}, \bar{\mathbf{r}}), \quad (38)$$

with

$$\bar{\mathbf{r}}_i = \begin{cases} r_i & i \in \mathcal{F}, \\ c_\ell & i = \ell. \end{cases} \quad (39)$$

**Lemma 1** (Input-to-state stability of population dynamics for time-varying fitness design). Given the fitness function proposed in (33) with the B-spline reference given in the offline phase for the leader and the gaps in (37), and assuming that the graph is connected along the formation, the distributed population dynamics (38) for DRDs, DSDs, and DPDs are ISS.

**Proof.** To establish the proof of ISS, we aim to verify the two conditions outlined in Theorem (2). We choose the Lyapunov candidate inspired from the potential function in (23) as follows:

$$\begin{aligned} E(\mathbf{p}) &= -\Phi(\mathbf{p}), \\ &= \frac{1}{2} \sum_{i \in \mathcal{F}} (r_i(t) - p_i)^2 + c_\ell(t) p_\ell. \end{aligned} \quad (40)$$

First, the condition (31) ensures that the Lyapunov function is positive definite and radially unbounded. Since the fixed version remains unchanged, this condition is already satisfied.

Second, (32) corresponds to following condition:

$$\frac{\partial E}{\partial \mathbf{p}} g(\mathbf{p}, \bar{\mathbf{r}}) \le -W_3(\mathbf{p}), \quad \forall ||\mathbf{p}|| \ge \zeta(||\bar{\mathbf{r}}||) > 0, \quad (41)$$

$\forall \mathbf{p} \in \mathbb{R}^n, \bar{\mathbf{r}} \in \mathbb{R}^n, \zeta \in \mathcal{K}$ and $W_3(\mathbf{p})$ positive definite. Consider $||\mathbf{p}|| = \sqrt{\sum_{i \in \mathcal{F}} p_i^2}$. Applying the Cauchy-Schwarz inequality for $p_i \in \mathbb{R}_+$ yields

$$||\mathbf{p}|| \ge \sqrt{\frac{(\sum_{i \in \mathcal{F}} p_i)^2}{\nu}}. \quad (42)$$

The states are constrained to the condition given (3) since the resource of the population is finite, thus $\sum_{i \in \mathcal{F}} p_i = m$. Here, $m$ can be chosen such that

$$m = \nu \max(\|\underline{\mathbf{r}}\|, \|\bar{\mathbf{r}}\|) + \sqrt{\nu} \sup_{q \in \{\text{x}, \text{y}, \text{z}\}, j \in \{0, \dots, n-1\}} |P_{\ell,j}^q|, \quad (43)$$

where $P_{\ell,j}^q \in \mathbb{R}_+$ denotes the $j$-th control point for the leader $\ell$ in the axis $q$, the values $\underline{\mathbf{r}}$ and $\bar{\mathbf{r}}$ denote the lower and upper bounds for the relative distance between followers and leader. Note that this selection for $m$ in (43) is suitable across all the populations (all axes $x, y, z$) as it uses the maximum leader's control point for all the axes, and the norm for the B-spline curves' bounds as shown in (36).
Describing the trajectories of the agents using B-spline curves allows the inter-agent distance to also be represented as a scalar B-spline curve. Therefore, the control points of this curve can be used to lower and upper bound the actual distance since we are using the convexity property of B-spline curves.

Additionally, we can choose the gain function simply as $\zeta(||\bar{r}||) = ||\bar{r}||$. From the constraints of the B-spline given in (37) and the choice of $m$ above (43), the following inequality can be readily deduced:

$$ \begin{aligned} \|\mathbf{p}\| &\ge \frac{m}{\sqrt{v}} \\ &= \sqrt{v \left( \max(||\bar{r}||, ||\bar{r}||) + \frac{1}{\sqrt{v}} \sup_{q \in \{x,y,z\}, j \in \{0,...,n-1\}} |P_{\varepsilon,j}^q| \right)^2} \\ &\ge \sqrt{v \left( \max(||\bar{r}||, ||\bar{r}||) \right)^2 + \left( \sup_{q \in \{x,y,z\}, j \in \{0,...,n-1\}} |P_{\varepsilon,j}^q| \right)^2} \\ &\ge \sqrt{\sum_{i \in \mathcal{V}} \bar{r}_i^2}, \end{aligned} \quad (44) $$

by evoking (39). Regarding the derivative of the Lyapunov function, it is equivalent to the following potential function when $\bar{r}$ is unvaried:

$$ \frac{\partial E}{\partial \mathbf{p}} g(\mathbf{p}, \bar{\mathbf{r}}) = -\mathbf{f}(\mathbf{p})^{\top} \mathbf{p}. \quad (45) $$

Since the three dynamics considered are DRDs, DSDs, and DPDs, they can all be represented in Laplacian form as presented in (6), (9) and (12).

Without losing generality, we can denote the Laplacian of these dynamics as $L^{(\mathbf{p})}$. This fact leads to

$$ \begin{aligned} \frac{\partial E}{\partial \mathbf{p}} g(\mathbf{p}, \bar{\mathbf{r}}) &= -\mathbf{f}(\mathbf{p})^{\top} L^{(\mathbf{p})} \mathbf{f}(\mathbf{p}) \\ &= -\frac{1}{2} \sum_{i,j=1}^{n} a_{ij}^{(\mathbf{p})} (f_i - f_j)^2 \\ &= -W_3(\mathbf{p}). \end{aligned} \quad (46) $$

As the entries $a_{ij}^{(\mathbf{p})} \ge 0 \forall i, j \in \mathcal{V}$ for all the dynamics as soon as the graph is connected, the second condition is thus proved. $\square$

**Remark 3.** Proposition 2 and Lemma 1 are quite meaningful for many practical scenarios where the formation can be varied. The idea of using B-spline basis functions and control points to parameterize the fitness function is the key point here since the ISS of the population dynamics depends a lot on the trajectories constraints. B-splines provide a convenient way to formalize these constraints, enabling their inclusion in an optimization problem solved at an offline stage. The optimization problem in (37) can include additional constraints related to either the formation collision avoidance or the reconfiguration timing in order to improve the flight performance and safety.

### 3.3. Tracking control

Various stabilizing position controllers can be integrated with the B-spline-based evolutionary dynamics design during the online phase. These controllers can range from simple approaches such as Proportional-Integral-Derivative (PID) and Linear-Quadratic Regulator (LQR) to more advanced methods like MPC (Hu et al., 2020), Pan et al. (2021), flatness-based MPC (Do & Prodan, 2025) or Control Lyapunov Function (CLF) (Do et al., 2025). In this work, we implement a tracking controller for formation control while incorporating a Control Barrier Function (CBF) to ensure safe trajectory tracking and to enforce state constraints, as shown in Fig. 1.

We recall the control barrier function used in the implementation from our previous work in Marguet et al. (2024b). For each agent $i$, after computing the nominal control law from feedback $\xi_i$ and reference

position generated from population dynamics $EV_i$, a safety Control Barrier Function - Quadratic Programming (CBF-QP) filter

$$ \mathbf{u}_i \triangleq \mathbf{k}(u_{des,i}, \xi_i) = \Psi^{-1}(\mu_i) \quad (47) $$

is used to accomplish the safety constraints, i.e., the tracking errors on every axis are bounded as follows:

$$ \min_{\mu_i} \|\Psi(\mathbf{u}_{des,i}(\mathbf{z}_i)) - \mu_i\|_2^2 \quad (48) $$

s.t.

$$ |\mu_{q,i} - \dot{q}_i^{\text{ref}} - a_1(\dot{q}_i^{\text{ref}} - \dot{q}_i) - a_2(q_i^{\text{ref}} - q_i)| \le a_2\delta, \quad (49) $$

where $q \in \mathcal{Q} = \{x, y, z\}$, $\mu$ denotes the virtual control input of system (21), and $\Psi(\mathbf{u}_{des}(\mathbf{z}))$ is the desired virtual stabilizing controller, with $\mathbf{u}_{des}(\mathbf{z})$ the nominal control. $a_1, a_2 > 0$ are chosen such that roots of $\lambda^2 + a_1\lambda + a_2$ are negative reals, according to Freire and Xu (2023). $\delta$ is the safety margin, defining how far the UAV can deviate from the reference.

Algorithm (1) summarizes the online implementation of the proposed distributed scheme.

### Algorithm 1 Distributed formation control.

**Input:** $i$-th agent state, $\xi_i$ in (21), leader's trajectory, $c_i$ in (33), desired gaps, $r_i$ in (36) $\forall i \in \mathcal{V}$

**Output:** Control input with safety assurance, $\mathbf{u}_i$ in (47), actual motion, $z_i$ in (17)

1. **for all** $i \in \mathcal{V}$, with $\mathcal{E}(t)$ in (1) **do**

2. Consider the population dynamics (38) for DRDs in (6), DSDs in (9), or DPDs in (12), with fitness function $f$ in (30)

3. **if** $i \in \mathcal{F}$ **then**

4. $\mathbf{p}_i^{\text{ref}}(t) \leftarrow \mathbf{p}_i(t)$

5. **else**

6. $\mathbf{p}_i^{\text{ref}}(t) \leftarrow c_i(t)$

7. **end if**

8. Compute nominal control from feedback $u_{des,i}(\xi_i, \mathbf{p}_i^{\text{ref}})$

9. Apply CBF-QP safety filter defined in (48) $\mathbf{u}_i \leftarrow \mathbf{k}(u_{des,i})$

10. Apply control to system in (18)

11. **end for**

### 4. Validation

To validate the proposed framework, we assess its effectiveness in ensuring the online convergence of population dynamics, considering the generated B-spline trajectories and the controller's tracking performance. Three scenarios are considered:

**Scenario 1:** We compare the convergence behavior of different population dynamics, DRDs, DSDs and DPDs in (6), (9), and (12) for the formation control task. The goal is to evaluate the framework's scalability and real-time applicability for multiple agents.

**Scenario 2:** Compare various ways for switching the formation, including target points, interpolated line and B-spline references for the agents. The goal is to analyze the benefits of using B-splines to facilitate smooth switching operations.

**Scenario 3:** Test a plug-and-play formation mission where an agent integrates an existing configuration. The goal is to validate the distributed population dynamics for formation convergence, investigating the ease of incorporating additional agents.

The proposed scenarios are first evaluated through simulation. The final two scenarios are additionally validated in our experimental platform. Simulations are performed on a PC equipped with an Intel(R) Core(TM) i7-13700H processor using MATLAB 2023a. Experimental tests are conducted using Crazyflie 2.1 nanodrones, operated in swarm configuration via a radio communication dongle.
<table><caption>**Table 1**
Specifications on the validation scenarios.</caption><thead><tr><th>Scenario</th><th>1</th><th>2</th><th>3</th></tr></thead><tbody><tr><td>Population dynamics</td><td>DRDs, DSDs, and DPDs in (6), (9) and (12)</td><td>DRDs</td><td>DRDs</td></tr><tr><td>No. of drones</td><td>{5, 10, 20, 50}</td><td>4</td><td>4</td></tr><tr><td>No. of control points (13)</td><td>n = 8</td><td>n = 11</td><td>n = 20</td></tr><tr><td>Formation</td><td>Convergence to a circle configuration</td><td>Varying configuration, from square to line</td><td>Varying configuration, from triangle to square</td></tr><tr><td>Simulation / Experiment</td><td>Only simulation</td><td>Both</td><td>Both</td></tr><tr><td>Solver CBFQP</td><td>Matlab</td><td>OSQP (Stellato et al., 2020)</td><td>Matlab and Python</td></tr><tr><td>(a<sub>1</sub>, a<sub>2</sub>, δ) in (48), (49)</td><td></td><td>(6.0, 8.0, 0.1)</td><td></td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Multi-UAV tracking (Scenario 1).</figcaption></figure>

Specifications on the controller parameters and formation configurations for simulations and experiments are given in Table 1.

## 4.1. Simulation scenarios

### Scenario 1. Multi-drone tracking formation

This scenario is designed to assess the real-time applicability of the proposed formation control framework in a multi-drone setting. The leader follows a sequence of predefined waypoints, while the followers track its trajectory to maintain a specified formation. To introduce variability, the initial positions of all agents are perturbed probabilistically by adding zero-mean Gaussian noise with mean $\mu = 0_{3 \times 1}$, and covariance $\sigma = 0.5I_3$. The scenario is evaluated with varying team sizes of 5, 10, 20, and 50 agents to examine the scalability and convergence properties of the distributed evolutionary dynamics under a time-varying fitness function generated using B-spline parameterization.

A circular formation is used to evaluate the performance of the distributed control framework. An illustrative scenario involving 10 agents is presented in Fig. 2, where the formation transitions through the waypoints [6.0, 4.0, 0.7], [9.0, 3.0, 0.7], and ultimately reaches [12.0, 5.0, 0.7]. The trajectories of the follower agents are depicted in black, while the leader's trajectory is shown in green.

The convergence properties of the three distributed dynamics DRDs, DPDs, and DSDs are illustrated in Fig. 3. As shown in Fig. 3(d), DRDs achieve the fastest convergence and lowest computation time among the three methods. In contrast, DPDs and DSDs exhibit slower convergence and higher computational demands. These results highlight the practical trade-offs in performance and efficiency between the different distributed dynamics.

### Scenario 2. Formation switching for four drones

In order to prevent collisions during flight, switching operations for formation control often require either an online motion planning system or safe offline trajectory generation. Moreover, every vehicle must adhere to input limits in addition to maintaining state safety. Our methodology establishes the switching gaps for followers in the formation and the leader's a priori B-spline reference, ensuring constraint compliance and collision avoidance.

In this scenario, the four drones switch between a square and a line formation. We will compare the advantages of B-spline trajectory generation with abrupt target formation changes and linear path switching. After comprehensively comparing the three cases in simulation, the next section will validate the B-spline approach in our experimental platform. This comparative simulation investigates the differences in formation tracking performance between three planning approaches in the offline stage. The reference of the leader is illustrated in Fig. 4 and obtained with the algorithm presented in Marguet et al. (2024a).

Consider a scenario with four UAVs initially arranged in a square formation, where $g_0$ denotes the length of each side. The formation then transitions to a line configuration, with $g_0$ also representing the spacing between consecutive UAVs, before returning to the original square formation prior to landing. In this scenario, the altitude of each follower is assumed to match that of the leader at every time step, implying that all followers maintain a zero gap along the z-axis. The time evolution of the followers' gaps along the x- and y-axes is shown in Fig. 5.

One can notice that the gap is kept constant in the intervals of time [0 8] s, [16.4 21.2] s and [29.6 44] s as those intervals correspond to the time when the formation is supposed to be kept. The two remaining
<figure><img src="image_5.png" /><figcaption>Fig. 3. Convergence of (a) DRDs, (b) DPDs, (c) DSDs over the sampling time $h = 0.1s$. (d) Comparison between DSDs, DPDs and DRDs (Scenario 1).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. B-spline reference trajectory of the leader (Scenario 2).</figcaption></figure>

intervals [8 16.4] s and [21.2 29.6] s corresponds to the transition phases modeled with B-splines.

The degree was set to five and the optimization problem (36) was solved using 11 control points and a uniform clamped knot vector starting at $t_0 = 8$ s and ending at $t_2 = 16.4$ s. The intermediate way-point was associated to the mean of the time interval $t_1 = (t_0+t_2)/2 = 12.2$ s. This way-point is not mandatory but it helps to design trajectories that are not too close from one agent to the other. The way-points for the gaps of all the followers are represented in Table 2 for the transition from the square formation to the line formation. For the other transition (from the line formation to the square formation), we only needed to inverse the way-points (start from $k=2$ until $k=0$) with the associated times $t_0 = 21.2$ s, $t_1 = 25.4$ s and $t_2 = 29.6$ s.

The rigid body associated with the leader is called “DroneE9”, while the 3 followers have a rigid body called, respectively, “DroneE1”,

<figure><img src="image_7.png" /><figcaption>Fig. 5. B-spline gaps on axes x and y of the followers (Scenario 2).</figcaption></figure>

Table 2

Numerical values of the way-points for the gaps (36d) of the three followers.

<table><thead><tr><th>UAV</th><th colspan="3">DroneE1</th><th colspan="3">DroneE2</th><th colspan="3">DroneE3</th></tr><tr><th>k</th><th>0</th><th>1</th><th>2</th><th>0</th><th>1</th><th>2</th><th>0</th><th>1</th><th>2</th></tr></thead><tbody><tr><th>p<sub>k</sub><sup>wp</sup></th><td>g<sub>0</sub></td><td>g<sub>0</sub></td><td>g<sub>0</sub></td><td>0</td><td>g<sub>0</sub></td><td>2g<sub>0</sub></td><td>g<sub>0</sub></td><td>2g<sub>0</sub></td><td>3g<sub>0</sub></td></tr><tr><th></th><td>0</td><td>0</td><td>0</td><td>-g<sub>0</sub></td><td>-g<sub>0</sub></td><td>0</td><td>-g<sub>0</sub></td><td>-g<sub>0</sub></td><td>0</td></tr><tr><th></th><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr></tbody></table>

“DroneE2”, and “DroneE3”. Taking the sum of the reference of the leader and the predefined gap of a given UAV, its expected reference trajectory is known as illustrated in Fig. 6. The reference trajectories are only expected but may differ online if at any time step one or more agents lose the (direct or indirect) communication with the leader or if the replicator dynamics do not have time to converge. The desired
<figure><img src="image_8.png" /><figcaption>Fig. 6. Expected reference trajectory for the followers (Scenario 2).</figcaption></figure>

communication links to be preserved during the whole experiment are between the leader “DroneE9” and the first follower “DroneE1”, between the first follower “DroneE1” and the second follower “DroneE2” and between the second follower “DroneE2” and the third follower “DroneE3”. In fact, maintaining these links guarantees that any UAV can communicate directly or indirectly with the leader to obtain its next reference point, as determined by the evolutionary dynamics and the adjacency matrix. The communication range was set to $ρ = 1.1$ m, the square formation has a side of length 0.5 m and the line formation is obtained with each follower having an offset of 0.5 m on the x-axis compared to the previous agent. These values allow to have a connected graph for each formation, but also to have a graph that varies as in the line formation, the leader cannot communicate directly with the last follower, as their distance is supposed to be 1.5 m, which is bigger than the communication range $ρ$.

The generated B-spline gap references for followers will be utilized in the distributed population game and tracking control during the online stage. The objective is to analyze the tracking errors in conjunction with the satisfaction of state and input constraints throughout the trajectory.

Figs. 7–9 compare our spline-based reference with abrupt switching and linear interpolation. Fig. 7 shows the tracking formations for the three planning approaches. The abrupt switch method shows large tracking errors due to sudden changes. Linear interpolation smooths the tracking but cannot guarantee continuity between switches. Our B-spline approach achieves smoother trajectories and smaller tracking errors, enabling accurate tracking even in multiple switching scenarios.

Fig. 8 shows the acceleration profiles of the drones. The abrupt switching requires aggressive, high-magnitude control actions, sometimes reaching 1 m/s². Linear interpolation reduces this but still causes peaks at switch points. Our spline-based method keeps actuation below 0.2 m/s², resulting in smoother and safer control.

Fig. 9 shows the inter-agent distances during formation flight. The abrupt switching can cause agents to get too close. Both interpolation and spline approaches better manage inter-agent distances during transitions, with the spline method providing the smoothest and safest spacing.

### Scenario 3. Plug-and-play formation

This scenario explores the advantages of distributed population dynamics in our formation control framework. A group of three agents, including a leader, maintains a triangular formation. The leader's trajectory is generated using a B-spline to ensure it passes through predefined way-points. An external agent, initially outside the group's communication range, follows a pre-defined trajectory. Once it enters the communication region, it establishes connectivity with the group and integrates

into the formation as the fourth vertex, transitioning the structure into a square.

Fig. 10 illustrates the generated B-spline trajectories and the timing of way-points for both the leader agent and the external agent. The set of way-points chosen for leader is denoted as $p_{\ell}^{wp}$ while that for the external agent is $p_{ex}^{wp}$.

The total simulation time is 30 s, with way-points occurring at $t_1 = 0.0$ s, $t_2 = 7.5$ s, $t_3 = 15$ s, $t_4 = 22.5$ s, and $t_5 = 30$ s. Followers 1 and 2 maintain a constant distance of $d_0 = 0.6$ m from the leader, forming a triangular formation, as shown in the tracking simulation in Fig. 11. At $t_2$, the external agent enters the communication range $ρ_0 = 0.8$ m, allowing it to join the connectivity graph with the two followers. From that point onward, it changes from following its own trajectory to joining a square formation with the other two followers. Fig. 12 shows the inter-agent distances for the four agents. The formation objective of the incoming agent is easily and effectively achieved through the distributed population dynamics, which uses only local information about neighborhood agents' references. Moreover, the integration of the new agent into the connectivity graph is accomplished seamlessly, without affecting the behavior of the existing agents.

## 4.2. Experiments

The link for the experimental video is GTBspline. The experiments are conducted in our localization platform with four nanodrones (see Fig. 13). Accurate tracking of agent positions during the physical experiments was achieved using a system of 12 motion capture cameras covering a calibrated tracking volume of approximately $4 \times 4 \times 2.9$ m³. Passive reflective markers were affixed to each drone in uniquely asymmetric configurations to ensure robust identification and pose estimation. The system provides high-precision 6-DoF (Degrees of Freedom) measurements with sub-millimeter positional accuracy and operates at a sampling frequency of 100 Hz. These measurements were used both for real-time distance estimation required in the control loop and for validating the recorded trajectories.

### Experiment 1. Formation switching scenario for four drones, Scenario 2

In our experimental platform, ground effects can be observed during takeoff and landing. To avoid the complexity of modeling these notoriously nonlinear interactions and to ensure conservative safety margins, the trajectories generated using our approach are decomposed as follows:

* the first eight seconds are associated with the vertical take-off phase, reaching an altitude of 80 cm. During this phase, the square formation is supposed to be kept;

* the leader follows its planned trajectory and the followers track the leader while switching from the square formation to the line formation;

* the line formation is kept from 16.4 s to 21.2 s while the agents are moving;

* the formation switches back to the square formation;

* the square formation is kept and the UAVs are coming back to their final position which corresponds to their position at the end of the take-off phase; and

* the last six seconds correspond to the landing phase.

The replicator dynamics are applied during all these phases.

We remind that the main objective of the experiment is the control of the formation by giving the appropriate reference point to each UAV, which means that the contributions are focusing more on motion planning than the individual tracking. The controller chosen for each UAV, is the same LQR (Linear Quadratic Regulator) controller with CBF (Control Barrier Functions) used in Marguet et al. (2024a).

Figs. 14 to 22 allow to analyze the experimental results. In particular, Figs. 14 and 15 confirm that the actual motion of the UAVs follows the
<figure><img src="image_9.png" /><figcaption>Fig. 7. Tracking formation for target points, interpolated line, B-spline references (Scenario 2).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Acceleration of drones for target points, interpolated line, B-spline references (Scenario 2).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9. Inter-distance among drones for target points, interpolated line, B-spline references (Scenario 2).</figcaption></figure>

expected trajectories and that the leader is indeed passing in the neighborhood of the way-points. Fig. 16 confirms that the UAVs are flying inside the allowed box. The magnitude of the velocity is always lower than the upper bound, in Fig. 17. The bounds imposed on thrust, roll and pitch angles are respected in Fig. 18. Fig. 19 shows that the angular velocities are satisfied during the whole experiment for all the UAVs excepting a spike observed in the take-off phase of the first follower.

Fig. 20 confirms that the desired communication links were maintained throughout the experiment, as the distance between UAVs associated with these links is always smaller than $\rho$. All phases of the experiment are represented with colored boxes in this figure. Fig. 21 augments the results of the previous figure by representing all the inter-UAV distances during the experiment. This fact allows first to show that there was no collision as the minimal measured value is 0.32 m. We remind the reader that the half size of the nanodrones used is 0.05 m, meaning that if the center of two nanodrones has a distance higher than 0.1 m, then there is no collision between those two nanodrones. The figure also highlights the fact that the communication link between the

leader “DroneE9” and the third follower “DroneE3” is indeed broken from 12.4 s to 25.6 s (around the time associated to the line formation) as expected. The consequence is that the communication graph was varying as depicted in the top of the figure, which leads to a varying adjacency matrix in the replicator dynamics. The loss of this communication link did not prevent the reconfiguration in the square formation as the third follower succeeded to communicate with the other followers that were actually receiving the information from the leader. One may observe that when the formation is supposed to be kept (from 16.4 to 21.2 s for the line formation and from 29.6 to 38 s for the square formation), the inter-UAV distances are almost kept constant, confirming that the formation was kept during the associated time interval.

Fig. 22 shows the tracking performance of the controller used: each color represents a UAV and each subplot corresponds to a direction (upper or lower) on each axis. Having all curves of the same color positive means that the actual position of the UAV is staying during the whole experiment inside the corridor centered in the desired position of the given agent with a size of $\delta = 0.1$ m for each axis. It can be seen that,
<figure><img src="image_12.png" /><figcaption>Fig. 10. Reference trajectories of leader and external agents (Scenario 3).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. Simulation: Tracking formation for plug-and-play scenario (Scenario 3).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12. Simulation: Inter-distance between agents (Scenario 3).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 13. Experimental platform for testing.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14. Experiment: Actual motion of the agents (Scenario 2).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 15. Top view of Fig. 14.</figcaption></figure>

except the phases of take-off and landing that are not the point of this paper, all curves are indeed positive.

To conclude on this experimental section, our approach imposing many constraints for the leader also succeeded to satisfy these constraints for the followers without imposing them directly. The replicator dynamics combined with smooth B-spline gaps succeeded to satisfy the
<figure><img src="image_18.png" /><figcaption>Fig. 16. Experiment: Position on each axis of the agents (Scenario 2).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 17. Experiment: Magnitude of the velocity of each agent (Scenario 2).</figcaption></figure>

objective of reconfiguring the formation while respecting the dynamical constraints of all the agents.

## Experiment 2. Plug-and-play formation, Scenario (3)

Similar to the **Experiment 1**, we design taking-off and landing phases for this experiment. The first $t = 7$ s, the four drones take off to a height of 0.8 m and then start following the desired formation online (depicted in Fig. 11 using the distributed replicator dynamics with the designed fitness function.

The rigid bodies in this experiment are designated “DroneE1” for the external agent, “DroneE9” for the leader, and “DroneE2” and “DroneE3” for the two followers, respectively.

Between $t = 7$ s and $t = 14.5$ s, the three initial agents maintained a triangular formation, each separated by a fixed distance of 0.6m, as illustrated in Figs. 23 and 25. The communication range for this experiment was set to $ρ = 0.8$ m, allowing direct communication between adjacent agents while preventing connectivity with the agent positioned diagonally. Meanwhile, an external agent started from a distant position and gradually approached the formation. At $t = 14.5$ s, this agent came within range of the two followers and was able to establish communication links with them.

After integration, the full group of four agents continued to move in square formation, with the leader following its precomputed trajectory generated by a B-spline. The others adjusted accordingly to maintain the geometric configuration. At the end of the experiment, from $t = 37$ s

<figure><img src="image_20.png" /><figcaption>Fig. 18. Experiment: Thrust, roll angle and pitch angle of each agent (Scenario 2).</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 19. Experiment: Angular velocities of each agent (Scenario 2).</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 20. Experiment: Inter-UAV communication to keep (Scenario 2).</figcaption></figure>

to $t = 44$ s, the four agents performed a landing from their current formation.

As shown in Fig. 24, all agents successfully executed the formation control task. The results confirm that the proposed online distributed B-spline-based population dynamics approach is robust and flexible, enabling the system to incorporate new agents on the fly without degrading performance.
<figure><img src="image_23.png" /><figcaption>Fig. 21. Experiment: Communication graph evolving in function of the inter-UAV distances (Scenario 2).</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 22. Experiment: Tracking performance on each axis for each UAV (Scenario 2).</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 23. Experiment: Actual motion of the agents for plug-and-play formation (Scenario 3).</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 24. Experiment: Position on each axis of the agents for plug-and-play formation (Scenario 3).</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig. 25. Experiment: Inter-agent distance for plug-and-play formation (Scenario 3).</figcaption></figure>

## 5. Conclusion

This paper presents a new method to reconfigure the formation of UAVs taking into consideration online the communication links between agents. The trajectory of the leader is computed offline as a B-spline curve in order to be able to take into account plenty of constraints in a smooth way. Some evolutionary dynamics are then used in an online manner to update the target position of the followers according to the current adjacency matrix and the predefined B-spline gaps representing the expected formation variation. Experiments confirm the expected behavior of the formation and highlight the fact that the smoothness of the B-spline allow respecting several constraints for the followers without imposing them in an optimization problem. The presented method assumed that the pattern for each agent was known a priori, meaning that each agent had a predefined position in all the formations. This assumption allowed to design the transition phases such that the UAVs do not collide. We have shown that the work presented in this paper outperforms already reported approaches using evolutionary dynamics for the distributed formation control.

Future work will concern online methods to find the optimal pattern for each agent to achieve the desired formations, which will imply to consider as well the collision avoidance during the transition phases.

## CRediT authorship contribution statement

**Cong Khanh Dinh:** Writing – review & editing, Writing – original draft, Validation, Software, Investigation, Formal analysis; **Vincent Marguet:** Writing – review & editing, Writing – original draft, Vali-
dation, Software, Formal analysis; Ionela Prodan: Writing – review & editing, Writing – original draft, Supervision, Resources, Methodology, Funding acquisition, Formal analysis, Conceptualization; Carlos Ocampo-Martinez: Writing – review & editing, Writing – original draft, Visualization, Formal analysis; Nicanor Quijano: Writing – review & editing, Writing – original draft, Visualization, Formal analysis; Julian Barreiro-Gomez: Writing – review & editing, Writing – original draft, Methodology, Conceptualization.

## Data availability

Data will be made available on request.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgements

This work has benefited from a French government grant managed by the Agence Nationale de la Recherche under the France 2030 program, reference ANR-23-IACL-0006; the La Région Auvergne-Rhône-Alpes, Pack Ambition Recherche 2021 - PlanMAV, reference REPLAMALCIR; International Research Booster Univ. Grenoble Alpes, reference C900IDX24V30. The work of J. Barreiro-Gomez is supported by the Khalifa University Faculty Start-Up (FSU) Grant 8474000774. The work of C. Ocampo-Martinez is supported by the Spanish project SEAMLESS: “Sustainable learning-based Management of Multi-resource Large-scale Systems” (ref. PID2023-1488400OB-100), funded by MCIN/AEI/10.13039/501100011033/FEDER, UE. J. Barreiro-Gomez is profoundly grateful to God and to the intercession of Our Lady of Lourdes for the blessings of health and life, without which his contributions to this work would not have been possible.

## Appendix A. The choice of m in (39)

Let m denote the mass of the population, as defined in the context of the simplex. In the population dynamics under consideration, m can be selected such that the state vector remains within the simplex, i.e., all components are strictly positive and their sum is equal to m.

To establish input-to-state stability (ISS), we utilize the simplex structure to prove the inequality

$$ \|p\| \ge \|\bar{r}\| = \sqrt{\sum_{i=1}^{v} \bar{r}_i^2}, \qquad (A.1) $$

where $\bar{r}_i$ represents the deviation components of the generated trajectories. In Eq.(40), this is equivalent to showing:

$$ m^2 \ge \sum_{i=1}^{v} \bar{r}_i^2. \qquad (A.2) $$

Note that the selection of this m satisfies the condition in (27) provided that $v \sum_{i \in F} r_i > \sum_{i \in F} r_i$. For a system with v agents, consisting of v-1 followers and one leader, the summation $\sum_{i=1}^{v} \bar{r}_i^2$ can be partitioned as follows:

### • Followers:

$$ \sum_{i=1}^{v-1} \bar{r}_i^2 \le (v-1) \cdot \max_i(\bar{r}_i^2) \le (v-1) \cdot \max (\|\bar{r}\|, \|\bar{r}\|), \qquad (A.3) $$

where $\bar{r}$ and $\bar{r}$ denote the lower and upper bounding trajectories for the followers.

### • Leader:

$$ \bar{r}_i^2 = c_i^2 \le \sup (P_i^{(q)})^2, \qquad (A.4) $$

where $c_i$ is bounded by the convexity property of the B-spline trajectory. As $c_i$ corresponds to a cost-related control signal in the game-theoretic formulation and is assumed positive, its upper bound is given by the largest control point of the B-spline.

Therefore, by selecting m as in Eq. (43), we ensure that the inequality in (40) is satisfied, thereby fulfilling the ISS condition.

## References

* Barreiro-Gomez, J., Mas, I., Giribet, J. I., Moreno, P., Ocampo-Martinez, C., Sánchez-Peña, R., & Quijano, N. (2021). Distributed data-driven UAV formation control via evolutionary games: Experimental results. *Journal of the Franklin Institute*, 358 (10), 5334–5352.

* Barreiro-Gomez, J., Mas, I., Ocampo-Martinez, C., Sanchez-Peña, R., & Quijano, N. (2016a). Distributed formation control of multiple unmanned aerial vehicles over time-varying graphs using population games. In *2016 IEEE 55th conference on decision and control (CDC)* (pp. 5245–5250). IEEE.

* Barreiro-Gomez, J., Obando, G., & Quijano, N. (2016b). Distributed population dynamics: Optimization and control applications. *IEEE Transactions on Systems, Man, and Cybernetics: Systems*, 47(2), 304–314.

* Barreiro-Gomez, J., & Park, S. (2025). Optimal strategy revision in population games: A mean field game theory perspective. *arXiv:2501.01389*.

* Barreiro-Gomez, J., & Poveda, J. I. (2023). Recurrent neural network ODE output for classification problems follows the replicator dynamics. *IEEE Control Systems Letters*, 7, 3783–3788.

* Barreiro-Gomez, J., & Tembine, H. (2018). Constrained evolutionary games by using a mixture of imitation dynamics. *Automatica*, 97, 254–262.

* Cao, N., & Lynch, A. (2021). Predictor-based control design for UAVs: Robust stability analysis and experimental results. *International Journal of Control*, 94(6), 1529–1543.

* Do, H.-T., Blanchini, F., Miani, S., & Prodan, I. (2025). Reducing real-time complexity via sub-control Lyapunov functions: from theory to experiments. *arXiv:2501.09143*.

* Do, H.-T., & Prodan, I. (2025). Explicit flatness-based MPC design for indoor nanoquadcopter position tracking. *International Journal of Control*, 98(2), 402–413.

* Do, H.-T., Prodan, I., & Stoican, F. (2021). Analysis of alternative flat representations of a UAV for trajectory generation and tracking. In *25th International conference on system theory, control and computing* (pp. 58–63). IEEE.

* Freire, V., & Xu, X. (2023). Flatness-based quadcopter trajectory planning and tracking with continuous-time safety guarantees. *IEEE Transaction on Control Systems Technology*, 31(6), 2319–2334.

* Hu, J., Wang, M., Zhao, C., Pan, Q., & Du, C. (2020). Formation control and collision avoidance for multi-UAV systems based on Voronoi partition. *Science China Technological Sciences*, 63(1), 65–72.

* Khalil, H. K., & Grizzle, J. W. (2002). Nonlinear systems (vol. 3). Prentice hall Upper Saddle River, NJ.

* Lyche, T., Manni, C., & Speleers, H. (2018). Foundations of spline theory: B-splines, spline approximation, and hierarchical refinement. In *Splines and PDEs: From approximation theory to numerical linear algebra* (pp. 1–76). Springer.

* Marguet, V., Dinh, C. K., Prodan, I., & Stoican, F. (2024a). Constrained PSO-splines trajectory generation for an indoor nanodrone. In *2024 International conference on unmanned aircraft systems, IICUAS’24*.

* Marguet, V., Dinh, C. K., Stoican, F., & Prodan, I. (2024b). Indoor formation motion planning using b-splines parametrization and evolutionary optimization. *Control Engineering Practice*, 152, 106066.

* Martinez-Piazuelo, J., Quijano, N., & Ocampo-Martinez, C. (2022a). Nash equilibrium seeking in full-potential population games under capacity and migration constraints. *Automatica*, 141, 110285.

* Martinez-Piazuelo, J., Quijano, N., & Ocampo-Martinez, C. (2022b). A payoff dynamics model for generalized nash equilibrium seeking in population games. *Automatica*, 140, 110227.

* Miele, A., Lippi, M., & Gasparri, A. (2024). A framework based on control barrier functions for time-varying connectivity maintenance. In *2024 IEEE 20th international conference on automation science and engineering (CASE)* (pp. 84–89). IEEE.

* Nguyen, N. T., Schilling, L., Angern, M. S., Hamann, H., Ernst, F., & Schildbach, G. (2021). B-spline path planner for safe navigation of mobile robots. In *2021 IEEE/RSJ International conference on intelligent robots and systems (IROS)* (pp. 339–345). IEEE.

* Obando, G., Pantoja, A., & Quijano, N. (2013). Building temperature control based on population dynamics. *IEEE Transactions on Control Systems Technology*, 22(1), 404–412.

* Obando, G., Quijano, N., & Ocampo-Martinez, C. (2022). Decentralized control for urban drainage systems using replicator dynamics. *IEEE access*, 10, 56740–56762.

* Pan, Z., Sun, Z., Deng, H., & Li, D. (2021). A multilayer graph for multiagent formation and trajectory tracking control based on MPC algorithm. *IEEE Transactions on Cybernetics*, 52(12), 13586–13597.

* Park, S., & Leonard, N. E. (2025). Learning with delayed payoffs in population games using Kullback-Leibler divergence regularization. *IEEE Transactions on Automatic Control*, (pp. 1–16). DOI: 10.1109/TAC.2025.3561108.

* Quijano, N., Ocampo-Martinez, C., Barreiro-Gomez, J., Obando, G., Pantoja, A., & Mojica-Nava, E. (2017). The role of population games and evolutionary dynamics in distributed control systems: The advantages of evolutionary game theory. *IEEE Control Systems Magazine*, 37(1), 70–97.

* Sandholm, W. H. (2010). Population games and evolutionary dynamics. MIT press.
Sforni, L., Carnevale, G., & Notarstefano, G. (2024). A distributed feedback-based framework for nonlinear aggregative optimal control. *IEEE Transactions on Automatic Control*.

Sontag, E. D., & Wang, Y. (1995). On characterizations of the input-to-state stability property. *Systems & Control Letters, 24*(5), 351–359.

Sontag, E. D. (2008). Input to state stability: Basic concepts and results. In Agrachev, A. A., Morse, A. S., Sontag, E. D., Sussmann, H. J., & Utkin, V. I., *Nonlinear and optimal control theory: lectures given at the CIME summer school held in Cetraro, Italy June 19–29, 2004*, (pp. 163–220). Springer.

Stellato, B., Banjac, G., Goulart, P., Bemporad, A., Boyd, S.(2020). OSQP: An operator splitting solver for quadratic programs. *Mathematical Programming Computation, 12*(4), 637–672.

Vargas, S., Becerra, H. M., & Hayet, J.-B. (2022). MPC-based distributed formation control of multiple quadcopters with obstacle avoidance and connectivity maintenance. *Control Engineering Practice, 121*, 105054.

Vidal, J. M. (2001). Learning in multiagent systems: An introduction from a game-theoretic perspective. In *Symposium on adaptive agents and multi-agent systems* (pp. 202–215). Springer.

Zhang, Z., Wu, Y., Jiang, J., Zheng, N., & Meng, W. (2024). Realization of robust formation for multi-UAV systems using control barrier functions. *Unmanned Systems, 13*(03), 1–15.