Delayed Interactions in Active Agents: Stability and Formations

TL;DR

This paper analyzes the stability of pattern formations in active agents with inertia and time delays, deriving conditions for stability and exploring the effects of delays in both uncoupled and coupled systems.

Contribution

It introduces a comprehensive high-dimensional DDE model for analyzing stability of formations with delays, including new criteria for delay-independent and delay-dependent stability.

Findings

Delay-independent stability conditions identified

Spectrum classification for coupled agents provided

Stability region converges to a circle in large-delay limit

Abstract

Active agents with time-delayed interactions arise naturally in various real-world systems, such as biological systems, transportation networks and robotic swarms. Such systems are typically modeled as Delay Differential Equations (DDEs) that incorporate inertial effects. In this paper, we investigate the stability of pattern formation of active agents with inertia and time delays, in both uncoupled and coupled scenarios. We derive and analyze a high-dimensional linear DDE model that characterizes the stability of such formations. Starting with the uncoupled scenario, where agents are driven only by a virtual leader, we describe the stability spectrum and provide conditions for the delay-independent (absolute) stability of the formations, as well as delay-dependent stability and unstable hyperbolic behavior. Different cases correspond to distinct universality classes of the corresponding spectrum. For the coupled scenario, where agents are driven by both the virtual leader and inter-agent interactions, we consider both symmetric and non-symmetric coupling topologies. Here we also provide an explicit spectrum classification, including the absolute stability criterion. Additionally, we investigate interactions in the large-delay limit, where delays affect inter-agent coupling, while local feedback remains instantaneous. In this limit, we prove rigorously that the stability region in the complex plane of the eigenvalues of the Laplacian matrix converges to a circle centered at the origin, a phenomenon previously observed in delay-coupled networks. Our findings provide a universal framework for understanding stable formations and motions of active agents with delayed interactions.

Equations

$$\dot{\mathbf{R}}(t)$$

$$\dot{\mathbf{V}}(t)$$

$$\mathbf{e}$$

$$\mathbf{\xi }$$

$$U_i=$$

$$+k^{\tau }\left[\left(R_i(t-\tau )-R_j(t-\tau )\right)-\left(s_i-s_j\right)\right]$$

$$\displaystyle+h\left[V_{i}(t)-V_{j}(t)\right]+h^{\tau}\left[V_{i}(t-\tau)-V_{j}(t-\tau)\right]\Bigr{\}}$$

$$-k_0\left[R_i(t)-\left(R_0(t)+s_i\right)\right]-k_0^{\tau }\left[R_i(t-\tau )-\left(R_0(t-\tau )+s_i\right)\right]$$

$$-h_0\left[V_i(t)-V_0(t)\right]-h_0^{\tau }\left[V_i(t-\tau )-V_0(t-\tau )\right]+U_0(t),$$

$$i=1,2,\dots ,N.$$

$$U_i=$$

$$\displaystyle\quad\quad\quad\quad\quad\quad+h\left[\xi_{i}(t)-\xi_{j}(t)\right]+h^{\tau}\left[\xi_{i}(t-\tau)-\xi_{j}(t-\tau)\right]\Bigr{\}}$$

$$-k_0{e}_i(t)-k_0^{\tau }{e}_i(t-\tau )-h_0{\xi }_i(t)-h_0^{\tau }{\xi }_i(t-\tau )+U_0(t).$$

$$\mathbf{U}=$$

$$-k_0\mathbf{e}(t)-k_0^{\tau }\mathbf{e}(t-\tau )-h_0\mathbf{\xi }(t)-h_0^{\tau }\mathbf{\xi }(t-\tau )+{\mathbf{U}}_0(t),$$

$$\displaystyle\left[\mathbf{L}\right]_{ij}=\left\{\begin{array}[]{cc}\sum\limits^{N}_{l=1,l\neq i}a_{il},&i=j;\\ -a_{ij},&i\neq j;\quad i=1,2,\cdots N,\end{array}\right.$$

$$\dot{\mathbf{e}}=$$

$$\dot{\mathbf{\xi }}=$$

$$-k_0\mathbf{e}(t)-k_0^{\tau }\mathbf{e}(t-\tau )-h_0\mathbf{\xi }(t)-h_0^{\tau }\mathbf{\xi }(t-\tau ).$$

$$\dot{Z}(t)=$$

$$+\left[\left(1_N\otimes M^{\tau }-\text{L}\otimes P^{\tau }\right)\otimes 1_3\right]Z(t-\tau ),$$

$$M=$$

$$P=$$

$${\dot{Z}}_{ijs}=({\delta }_{im}{M}_{jl}-[\mathbf{L}{]}_{im}{P}_{jl})Z_{mls}(t)+({\delta }_{im}{M}_{jl}^{\tau }-[\mathbf{L}{]}_{im}{P}_{jl}^{\tau })Z_{mls}(t-\tau ).$$

$${\mathbf{H}}^{-1}\mathbf{LH}=\mathbf{\Lambda }=\mathrm{diag}\left\{{\lambda }_1,\cdots ,{\lambda }_N\right\},$$

$$\dot{X}(t)=$$

$$+\left[\left(1_N\otimes M^{\tau }-\mathbf{\Lambda }\otimes P^{\tau }\right)\otimes 1_3\right]X(t-\tau ).$$

$${\dot{x}}_{\ell }(t)=M{x}_{\ell }(t)+M^{\tau }{x}_{\ell }(t-\tau )-{\lambda }_{\ell }\left[P{x}_{\ell }(t)+P^{\tau }{x}_{\ell }(t-\tau )\right],$$

$$\dot{x}(t)=(M-\lambda P)x(t)+(M^{\tau }-\lambda P^{\tau })x(t-\tau ),$$

$$\dot{x}(t)=Ax(t)+Bx(t-\tau ).$$

$$\det \left[\mu I-A-B{e}^{-\mu \tau }\right]=0.$$

$${\chi }_{\omega }(Y):=\det \left[i\omega I-A-BY\right]$$

$${\Lambda }_{\text{ACS}}:=\left\{\frac{1}{\tau }{\gamma }_j(\omega )+i\omega \in \mathbb{C}:\omega \in \mathbb{R},j=1,\dots ,m\right\},$$

$$\text{where}{\gamma }_j(\omega )=-\ln \left|Y_j(\omega )\right|.$$

$${\tau }_j=\frac{1}{{\omega }_H}({\varphi }_H+2\pi j),j=0,1,\dots$$

$${\tau }_j^{(1)}=\frac{1}{{\omega }_{H_1}}({\varphi }_{H_1}+2\pi j),$$

$${\tau }_j^{(2)}=\frac{1}{{\omega }_{H_2}}({\varphi }_{H_2}+2\pi j),$$

$$\mathbf{U}=-k_0\mathbf{e}(t)-k_0^{\tau }\mathbf{e}(t-\tau )-h_0\mathbf{\xi }(t)-h_0^{\tau }\mathbf{\xi }(t-\tau )+{\mathbf{U}}_0(t),$$

Full text

1 Introduction

Active agents possess the capability to move autonomously by perceiving their environment and information about interactions accordingly. They often exhibit collective behaviors when such agents interact with each other [42, 39]. Such systems are prevalent across many fields, ranging from bacterial colonies [41, 39] and bird flocks [31] to robotic swarms [32, 40, 27, 33]. The processes of perception, feedback, and decision-making of information in these systems inherently involve time delays, resulting in delayed active agent systems. Time delays significantly alter system dynamics, including behaviors such as oscillatory or chaotic motion, and complex behavior formation. Active agents with time delays can be found in various applications, cf. Refs. [12, 54, 23, 21, 38, 48, 52, 60, 6] for a recent overview.

In systems of active agents with time-delayed interactions, each agent adjusts its motion based on delayed information about its own state or the states of connected agents. These interactions are commonly modeled using Delay Differential Equations (DDEs) with inertial effects, which emphasize the influence of past states on current dynamics [5, 17, 21, 9, 60, 56, 44, 22]. The inertial delay observed in active agents arises from the finite relaxation time in their motion, which deeply influences long-term dynamics and therefore typically requires consideration [40, 2]. In this study, we focus on a class of models in which the dynamics of agents are governed by a high-dimensional linear system of DDEs, incorporating time delays into both local feedback and interactions. A key property of such systems is that the stability of collective formations can be explained by the stability of the linear high-dimensional DDE system. By analyzing the spectrum of this DDE system and in particular by employing the Asymptotic Continuous Spectrum (ACS) technique [53, 57, 37, 59, 51], we can determine the conditions whether a desired formation will remain stable, persists, or become unstable in the presence of time delays.

For a general linear time-delayed system, absolute stability guarantees stability for all delay values. In Ref. [59], we derived explicit criteria for the absolute stability in linear DDEs. Further development of these results in Ref. [51] provides a universality classification of linear DDEs such that each class has the same sequence of stabilizing or destabilizing bifurcations when time delay is changing. Thus ACS is a useful tool for analyzing the stability and bifurcation behavior of DDEs. Originally introduced for DDEs with large delays [8, 53, 58, 25, 43, 57, 56, 20, 37, 59, 35, 51], the ACS provides a computationally rather basic framework for characterizing desctabilization scenarios.

In this study, we examine thoroughly two cases: a system of uncoupled agents driven only by a virtual leader (prescribed target function), and a coupled system of active agents not only driven by a virtual leader but also interacting with each other via a time delay. We employ the ACS and the Master Stability Function (MSF) [34, 18, 1, 36, 15, 7] approaches to provide theoretical analysis of these two cases. Using the universal classification introduced in [51], we provide explicit conditions for the agent formation to be: (i) stable for all delays (absolutely stable or universality class 0); (ii) to undergo an explicitly determined destabilization sequence as time delay increases (universality class I); (iii) to undergo explicitly determined destabilization and stabilization sequences as time delay increases (universality class II); or (iv) unstable and hyperbolic for all time delays (universality class $U$). In particular, when delayed coupling is introduced, the stability of the formation can be described by a combination of the above cases.

This work is organized as follows. In Section 2, we derive a high-dimensional linear DDE model that describes the dynamics of active agents. The stability of this model directly determines the stability of the resulting formations. Section 3 describes two analytical tools for studying formations of interacting agents: the master stability function and the asymptotic continuous spectrum. Section 4 focuses on the uncoupled case, where agents are driven solely by the virtual leader. We identify four universality classes (0, I, II and U) of ACS of the high-dimensional DDE system and illustrate them with a bifurcation diagram. We then analyze the motion and pattern formation of active agents within each class. Section 5 extends this analysis to the coupled case, where active agents are driven by a virtual leader and interact with others incorporating time delays. We then examine the motion and formation behaviors of the agents under both symmetric and asymmetric interaction topologies. We show how the spectrum of the whole coupled system can be expressed as a combination of the spectral branches of types 0, I, II, and U. Furthermore, we examine the large delay limit, in which delays are present only in the interactions. Notably, in this limit, the stability region in the plane of the complex eigenvalues of the Laplacian (describing the coupling topology) asymptotically approaches a circle centered at the origin. This circular property of the “master stability function” has previously been observed in delay-coupled networks [7], but here we provide a first rigorous proof of this property. Finally, Section 6 summarizes our main results and outlines potential directions for future research.

2 The linear model of active agents

We consider the following model of $N$ active agents [30, 46, 14, 40, 3, 24, 61, 50, 10, 55, 11] with inertia and time-delayed control

[TABLE]

where $\mathbf{R}:\mathbb{R}\to\mathbb{R}^{3N}$ is the position, $\mathbf{R}(t)=\left(R_{1}(t),\dots,R_{N}(t)\right)$, $R_{i}(t)\in\mathbb{R}^{3}$, $\mathbf{V}:\mathbb{R}\to\mathbb{R}^{3N}$ is the velocity, $\mathbf{V}(t)=\left(V_{1}(t),\dots,V_{N}(t)\right)$, $V_{i}(t)\in\mathbb{R}^{3}$ of moving agents. The system moves due to the control input force $\mathbf{U}:\left(\mathbb{R}^{3N}\right)^{4}\times\mathbb{R}\to\mathbb{R}^{3N}$, $\mathbf{U}=\left(U_{1},\dots,U_{N}\right)$, $U_{i}\in\mathbb{R}^{3}$.

The main goal of the control force $\mathbf{U}$ is to achieve a desired motion and formation of agents. More specifically, the target trajectory $R_{i}(t)$ of agent $i$ is prescribed by a time-dependent function $R_{0}(t)+s_{i}$, where $s_{i}\in\mathbb{R}^{3}$ determines the position of agent $i$ in the formation. Hence, the vector $\mathbf{s}=\left(s_{1},\dots,s_{N}\right)\in\mathbb{R}^{3N}$ determines the desired formation shapes of active agents. By denoting $\mathbf{R_{T}}=\left(R_{0}(t),\dots,R_{0}(t)\right)+\mathbf{s}\in\mathbb{R}^{3N}$, the control goal is to find conditions on the control force $\mathbf{U}$ such that $\mathbf{R}(t)\to\mathbf{R_{T}}(t)$ asymptotically and exponentially fast with $t\to\infty$. We assume that $R_{0}(t)$ (and hence $\mathbf{R_{T}}(t)$) is a two times continuously differentiable function.

We define the error variables for the deviations of the positions $\mathbf{e}(t)=\left(e_{1}(t),\dots,e_{N}(t)\right)$ and velocities $\boldsymbol{\xi}(t)=\left(\xi_{1}(t),\dots,\xi_{N}(t)\right)$ from the target functions as

[TABLE]

Here $\otimes$ denotes the Kronecker product and $\mathrm{1}_{N}$ is the $N$-dimensional vector with identical elements equal to 1.

We consider the Proportional Derivative (PD) controller with delayed terms, as studied in Refs. [47, 19, 13, 26, 16, 28, 49]. Adding delay terms to the standard PD controller can improve its response characteristics and ensures better performance. The considered PD controller has the following form:

[TABLE]

The coefficients $a_{ij}$, $j\neq 0$ determine the coupling between the agents $i$ and $j$; $a_{ij}>0$ if the information is transmitted from $j$ to $i$, otherwise, $a_{ij}=0$. $k$, $k^{\tau}$, $h$, $h^{\tau}$ are the corresponding control gains between the agents. $k_{0}$, $k^{\tau}_{0}$, $h_{0}$, and $h^{\tau}_{0}$ determine the coupling between the agents and the virtual leader. $V_{0}(t)=\dot{R}_{0}(t)$ and $U_{0}(t)=\ddot{R}_{0}(t)$ are the velocity and acceleration of the “virtual leader”, i.e. a prescribed position in space that should be followed by the active agents.

Using the error variables (3) and (4), the considered PD controller reads

[TABLE]

In vector form, the controller is expressed as

[TABLE]

where $\mathbf{L}$ is the Laplacian matrix

[TABLE]

and $\mathbf{U}_{0}(t)=1_{N}\otimes U_{0}(t)$.

Finally, using the relations $\mathbf{\dot{e}}=\mathbf{\dot{R}}-\mathbf{\dot{R}_{T}}$ and $\boldsymbol{\dot{\xi}}=\mathbf{\dot{V}}-\mathrm{1}_{N}\otimes\ddot{R}_{0}(t)=\mathbf{U}-\mathbf{U}_{0}(t)$, we obtain the following autonomous linear delay-differential equation describing the error dynamics

[TABLE]

The system defined by equations (11)–(2) governs the stability of the agent formation and the tracking of the virtual leader orbit. This will be the main object of this study.

System (11)–(2) can be rewritten in the following form

[TABLE]

where $Z(t)=[e_{1},\xi_{1}\cdots,e_{N},\xi_{N}]^{T}\in\mathbb{R}^{N\times 2\times 3}$ is the state vector. More precisely, $Z$ is a rank three tensor defined by $Z_{i1j}=e_{ij}$ and $Z_{i2j}=\xi_{ij}$ for $i=1,\dots,N$ and $j=1,2,3$. Matrices $M,~M^{\tau}$, $P,~P^{\tau}\in\mathbb{R}^{2\times 2}$ are given by

[TABLE]

The connectivity structure is given by the Laplacian matrix $\mathbf{L}\in\mathbb{R}^{N\times N}$ with zero row sum. For completeness, equations (13) can be written using explicit indices and the Einstein sum convention as follows

[TABLE]

3 Preliminaries

This section introduces two tools that will be used for studying the pattern formation of interacting agents: the Master Stability Function (MSF) and the Asymptotic Continuous Spectrum (ACS) approach for delay systems.

3.1 Master Stability Function for delay-coupled systems

We assume that the Laplacian matrix $\mathbf{L}$ can be diagonalized such that

[TABLE]

where $\lambda_{\ell}$ are the eigenvalues of $\mathbf{L}$. In this case, system (13) can be block diagonalized, similarly to the MSF approach in Ref. [34]. We define the new variable $X(t)$ using $X(t)=\left(\mathbf{H}\otimes{1_{2}}\otimes 1_{3}\right)^{-1}Z(t)$. Then, Eq. (13) has the block diagonal form with respect to the new variable

[TABLE]

This block structure leads to the set of $N$ independent equations

[TABLE]

where $x_{\ell}(t)\in\mathbb{R}^{2},~\ell=1,2,\dots,N$. Since the set of equations (16) only differs by the parameter $\lambda_{\ell}$, the stability problem is reduced to a single two-component delay equation

[TABLE]

where $x(t)\in\mathbb{R}^{2}$ and the parameter $\lambda$ takes the values $\lambda_{\ell}$ of the eigenvalues of the Laplacian matrix $\mathbf{L}$. System (17) allows a separation of the stability problem and the coupling topology given by $\mathbf{L}$. Solving the stability problem for (17) with $\lambda$ as a parameter provides stability conditions for an arbitrary coupling topology.

3.2 Delay-independent classification of linear DDEs

Using the notations $A=M-\lambda P$ and $B=M^{\tau}-\lambda P^{\tau}$, equation (17) can be rewritten in a more general form as a linear DDE with one delay

[TABLE]

The corresponding characteristic equation determining the stability and eigenvalues is

[TABLE]

To study the stability properties of (18) and the roots of the characteristic equation (19) as a function of time delay $\tau$, we will use the methods developed in [59] and [51]. These methods introduce classes of DDEs, which are either stable for all time delays or undergo different but universal destabilization scenarios. To classify these DDEs, we will introduce the concepts of the generating polynomial and the asymptotic continuous spectrum.

Definition 3.1**.**

[51, Definition 2.3]** The function

[TABLE]

*is called the generating polynomial. $\chi_{\omega}(Y)$ is a polynomial with respect to both $Y\in\mathbb{C}$ and $\omega\in\mathbb{R}$. We also denote the roots of the generating polynomial as $Y_{j}(\omega)$, i.e. $\chi_{\omega}(Y_{j}(\omega))=0$, where $j=1,\dots,m$. $m$ is the rank of matrix $B$. *

Using the roots of the generating polynomial, we introduce the following concept of the asymptotic continuous spectrum.

Definition 3.2**.**

[51*, Definition 2.4]**

The Asymptotic Continuous Spectrum (ACS) is given by*

[TABLE]

The ACS is shown to approximate the spectrum of DDEs of the form (18) for large delays [25, 43]. It consists of $m$ continuous curves in $\mathbb{C}$ in general. In our case, we have $B=M^{\tau}-\lambda P^{\tau}$, $\textnormal{rank}\,B=1$ and there is only one ACS curve.

Furthermore, we introduce the following classification of the ACS, which is closely related and will lead to the appropriate classification of the DDEs [51]:

Class 0 ACS shows no crossing with the imaginary axis, see Fig. 1(a).

Class I ACS has two points $\pm i\omega_{H}$ where it crosses the imaginary axis, see Fig. 1(b).

Class II ACS exhibits four crossing points $\pm i\omega_{H_{1}}$ and $\pm i\omega_{H_{2}}$ as shown in Fig. 1(c).

Additionally, a strongly unstable spectrum $\Lambda_{U}$ is defined as a set of eigenvalues of matrix $A$ with positive real parts [25]. Systems with a nonempty strongly unstable spectrum are always unstable for sufficiently large delays. Moreover, such systems become stable for small delays if the matrix $A+B$ is stable.

Based on the above discussion, we now introduce the concepts of absolute stability and universality classes I, II, and U of DDEs (18) to further describe the stability and various forms of instabilities of this system, see more details and proofs in [59, 51].

Definition 3.3**.**

[59*, Definition 1]** System (18) is absolutely stable if all roots $\mu$ of the characteristic equation (19) possess negative real parts $\Re(\mu)<0$ for all $\tau\geq 0$. *

Accordingly to [59], system (18) is absolutely stable iff the ACS is of class 0 and $A+B$ is stable.

Definition 3.4**.**

[51*, Definition 4.1]** We define the system (18) to be of universality class I if its ACS is of class I. *

The destabilization time delays for DDEs of class I are given by [51]

[TABLE]

with $\phi_{H}=-\arg\left[Y(\omega_{H})\right]\in[0,2\pi]$ and $\omega_{H}>0$, which are independent on $\tau$. Here $Y$ is defined by (20). Hence, there are two stability scenarios of class I DDEs: (i) stability for $0\leq\tau<\tau_{0}$ and instability for $\tau>\tau_{0}$ if $A+B$ is stable; and (ii) instability for all $\tau\geq 0$ if $A+B$ is unstable. Note that in case (ii) the DDE is also unstable for all positive time delays.

Definition 3.5**.**

[51*, Definition 5.1]** We define the system (18) to be of universality class II iff its ACS is of universality class II. *

The class II DDEs possess a sequence of destabilizing transitions [51] at

[TABLE]

and stabilizing transitions at

[TABLE]

as time-delay varies, where $0<\omega_{H_{2}}<\omega_{H_{1}}$, and $\phi_{H_{i}}=-\arg\left[Y(\omega_{H_{i}})\right]$ are delay-independent. Moreover, the regions of instability and stability can alternate (known as stability islands [4, 18, 29]), and for sufficiently large $\tau$ the system is unstable.

Definition 3.6**.**

*System (18) is said to be of universality class $U$, if for all $\tau>0$, the characteristic equation (19) has at least one root $\mu$ with a positive real part $\Re(\mu)>0$ (unstable) and no roots occur with $\Re(\mu)=0$ (hyperbolic). *

According to [59], the criterion for a linear DDE system to be of class $U$ is: ACS of type 0 and unstable $A$. In this case, the system remains unstable and hyperbolic for all positive $\tau$.

4 The case of uncoupled agents

4.1 Characteristic equation, weak and strong spectrum

We first analyze the case where all active agents receive information only from the virtual leader, which we call the case of uncoupled agents, see Fig. 2.

When the coupling is absent, we have $\mathbf{L}=\mathbf{0}$ and the control function (7) reads

[TABLE]

hence system (11)–(2) reduces to the following form

[TABLE]

This results in a set of uncoupled equations, each corresponding to an individual agent $i$, where $i=1,\dots,N$:

[TABLE]

Since the equations for all $i$ are identical, we can drop the index $i$. In the vector form, the error system becomes

[TABLE]

with $M$ and $M^{\tau}$ are defined by (14) and $y=(e_{i},\xi_{i})$. The stability of the error system (32) is governed by the characteristic equation

[TABLE]

The associated ACS curve described by (21) (refer to Sec. 3) reduces to a single expression for $\gamma(\omega)$, given by:

[TABLE]

The strongly unstable spectrum $\Lambda_{U}$ of system (31) is then given as the unstable spectrum of the instantaneous part of (31), which are the roots of the characteristic equation

[TABLE]

with positive real parts.

4.2 Stability diagrams for uncoupled agents

In this section, we obtain stability conditions for system (27)–(28) (equivalently, for (32)), which include delay-independent classifications, as well as diagrams for specific time delays.

The following theorem gives an explicit delay-independent classification. To shorten the notations, we denote the set of parameters of system (27)–(28) as $p_{0}:=(k_{0},h_{0},k_{0}^{\tau},h_{0}^{\tau})$.

Theorem 4.1**.**

The error time-delay system (27)–(28) has the following delay-independent classification:

(i)* The system is absolutely stable if and only if $p_{0}\in S_{0}$, where the parameter set $S_{0}$ is defined as*

[TABLE]

(ii)* The system belongs to universality class I if and only if $p_{0}\in S_{I}$ with*

[TABLE]

(iii)* The system belongs to universality class II if and only if $p_{0}\in S_{II}$ with*

[TABLE]

where $S_{II,\mathrm{part}}$ is defined as the set of all $p_{0}$ for which at least one of the following conditions holds

a) $k_{0}>|k^{\tau}_{0}|$,

b) $-\left[\left(\frac{h^{\tau}_{0}}{2}\right)^{2}+\left(\frac{k^{\tau}_{0}}{h^{\tau}_{0}}\right)^{2}\right]\leq k_{0}<-|k^{\tau}_{0}|$.

(iv)* The system belongs to universality class $U$ if and only if $p_{0}\in S_{U}$, where $S_{U}$ is defined as a the set of all $p_{0}$ for which at least one of the following conditions holds*

a) $k_{0}<-\left[\left(\frac{h^{\tau}_{0}}{2}\right)^{2}+\left(\frac{k^{\tau}_{0}}{h^{\tau}_{0}}\right)^{2}\right]$,

b) $-\left[\left(\frac{h^{\tau}_{0}}{2}\right)^{2}+\left(\frac{k^{\tau}_{0}}{h^{\tau}_{0}}\right)^{2}\right]\leq k_{0}<-|k^{\tau}_{0}|$, and $\left|h_{0}\right|>|h^{-}_{0}|$,

*c) $k_{0}>|k^{\tau}_{0}|$ and $h_{0}<-|h^{-}_{0}|$. *

Proof 4.2**.**

Firstly, we note that the strongly unstable spectrum is not empty (the characteristic equation (37) for the instantaneous part of the system has at least one root with positive real part) if and only if $k_{0}<0$ or $h_{0}<0$.

The absolute stability criterion is that the strongly unstable spectrum is absent (i.e., $k_{0}\geq 0$ and $h_{0}\geq 0$), and the ACS belongs to class 0, i.e, it does not cross the imaginary axis [59]. The following equation determines the points, at which the curve of the ACS (34) crosses the imaginary axis and $\gamma(\omega)=0$:

[TABLE]

With $\nu=\omega^{2}$, equation (41) becomes quadratic in $\nu$, and

[TABLE]

The ACS is of class 0 if this quadratic equation has no real non-negative roots, which corresponds to the condition that the spectrum does not intersect the imaginary axis. This can occur in two cases:

[TABLE]

By combining the conditions in (43c) and (43g), the ACS is of class 0 if and only if

[TABLE]

which can be simplified as

[TABLE]

where

[TABLE]

This proves (38) and statement (i) of the theorem.

Next, the condition for the system (27)–(28) to be of universality class I is that the polynomial (41) has exactly two real roots $\pm\omega_{H}$ and $\gamma(0)>0$. Then, from equation (42), it follows that the criterion for satisfying the inequality $0<\nu_{+}=\omega_{H}^{2}$ is

[TABLE]

or, equivalently, $\left|k_{0}\right|<\left|k^{\tau}_{0}\right|$. This proves statement (ii) of the theorem.

The system belongs to the universality class II if the polynomial (41) has two pairs of real roots $\pm\omega_{H_{1}}$ and $\pm\omega_{H_{2}}$. These exist if the discriminant is positive and if both roots $\nu_{\pm}=\omega_{H_{1,2}}^{2}$, given by equation (42), are positive. This holds under the following conditions:

[TABLE]

Straightforward calculations lead to one of the following conditions holding:

[TABLE]

Combining the conditions in (44c) and (44g) yields

[TABLE]

and either of the following conditions is satisfied:

[TABLE]

Note that in the case $h_{0}^{\tau}=0$, the left-hand side can be considered to be $-\infty$ and $k_{0}$ is not bounded from below. This proves statement (iii) of the theorem.

Finally, for the system to be of universality class U, we require that the strongly unstable spectrum is not empty (unstable), i.e., $k_{0}<0$ or $h_{0}<0$, and the ACS does not cross the imaginary axis (hyperbolic). From the no-crossing conditions in Eqs. (43c) and (43g), we first consider

[TABLE]

which leads to

[TABLE]

Further, if

[TABLE]

we get

[TABLE]

or

[TABLE]

*This confirms statement (iv) of the theorem. *

Figure 3 illustrates the classification provided by Theorem 4.1 in the $(k_{0},h_{0})$-plane for fixed parameter values $k^{\tau}_{0}=1.5$ and $h^{\tau}_{0}=-3$.

The right panel of Fig. 3 also sketches the asymptotic continuous spectrum (blue curves) and the strongly unstable spectrum (red symbols) for the different classes. This classification is independent of delay, with each region corresponding to either a stable state or a specific bifurcation scenario as the time delay varies [51], as we will illustrate below.

Figure 4 illustrates the spectrum of Eq. (31) for fixed time delays $\tau=20$ and different parameter regions corresponding to the classification from Fig. 3.

So one can observe the stable spectrum in the region $S_{0}$, the unstable and hyperbolic spectrum in region $S_{U}$, and the multiple unstable characteristic roots in the regions $S_{I}$ and $S_{II}$ with relatively small real parts.

The delay-independent stability region is given by $S_{0}$. However, the system may also be stable in the other regions (except $S_{U}$) for small or intermediate delays. The stabilization mechanisms in the other regions are determined by specific sequences of bifurcations as described in [51]. To find the exact values of these bifurcations, we substitute $\mu=i\omega$ into the characteristic equation (33),

[TABLE]

By separating the real and imaginary parts of (45), we have

[TABLE]

which can be solved for $k_{0}(\omega)$ and $h_{0}(\omega)$ to obtain the stability boundary parametrically in the $(k_{0},h_{0})$ parameter plane:

[TABLE]

The stability boundaries for fixed time delays $\tau=2$ and $\tau=20$ are shown in Fig. 5.

The stability regions are indicated by the cross-hatched pattern. We observe that the stability region includes the whole region $S_{0}$ and also some parts of the other regions, depending on the time delay. For larger delays, the stability region shrinks to $S_{0}$, while for smaller delays it becomes larger. We also see how the bifurcation lines intersect in the region $S_{II}$, indicating the presence of two pairs of purely imaginary characteristic roots [51]. Note that degenerate points can only appear in the parameter region $S_{II}$.

4.3 Examples with vanishing coupling coefficients

Here we analyze four specific scenarios that illustrate the impact of incomplete communication on the dynamics of active agents. Each example highlights some communication constraints and their influence on the stability of motion and formation of agents.

Example 1: The velocity information is not transmitted from the virtual leader to the active agents, i.e., $h_{0}=h^{\tau}_{0}=0$. The other parameters are fixed as $k_{0}=6$, $k^{\tau}_{0}=0.3$. According to Theorem 4.1, under these conditions, the system (31) falls into universality class II. Figure 6(a) illustrates the stability of this system, represented by the maximal real parts of the eigenvalues, as the delay $\tau$ increases. We observe the intervals in $\tau$, where the system gains stability, while it is unstable for all $\tau\gtrsim 25$. This is an example of stability islands appearing in the parameter space, similar to those in [4, 18, 29].

Example 2: Consider the case $k_{0}=h_{0}=0$, i.e., the active agents communicate with time delays only. Theorem 4.1 implies that then the system (31) is classified as belonging to class I. Figure 6(b) shows a stability region for small delays and an increasing number of unstable eigenvalues with increasing $\tau$.

4.4 Numerical simulation for an example with three agents

We now present numerical simulations of a system with three active agents, which communicate only with the virtual leader, see Fig. 7.

We consider four cases. For cases 1, 2, and 3 we choose the parameters as in example 1 from Fig. 6(a), and for case 4 we choose $k_{0}=2$, $h_{0}=3$, $k_{0}^{\tau}=1.5$, and $h_{0}^{\tau}=1.2$. For all cases, the target trajectory can be set as follows:

[TABLE]

which generates a parabola trajectory for all agents, while a desired isosceles triangular formation configuration. To clarify the motion trajectory, we consider that cases 1, 2, and 3 have $c=1$, and case 4 has $c=10$.

Figure 7 presents the simulation results, which illustrate the formation dynamics of the active agents governed by Eqs. (1)–(2). The control input applied to each agent is defined as

[TABLE]

Additionally, the figure depicts the time evolution of the formation error, as characterized by Eqs. (27)–(28), and the ACS for the error system (21).

For cases 1 to 3, we find that with increasing $\tau$ the system makes an interesting transition from a stable formation configuration to an unstable one and then back again to a stable formation. More explicitly, for $\tau=4.5$ (case 1), all agents asymptotically reach the desired formation configuration, see Fig. 7(a). As $\tau$ increases and crosses a critical instability threshold and at $\tau=5.7$ (case 2), the desired formation configuration is not reached and the agents diverge. The agent trajectories transition to oscillations with exponentially increasing amplitudes and become increasingly unstable, see Fig. 7(b). Then, as $\tau$ continues to increase, the agents regain the ability to achieve the desired formation configuration after crossing a critical stability threshold as shown for $\tau=6.8$ (case 3) in Fig. 7(c). This cycle repeats as $\tau$ continues to change. For case 4, the spectrum is in the absolutely stable region (i.e., class 0), it maintains a stable motion, and the desired formation can always be obtained, for all values of $\tau$, as shown in Fig. 7(d). The converge of the trajectory has a square-wave pattern, often observed in DDEs with the type I continuous spectrum, see e.g. [45].

5 Stability of coupled agents

We now study the case in which agents exchange information with each other, as well as with the virtual leader, see Fig. 8.

In this case, the Laplacian $\mathbf{L}$ in the control function (7) is not vanishing, leading to the the error system (11)–(2). Hence, the MSF approach from Sec. 3.1 can be applied to assess the stability, and the stability problem is reduced via equation (17) to

[TABLE]

where $\lambda_{\ell}$, $\ell=1,\dots,N$ are the eigenvalues of $\mathbf{L}$. In this way, the stability problem is reduced to the stability of each individual mode corresponding to different $\lambda_{\ell}$.

The reduced system (49) allows to describe the structure of the spectrum of the full error system based on the spectra of individual modes.

5.1 Stability of formation with symmetric coupling

We first consider the case of the reciprocal coupling (symmetrical topology) leading to real eigenvalues $\lambda_{\ell}$ of Laplacian $\mathbf{L}$. The main observation here is that equation (49) has the same form as equation (32) for the system coupled only to the virtual leader, but the matrices $M$ and $M^{\tau}$ are replaced by $M-\lambda_{\ell}P$ and $M^{\tau}-\lambda_{\ell}P^{\tau}$, respectively. As a result, Theorem 4.1 can be applied to each coupling mode separately, and each mode can be classified as being in one of the universality classes 0, I, II, or U. Therefore, the spectrum of the whole system (11)–(2) is a union of the spectra of individual modes, and it can be described by the following notation $\text{0}_{i}\text{I}_{j}\text{II}_{m}\text{U}_{s}$, where $i,j,m,$ and $s$ denote the number of modes having the spectrum of class 0, I, II, or U, respectively, see Fig. 9 for illustration.

Let us analyze the special structure of the parameter dependency in system (49). For this, we use our notation $p_{0}=(k_{0},h_{0},k_{0}^{\tau},h_{0}^{\tau})$ and additionally define $\overline{p}:=(k,h,k^{\tau},h^{\tau})$. Since the matrices in (49) are

[TABLE]

we observe that system (49) only depends on $p_{0}+\lambda_{\ell}\overline{p}$. This means that the changes in the parameters $\lambda_{\ell}$ correspond to a line in the parameter space through the point $p_{0}$ and a direction given by $\overline{p}$. More specifically, the following corollary holds.

Corollary 5.1**.**

The following statements holds true for the coupled system (11)–(2).

(i)* The spectrum of eigenvalues of (11)–(2) has the structure $\text{0}_{i}\text{I}_{j}\text{II}_{m}\text{U}_{s}$ with $i,j,m,s\geq 0$ and $i+j+m+s=N$.*

(ii)* The system is absolutely stable if and only if $p_{0}+\lambda_{\ell}\overline{p}\in S_{0}$ for all $\ell=1,\dots,N$, where $S_{0}$ is defined in (38).*

(iii)* The system is unstable and hyperbolic if there exist $\ell$ such that $p_{0}+\lambda_{\ell}\overline{p}\in S_{U}$, where $\ell=1,\dots,N$, and the set $S_{U}$ is defined in Theorem 4.1(iv).*

Figure 10 illustrates how the stability and the spectrum structure in the coupled system can be obtained using Corollary 5.1 and the results for the uncoupled system.

The figure shows two cases: large and intermediate delays. Each dashed line is determined by the point $p_{0}=(k_{0},h_{0},1.5,-3)$ and the direction $\overline{p}=(k,h,0,0)$. The parameter values for each transverse mode $\ell$ are given by $p_{0}+\lambda_{\ell}\overline{p}$ and are indicated by white points in Fig. 10. $\lambda_{\ell}$ are exemplary eigenvalues of the Laplacian matrix $\mathbf{L}$ (10), which correspond to the transverse modes.

The points on the green dashed line of Fig. 10 correspond to the spectrum structure $\text{0}_{3}\text{I}_{1}\text{II}_{0}\text{U}_{0}$, i.e., all transverse modes are stable, but the mode with $\lambda_{0}=0$ is unstable. Hence, the agent formation is stable, while the trajectory can deviate from the virtual leader trajectory for large delay. As the delay becomes smaller ($\tau=2$) in Fig. 10(b), all points including $p_{0}$ fall into the cross-hatched region, and therefore the formation configuration is stable and follows the virtual leader stably. Similarly, the blue dashed line illustrates how the transverse stability changes as the delay increases or decreases. The red dashed line provides an example where the formation configuration cannot be stabilized by any delay or coupling.

To illustrate the effect of different coupling formation configurations, we compute the bifurcation diagram in $\lambda$ and $h_{0}$ parameter space, see Fig. 11.

Since $\lambda$ are the eigenvalues of the Laplacian matrix $\mathbf{L}$, this takes into account different possible coupling formation configurations, i.e., $\lambda$ is the general expression for $\lambda_{\ell}$. Figures Fig. 11(a) and (b) correspond to $k_{0}=0$ and $k_{0}=3$ respectively. The stable active agent formations correspond to the cross-hatched area in the parameter space for fixed delays, and the blue $S_{0}$ region when stability is desired for arbitrary positive delays.

The bifurcation lines create boundaries for the stability region, and are given by the explicit parametric expression

[TABLE]

which can be obtained from the following characteristic equation of the coupled system (49):

[TABLE]

by substituting $\mu=i\omega$ and solving it with respect to the real eigenvalue $\lambda$ (the coupling is symmetric) and $h_{0}$.

Similarly to the bifurcation diagram in Fig. 5, the stability region becomes larger for small delays and shrinks to $S_{0}$ for large delays. As a result, the system can also be stabilized in $S_{I}$ and $S_{II}$ for smaller delays.

5.2 Stability of formation of non-symmetrically coupled agents

We remind that stability of the active agent formation is described by system (49), where $\lambda_{\ell}$ are eigenvalues of the Laplacian matrix $\mathbf{L}$. The corresponding characteristic equation is (52).

For non-symmetrically coupled agents, $\lambda_{\ell}$ are generally complex. Therefore, it is convenient to represent the stability region in the complex plane for $\lambda\in\mathbb{C}$. From the characteristic equation (52), we find

[TABLE]

Then, the explicit expression for $\lambda$ is

[TABLE]

which can be used to find the bifurcation curves and the stability region in the $\lambda$ complex plane parametrically. Specifically, when the eigenvalues $\mu$ are purely imaginary, this corresponds to $\mu=i\omega$, $\omega\in\mathbb{R}$, and the corresponding curves in the parameter space are given parametrically as $\lambda(i\omega)$, $\omega\in\mathbb{R}$.

By plotting the bifurcation curves as outlined above, we obtain the boundary of the stability region, which is illustrated by the green lines in Figs. 12(a) and (b) for two sets of parameters.

Additionally, Figs. 12(a) and (b) show the largest Lyapunov exponent, defined as $\Lambda_{\max}=\max_{j}(\Re(\mu_{j}))$, as a color map. This is computed numerically directly from equation (52). It is clear that the stability boundary $\Lambda_{\max}=0$ corresponds to the analytical green bifurcation line. Furthermore, the dark blue colored areas correspond to regions of absolute stability. Different parameters lead to different shapes of the stable regions, as can be seen in Figs. 12(a) and (b). Furthermore, the red and gray dashed lines in the figures show the contour lines for $\Lambda_{\max}$, which are obtained by setting $\mu=\text{const}+i\omega$ into Eq. (54) and plotting $\lambda(\text{const}+i\omega)$. Detailed parameterized curves for the contour lines on the contour map are shown in Appendix A.

Figures 12(c) and (d) show how the largest Lyapunov exponent $\Lambda_{\max}$ varies with $\Re(\lambda)$, where $\Im(\lambda)=0$ is fixed.

In the remaining part of this section, we will consider the case of $M^{\tau}=P=0$ and a large delay limit. This corresponds to a scenario in which the delay terms occur exclusively in the interaction, and the local feedback is non-delayed ($k=h=k_{0}^{\tau}=h_{0}^{\tau}=0$). The properties of the MSF in this case are reported in [7], and a notable feature is that the stability region asymptotically becomes a circle around the origin in the $\lambda$-plane as the delay increases. We demonstrate this for our system when $M^{\tau}=P=0$, providing the first rigorous proof of this property. Moreover, we calculate the first-order correction terms that cause the bifurcation curve to deviate from a circular shape. To the first approximation, it has a rotated teardrop shape as in Fig. 13.

Theorem 5.2**.**

Let $M^{\tau}=P=0$, $k_{0}>0$, $h_{0}>0$, and $k^{\tau}\neq 0$. Then the bifurcation curve $\lambda(i\omega)$ defined by Eq. (54) has the following properties for sufficiently large $\tau$ and $|\omega\tau|<T$ with some $T>0$ independent on $\tau$.

(i)* The following asymptotic representation holds*

[TABLE]

where $\lambda_{0}=\lambda(0)=-k_{0}/k^{\tau}$.

(ii)* The curve $\lambda(i\omega)$ has self-intersection points at $\lambda(i\omega_{j})=\lambda(-i\omega_{j})\in\mathbb{R}$, where*

[TABLE]

(iii)* The angles $\theta_{j}$ between tangential vectors of the curve at the self-intersection points $\lambda\left({i\omega_{j}}\right)$ are given by*

[TABLE]

(iv)* The stability region of the MSF defined as $\Lambda(\lambda):=\max(\Re(\mu(\lambda)))$ with $\max$ taken over all characteristic roots $\mu$ of equation (52) is confined to the connected region containing the origin.*

The geometric structure described by the theorem is illustrated in Fig. 13.

Proof 5.3**.**

Under the condition $M^{\tau}=P=0$, the characteristic equation (52) is reduced to

[TABLE]

For $\mu=0$, equation (57) reads

[TABLE]

Let $\lambda_{0}$ and $\boldsymbol{\upsilon}_{0}$ be the corresponding eigenvalue and the eigenvector so that

[TABLE]

A direct calculation gives

[TABLE]

Now consider the case $\mu=i\omega~(\omega\neq 0)$. Then condition (57) is equivalent to

[TABLE]

where $\boldsymbol{\upsilon}_{\omega}$ is the right null eigenvector of matrix $D_{\omega}=i\omega I-M+\lambda_{\omega}P^{\tau}e^{-i\omega\tau}$, and $\lambda_{\omega}=\lambda\left(i\omega\right)$. We apply the following transformations

[TABLE]

[TABLE]

where $\beta_{\omega}$ and $\tilde{\boldsymbol{\upsilon}}_{\omega}$ are new unknown function and a vector, which describes the perturbation terms for $\lambda_{\omega}$ and $\boldsymbol{\upsilon}_{\omega}$. By substituting Eqs. (63) and (64) into Eq. (62), we obtain

[TABLE]

Considering the leading order terms $\mathcal{O}(1/\tau)$ of (65), we get

[TABLE]

Equation (66) is solvable if and only if $\left(i\omega\tau I+\beta_{\omega}P^{\tau}\right)\boldsymbol{\upsilon}_{0}\in\text{Im}(D_{0})$, or, equivalently

[TABLE]

The condition (67) can be rewritten using the left null eigenvector $\boldsymbol{\upsilon}^{+}_{0}=\left[h_{0}-h^{\tau}\frac{k_{0}}{k^{\tau}},1\right]$ of the matrix $D_{0}$, which satisfies $\boldsymbol{\upsilon}^{+}_{0}D_{0}=0$. As a result, the solvability condition for (66) leads to the equation

[TABLE]

which can be solved with respect to $\beta_{\omega}$:

[TABLE]

By substituting (69) into (63), we yield (55). This completes the proof of the statement (i) of the theorem.

Next, we determine the offset angle at which the bifurcation curve self-intersects, i.e., at the point where $\lambda(i\omega)=\overline{\lambda(i\omega)}$. Due to $k=h=k_{0}^{\tau}=h_{0}^{\tau}=0$, Eq. (54) simplifies to

[TABLE]

Imposing the self-intersection condition $\lambda(i\omega)=\overline{\lambda(i\omega)}$ yields the requirement that $\lambda(i\omega)$ is real, i.e.,

[TABLE]

from which we obtain

[TABLE]

Defining $\Omega=\omega\tau$ and expanding in $1/\tau$ we arrive at

[TABLE]

which allows us to approximate the solutions using an asymptotic expansion. Let $\Omega=\Omega_{0j}+\frac{1}{\tau}\tilde{\Omega}_{j}$, where $\Omega_{0j}=j\pi$ is the leading-order approximation. Substituting it into Eq. (70) gives

[TABLE]

Taking into account the periodicity of the tangent, we can solve for $\omega_{j}=\frac{\Omega}{\tau}$, yielding

[TABLE]

The proof of the statement (ii) is complete.

We further compute the slope angle $\alpha_{j}$ of the bifurcation curve near the intersection point $\omega_{j}$. This requires evaluating the derivative of $\lambda\left(i\omega\right)$ with respect to $\omega$ at $\omega_{j}$:

[TABLE]

This derivative has the following form:

[TABLE]

where

[TABLE]

Then, the angle $\alpha_{j}$ is given by

[TABLE]

Substituting the expressions for $A_{-1}$ and $B_{0}$, we find

[TABLE]

and hence,

[TABLE]

Because the bifurcation curve is symmetric with respect to the imaginary axis, the offset angle $\theta_{j}$ at the self-intersection points is given by

[TABLE]

This completes the proof of the statement (iii).

*The bifurcation curves $\lambda(i\omega)$ determine the stability boundary where $\Lambda_{\max}=0$. Properties (i))–(iii) imply that these curves encircle the origin. Furthermore the spectrum of all characteristic roots $\mu$ for $\lambda=0$ is equal to the spectrum of $M$ and it is stable under the assumptions of the theorem. Therefore, the stable region is necessarily confined to the connected component containing the origin. This confirms that statement (iv) holds. *

Remark 5.4**.**

*According to Theorem 5.2, the stability region of the MSF is asymptotically circular for large-delay coupling. While this general result was announced in [7], a rigorous proof is given here for the first time, albeit for a particular class of delay-coupled systems. *

Figure 13 illustrates the key properties of the bifurcation curve $\lambda(i\omega)$ as discussed in Theorem 5.2 in the complex $\lambda$-plane for various time delays. Figures 13(a) and 13(b) focus on the case $\tau=3$, (a) zooms in on the stable region and the first self-intersection in (b). Figure 13(a) highlights the curve Eq. (54) (solid line), the asymptotic curve Eq. (55) (dashed line), the starting point $\lambda_{0}$, the first self-intersection point $\lambda_{1}$, tangents (arrowed line), and offset angle $\theta_{1}$, as calculated from Theorem 5.2. Figure 13(b) displays a larger piece of the bifurcation curve when $\tau=3$, emphasizing multiple self-intersection points $\lambda_{j}$ and the corresponding tangent directions and offset angles $\theta_{j}$. Figure 13(c) shows the bifurcation curves for different time delays $\tau$. For $\tau=1$, the curve forms a teardrop shape and has no rotational symmetry. As $\tau$ increases, for instance, at $\tau=5$, the curve starts to become circular. When $\tau\to\infty$, the curve converges to a circle around the origin.

5.3 Numerical simulation examples with three agents

To demonstrate how coupling can stabilize an otherwise unstable agent formation, we performed numerical simulations of a three-agent coupled formation system. Figure 14 presents the motion trajectories of agents, the time evolution of the motion error, and the asymptotic continuous spectrum of the error system.

The motion of agents are governed by Eqs. (1)–(2) with the control input (5). The error system is defined in Eqs. (11)–(2), which is evaluated with respect to the target trajectory. The target trajectory, the parameters $k_{0},h_{0},k^{\tau}_{0},h^{\tau}_{0}$, and time delays are the same as in Fig. 7. The coupling parameters in Eq. (5) are fixed to $k=3,~h=3,~k^{\tau}=-0.5,~\text{and}~h^{\tau}=0$. The coupling structure is determined by the following Laplacian matrix:

[TABLE]

The asymptotic continuous spectrum of the error system is described by (21).

The simulations followed a setup similar to the uncoupled case (see Fig. 7). This allows for a direct comparison of agent formation behavior before and after the introduction of coupling. We observe that a suitable coupling improves the formation stability in two key aspects: (i) the agents converge more reliably to the desired formation configuration; (ii) the amplitude of delay-induced oscillations in their trajectories is significantly reduced.

In Fig. 14, we show the stabilizing effect of inter-agent coupling, which significantly enhances the resilience to a communication delay of agents. Compared to Fig. 7, coupling can reduce delay-induced oscillations and form a more robust formation configuration. After introducing the coupling, the trajectories of agents remain following the virtual leader when delays are $\tau=4.5,6.8$, and $15$, see Figs. 14(a), (c), and (d). And agents can achieve the desired formation configuration for all delays, see Figs. 14(a)–(d). The middle row confirms the above coordinated behavior through the bounded decay of errors. The tracking error $\|\mathbf{e}(t)\|$ (black lines) has a significant reduction compared to Fig. 7, i.e., the oscillations caused by the delay are significantly attenuated. The formation error $\|(\mathbf{L}\otimes 1_{3})\mathbf{e}(t)\|$ (shown by the orange lines) represents the difference in the relative positions of the active agents. When this error asymptotically converges to [math], it indicates that the desired formation shape can be achieved, regardless of the delays. The bottom row shows the spectrum plots of the error system, reflecting the stability corresponding to the eigenvalues $\lambda=4,0,5$ of the coupling matrix $\mathbf{L}$ (74).

6 Conclusions

In summary, this study provides a comprehensive analysis of how time-delayed interactions affect the stability of motion and formation in active agent systems. We consider a general linear model incorporating delay effects into the dynamics, guided by a virtual leader. To quantify deviations in position and velocity of motion, we introduce error variables, which result in a high-dimensional linear DDE model with inertial effects describing the error dynamics. This equation serves as the basis for our stability analysis.

Our investigation focused on two fundamental scenarios: (1) agent motion driven only by a virtual leader (uncoupled case), and (2) agent motion influenced both by a virtual leader and mutual interactions among agents (coupled case). Furthermore, we applied recent results from [59, 51] to provide explicit parametric conditions for the formation to be stable for all delays (we called it absolutely stable), unstable for all delays, or possess explicitly defined stability domains. For this, we used the classification of the spectra of linear DDEs from [51]. All analytical findings were corroborated by numerical simulations, including a representative example of a pattern formation involving three active agents.

In the coupled case, we employ the MSF approach to analyze the formation stability of active agents, which is a powerful tool for studying the stability of coupled error systems by decoupling the dynamics into distinct coupling modes. For the case of complex Laplacian eigenvalues $\lambda$, we assess the stability of agent formation by mapping the largest Lyapunov exponent and identifying stability regions in the $(\Re(\lambda),\Im(\lambda))$-plane, as shown in Fig. 12. An interesting finding is that, as the interacting delay increases, the stability region gradually forms a circular shape centered at the origin. This work provides the first rigorous proof of this circular property, extending the results from [7]. For finite delays, it has a teardrop shape, the properties of which are described in Theorem 5.2.

A key feature of our results is that our theorems provide explicit conditions despite the large number of parameters involved in the control setup. These parameters include the elements of the control matrices $M$, $M^{\tau}$, $P$ and $P^{\tau}$, as well as the potentially arbitrary Laplacian matrix $\mathbf{L}$, which describes the interactions.

Our results could be useful for applications involving multiple agents interacting with time delays, such as in robotics, autonomous systems and networked control. A possible area for future research would be to investigate how adaptive interactions or time-varying delays affect the dynamical system presented here.

Appendix A Parameterized curves for contour lines of the master stability function

In Sec. 5.2 of the main text, we introduced the contour map of the MSF and identified the stable and unstable regions separated by the bifurcation curve. In the following, we provide more detailed parameterized curves for contour lines on the contour map (see Figs. 15(a) and (b)) and determine the main stable and unstable contour lines (see Figs. 15(c) and (d)).