# Liquid-Graph Time-Constant Network for Multi-Agent Systems Control

Antonio Marino, Claudio Pacchierotti, Paolo Robuffo Giordano

► To cite this version:

Antonio Marino, Claudio Pacchierotti, Paolo Robuffo Giordano. Liquid-Graph Time-Constant Network for Multi-Agent Systems Control. 2024. hal-04552893

HAL Id: hal-04552893

https://hal.science/hal-04552893

Preprint submitted on 19 Apr 2024

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 Attribution 4.0 International License
# Liquid-Graph Time-Constant Network for Multi-Agent Systems Control

Antonio Marino, Claudio Pacchierotti, Paolo Robuffo Giordano

► To cite this version:

Antonio Marino, Claudio Pacchierotti, Paolo Robuffo Giordano. Liquid-Graph Time-Constant Network for Multi-Agent Systems Control. CDC 2024 - 63rd IEEE Conference on Decision and Control, IEEE, Dec 2024, Milan (Italie), Italy. pp.793-800, 10.1109/CDC56724.2024.10886394. hal-04552893v2

HAL Id: hal-04552893

https://hal.science/hal-04552893v2

Submitted on 3 Sep 2024

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 Attribution 4.0 International License
# Liquid-Graph Time-Constant Network for Multi-Agent Systems Control

Antonio Marino<fnref n="1" />, Claudio Pacchierotti<fnref n="2" />, Paolo Robuffo Giordano<fnref n="2" />

**Abstract**— In this paper, we propose the Liquid-Graph Time-constant (LGTC) network, a continuous graph neural network (GNN) model for control of multi-agent systems based on the recent Liquid Time Constant (LTC) network. We analyse its stability leveraging contraction analysis and propose a closed-form model that preserves the model contraction rate and does not require solving an ODE at each iteration. Compared to discrete models like Graph Gated Neural Networks (GGNNs), the higher expressivity of the proposed model guarantees remarkable performance while reducing the large amount of communicated variables normally required by GNNs. We evaluate our model on a distributed multi-agent control case study (flocking) taking into account variable communication range and scalability under non-instantaneous communication.

***Index Terms***—Distributed Control, Graph Neural Network, Stability Analysis

## I. INTRODUCTION

In recent years, the interest in multi-agent systems (MAS) has significantly increased [1]. The effective coordination of multiple agents necessitates an understanding of group dynamics. Particularly in multi-robot scenarios, generating individual robot motions relies not only on local sensing but also on information about the group's state, often obtained through communication with a limited subset of nearby team members [2]. As a result, communication becomes a crucial element in achieving distributed solutions for MAS.

Neural networks have become prominent in data-driven control applications due to their outstanding approximation capabilities [3]. In distributed control, neural networks offer particular advantages by efficiently approximating complex distributed policies without the need for intricate optimizations and designs. The literature presents numerous data-driven approaches, with a focus on those rooted in reinforcement learning [4], [5]. However, these approaches operate on local sensing data without inter-agent communication. On the other hand, distributed machine learning employs communication to partition the learning process across multiple machines contributing to group knowledge [6].

Leveraging communication, recent trends in data-driven distributed control involve employing Graph Neural Networks (GNNs) to encode distributed control and planning solutions. In this respect, Gama et al. [7] extend GNNs to flocking control for large teams. Additional examples

of the use of GNNs for distributed control are found in space coverage [8], multi-robot path planning [9], and motion planning [10], including obstacle-rich environments [11]. GNNs also enhance multi-agent perception [12] and enable distributed active information acquisition [13], translating multi-robot information gathering into graph representation and formulating GNN-based decision-making.

In the recent literature, one of the objectives is to make the learning-based method robust and stable [14]. In this context, many works applied contraction analysis to demonstrate recurrent neural network stability [15], or directly closed-loop stability in continuous learning [16] and adaptive control [17]. Recently, the stability analyses presented in [18], [19] showcased the concepts of ISS (Input-to-State Stability) and incremental ISS ($\delta$ISS) [20], [21] in the context of LSTMs and GRUs, which are two of the most popular recurrent neural network models. Inspired by these last results, we proposed [22] the conditions for $\delta$ISS of the recurrent version of GNN, i.e. Gated Graph Neural Networks (GGNN) [23].

In this study, we introduce a continuous graph ordinary differential equation (ODE) network called Liquid-Graph Time-Constant (LGTC) network, inspired by its single-agent counterpart, the Liquid Time Constant (LTC) network [24]. The LTC network demonstrates a strong representation learning capability, empowering agents to extrapolate and make inferences even in unfamiliar scenarios. Additionally, its closed-form version [25] eliminates the need to solve the initial value problem of the ODE, a requirement in other ODE networks [26]–[28]. In the existing literature, there are few attempts to make graph ODEs for modelling dynamic network systems [29], predicting traffic flow [30], or for sequential recommendations [31]. These models typically employ an auto-regressive network over a vector field modelled by a Graph Neural Network (GNN), initializing the internal state with the previous layer's output. In contrast, the proposed LGTC aims to introduce a novel model that learns a dynamic (or “liquid”) time constant, influenced by both input data and the internal state communicated over the agent graph. This dynamic time constant enables each element of the state to capture specific dynamics, thereby enhancing the model's predictive capabilities.

In addition, our contribution includes proposing an approximated closed-form for the LGTC system, enhancing the computational efficiency and reducing the communication load. We proceed to evaluate the proposed models alongside GGNN and GraphODE in a flocking control example.

<footnote n="1">¹ A. Marino is with Univ Rennes, CNRS, Inria, IRISA - Rennes, France. E-mail: antonio.marino@irisafr.</footnote>

<footnote n="2">² C. Pacchierotti and P. Robuffo Giordano are with CNRS, Univ Rennes, Inria, IRISA - Rennes, France. E-mail: {claudio.pacchierotti.prg} @irisafr. This work was supported by the ANR-20-CHIA-0017 project "MULTI-SHARED".</footnote>
# Liquid-Graph Time-Constant Network for Multi-Agent Systems Control

Antonio Marino<fnref n="1" />, Claudio Pacchierotti<fnref n="2" />, Paolo Robuffo Giordano<fnref n="2" />

**Abstract**— In this paper, we propose the Liquid-Graph Time-constant (LGTC) network, a continuous graph neural network (GNN) model for control of multi-agent systems based on the recent Liquid Time Constant (LTC) network. We analyse its stability leveraging contraction analysis and propose a closed-form model that preserves the model contraction rate and does not require solving an ODE at each iteration. Compared to discrete models like Graph Gated Neural Networks (GGNNs), the higher expressivity of the proposed model guarantees remarkable performance while reducing the large amount of communicated variables normally required by GNNs. We evaluate our model on a distributed multi-agent control case study (flocking) taking into account variable communication range and scalability under non-instantaneous communication.

***Index Terms***—Distributed Control, Graph Neural Network, Stability Analysis

## I. INTRODUCTION

Communication is a crucial element in achieving distributed solutions for multi-agent systems (MAS) from control to planning [2]. Like many distributed control algorithms, also learning-based control can benefit from interagent communication to partition the prediction process across multiple machines contributing to task scalability and prediction accuracy [6].

Leveraging communication, recent trends in data-driven distributed control involve employing Graph Neural Networks (GNNs) to encode distributed control and planning solutions. In this respect, Gama et al. [7] extend GNNs to flocking control for large teams. Additional examples of the use of GNNs for distributed control are found in space coverage [8], multi-robot path planning [9], and motion planning [10], including obstacle-rich environments [11]. GNNs also enhance multi-agent perception [12] and enable distributed active information acquisition [13], translating multi-robot information gathering into graph representation and formulating GNN-based decision-making.

In the recent literature, one of the objectives is to make the learning-based method robust and stable [14]. In this context, many works applied contraction analysis to demonstrate recurrent neural network stability [15], or directly closed-loop stability in continuous learning [16] and adaptive control [17]. Recently, the stability analyses presented in [18], [19] showcased the concepts of ISS (Input-to-State Stability)

and incremental ISS ($\delta$ISS) [20], [21] in the context of LSTMs and GRUs, which are two of the most popular recurrent neural network models. Inspired by these last results, we proposed [22] the conditions for $\delta$ISS of the recurrent version of GNN, i.e. Gated Graph Neural Networks (GGNN) [23].

In this study, we introduce a continuous graph ordinary differential equation (ODE) network called Liquid-Graph Time-Constant (LGTC) network, inspired by its single-agent counterpart, the Liquid Time Constant (LTC) network [24]. The LTC network demonstrates a strong representation learning capability, empowering agents to extrapolate and make inferences even in unfamiliar scenarios. Additionally, its closed-form version [25] eliminates the need to solve the initial value problem of the ODE, a requirement in other ODE networks [26]–[28]. In the existing literature, there are few attempts to make graph ODEs for modelling dynamic network systems [29], predicting traffic flow [30], or for sequential recommendations [31]. These models typically employ an auto-regressive network over a vector field modelled by a Graph Neural Network (GNN), initializing the internal state with the previous layer’s output. In contrast, the proposed LGTC aims to introduce a novel model that learns a dynamic (or “liquid”) time constant, influenced by both input data and the internal state communicated over the agent graph. This dynamic time constant enables each element of the state to capture specific dynamics, thereby enhancing the model’s predictive capabilities. The specialized dynamics can be used to reduce the number of the network’s internal state (memory) and consequently the number of information each agent needs to communicate. The final aim of this study is to approach model-based algorithms’ communication efficiency while increasing prediction accuracy, which in control translates into a higher distributed control efficiency.

In addition, our contribution includes proposing an approximated closed-form for the LGTC system, enhancing the computational efficiency and reducing the communication load. We proceed to evaluate the proposed models alongside GGNN and GraphODE in a flocking control example.

## II. PRELIMINARIES

Let $\mathcal{G} = (\mathcal{V}, \mathcal{E})$ be an undirected graph where $\mathcal{V} = \{v_1, \dots, v_N\}$ is the vertex set (representing the $N$ agents in the group) and $\mathcal{E} \subseteq \mathcal{V} \times \mathcal{V}$ is the edge set. Each edge $e_k = (i, j) \in \mathcal{E}$ is associated with a weight $w_{ij} \ge 0$ such that $w_{ij} > 0$ if the agent $i$ and $j$ can interact and $w_{ij} = 0$ otherwise. As usual, we denote with $\mathcal{N}_i = \{j \in \mathcal{V} | w_{ij} > 0\}$ the set of neighbors of agent $i$. We also let $\mathbf{A} \in \mathbb{R}^{N \times N}$

<footnote n="1">¹ A. Marino is with Univ Rennes, CNRS, Inria, IRISA - Rennes, France. E-mail: antonio.marino@irisa.fr.</footnote>

<footnote n="2">² C. Pacchierotti and P. Robuffo Giordano are with CNRS, Univ Rennes, Inria, IRISA - Rennes, France. E-mail: {claudio.pacchierotti.prg}@irisa.fr. This work was supported by the ANR-20-CHIA-0017 project "MULTI-SHARED".</footnote>
## II. PRELIMINARIES

Let $\mathcal{G} = (\mathcal{V}, \mathcal{E})$ be an undirected graph where $\mathcal{V} = \{v_1, \dots, v_N\}$ is the vertex set (representing the $N$ agents in the group) and $\mathcal{E} \subseteq \mathcal{V} \times \mathcal{V}$ is the edge set. Each edge $e_k = (i, j) \in \mathcal{E}$ is associated with a weight $w_{ij} \ge 0$ such that $w_{ij} > 0$ if the agent $i$ and $j$ can interact and $w_{ij} = 0$ otherwise. As usual, we denote with $\mathcal{N}_i = \{j \in \mathcal{V} | w_{ij} > 0\}$ the set of neighbors of agent $i$. We also let $\mathbf{A} \in \mathbb{R}^{N \times N}$ be the adjacency matrix with entries given by the weights $w_{ij}$. Defining the degree matrix $\mathbf{D} = \text{diag}(d_i)$ with $d_i = \sum_{j \in \mathcal{N}_i} w_{ij}$, the Laplacian matrix of the graph is $\mathbf{L} = \mathbf{D} - \mathbf{A}$.

The graph signal $\mathbf{x} \in \mathbb{R}^N$, whose $i$-component $x_i$ is assigned to agent $i$, can be processed over the network by the following linear combination rule applied by each agent

$$ \mathbf{s}_i \mathbf{x} = \sum_{j \in \mathcal{N}_i} s_{ji} (x_i - x_j), \qquad (1) $$

where $\mathbf{s}_i$ is the $i$-th row of $\mathbf{S}$. The signal manipulation can be operated by means of any graph shift operator $\mathbf{S} \subseteq \mathbf{S}$, e.g., Laplacian, adjacency matrix, weighted Laplacian, which respects the sparsity pattern of the graph. Later, in Sect. V, we will use the Laplacian as support matrix, as it is commonly used in distributed control. However, the proposed techniques do not assume the use of a specific support matrix.

Performing $k$ repeated applications of $\mathbf{S}$ on the same signal represents the aggregation of the $k$-hop neighbourhood information. In analogy with traditional signal processing, this property can be used to define a linear graph filtering [32] that processes the multi features signal $\mathbf{x} \in \mathbb{R}^{N \times G}$ with $G$ features:

$$ H_S(\mathbf{x}) = \sum_{k=0}^{K} \mathbf{S}^k \mathbf{x} H_k \qquad (2) $$

where the weights $\mathbf{H}_k \in \mathbb{R}^{G \times F}$ define the output of the filter. Note that $\mathbf{S}^k = \mathbf{S}(\mathbf{S}^{k-1})$, so that it can be computed locally with repeated 1-hop communications between a node and its neighbors. Hence, the computation of $H_S$ is distributed on each node.

### A. Graph Neural Network

Although $H_S$ is simple to evaluate, it can only represent a linear mapping between input and output filters. GNNs increase the expressiveness of the linear graph filters by means of pointwise nonlinearities $\rho: \mathbb{R}^{N \times F_{l-1}} \to \mathbb{R}^{N \times F_{l-1}}$ following a filter bank. Letting $H_{S1}$ be a bank of $F_{l-1} \times F_l$ filters at layer $l$, the GNN layer is defined as

$$ \mathbf{x}_l = \rho(H_{S1}(\mathbf{x}_{l-1})), \quad \mathbf{x}_{l-1} \in \mathbb{R}^{N \times F_{l-1}}. \qquad (3) $$

Starting by $l=0$ with $F_0$, the signal tensor $\mathbf{x}_{l_n} \in \mathbb{R}^{N \times F_{l_n}}$ is the output of a cascade of $l_n$ GNN layers. This specific type of GNN is commonly referred to as a convolution graph network because each layer utilizes a graph signal convolution (2). By the use of imitation learning or similar techniques, the GNN can learn a distributed policy by finding the optimal filter weights $\mathbf{H}_k$ to propagate information among the agents and generate the desired output. Notably, each agent employs an identical version of the network and

exchanges intermediate quantities with the other agents in the network through the application of $\mathbf{S}$, resulting in an overall distributed neural network. GNNs inherit some interesting properties from graph filters, such as permutational equivariance [33] and their local and distributed nature, showing superior ability to process graph signals [9], [13], [34].

### B. Gated Graph Neural Network

Recurrent models of GNNs can solve time-dependent problems. These models, similarly to recurrent neural networks (RNNs), are known as graph recurrent neural networks (GRNNs). GRNNs utilize memory to learn patterns in data sequences, where the data is spatially encoded within graphs, regardless of the team size of the agents [35]. However, traditional GRNNs encounter challenges such as vanishing gradients, which are also found in RNNs. Additionally, they face difficulties in handling long sequences in space, where certain nodes or paths within the graph might be assigned more importance than others in long-range exchanges, causing imbalances in the graph's informational encoding. Forgetting factors can be applied to mitigate this problem, reducing the influence of past or new signal on the state. A Gated Graph Neural Network (GGNN) [23] is a recurrent Graph Neural Network that uses a gating mechanism to control how the past information influences the update of the GNN states. We can add a state and an input gate, $\hat{\mathbf{q}}, \hat{\mathbf{q}} \in \mathcal{Q} \subseteq [0, 1]^{N \times F}$, that are multiplied via the Hadamard product $\circ$ by the state and the inputs of the network, respectively. These two gates regulate how much the past information and the input are used to update the network's internal state. GGNNs admit the following state-space representation [36],

$$ \left\{ \begin{aligned} \tilde{\mathbf{q}} &= \sigma(\tilde{A}_S(\mathbf{x}) + \tilde{B}_S(\mathbf{u}) + \hat{\mathbf{b}}) \\ \hat{\mathbf{q}} &= \sigma(\hat{A}_S(\mathbf{x}) + \hat{B}_S(\mathbf{u}) + \tilde{\mathbf{b}}) \\ \mathbf{x}^+ &= \sigma_c(\hat{\mathbf{q}} \circ A_S(\mathbf{x}) + \tilde{\mathbf{q}} \circ B_S(\mathbf{u}) + \mathbf{b}) \end{aligned} \right. \qquad (4) $$

with $\sigma(x) = \frac{1}{1+e^{-x}}$ being the logistic function, and $\sigma_c(x) = \frac{e^x - e^{-x}}{e^x + e^{-x}}$ being the hyperbolic tangent. $\hat{A}_S, \hat{B}_S$ are graph filters of the forgetting gate, $\tilde{A}_S, \tilde{B}_S$ are graph filters (2) of the input gate, and $A_S$ and $B_S$ are the state graph filters (2). $\hat{\mathbf{b}}, \tilde{\mathbf{b}}, \mathbf{b} \in \mathbb{R}^{N \times F}$ are respectively the biases of the gates and the state built as $\mathbf{1}_N \otimes \mathbf{b}$ with the same bias for every agents. We consider the system under the following assumption

**Assumption 1.** The input $\mathbf{u}$ is unity-bounded: $\mathbf{u} \in \mathcal{U} \subseteq [-1, 1]^{N \times G}$, i.e. $||\mathbf{u}||_\infty \le 1$.

We identify the induced $\infty$-norm as $|| \cdot ||_\infty$. We used the following notation for the filters in (4):

$$ \begin{align} A_{0,K} &\triangleq [A_0, \dots, A_K]^T & S_{I,K} &\triangleq [I, S, \dots, S^K] \nonumber \\ \tilde{A}_{0,K} &\triangleq [\tilde{A}_0, \dots, \tilde{A}_K]^T & B_{0,K} &\triangleq [B_0, \dots, B_K]^T \nonumber \\ \tilde{B}_{0,K} &\triangleq [\tilde{B}_0, \dots, \tilde{B}_K]^T & \hat{A}_{0,K} &\triangleq [\hat{A}_0, \dots, \hat{A}_K]^T \tag{5} \\ \hat{B}_{0,K} &\triangleq [\hat{B}_0, \dots, \hat{B}_K]^T & \hat{B}_{0,K} &\triangleq [\hat{B}_0, \dots, \hat{B}_K]^T \nonumber \end{align} $$

where $K$ is the filters length. Then, in light of assumption (1) and knowing that $||\mathbf{x}||_\infty \le 1$, each gate feature $q_i$ satisfies:

$$ |\hat{q}_i| \le \sigma(||S_{0,K}||_\infty(||\hat{A}_{0,K}||_\infty + ||\hat{B}_{0,K}||_\infty) + ||\hat{b}||_\infty) \triangleq \sigma_{\hat{q}}. \qquad (6) $$
be the adjacency matrix with entries given by the weights $w_{ij}$. Defining the degree matrix $D = \text{diag}(d_i)$ with $d_i = \sum_{j \in N_i} w_{ij}$, the Laplacian matrix of the graph is $L = D - A$.

The graph signal $\boldsymbol{x} \in \mathbb{R}^N$, whose $i$-component $x_i$ is assigned to agent $i$, can be processed over the network by the following linear combination rule applied by each agent

$$ \boldsymbol{s}_i \boldsymbol{x} = \sum_{j \in N_i} s_{ji} (x_i - x_j), \qquad (1) $$

where $s_i$ is the $i$-th row of $S$. The signal manipulation can be operated by means of any graph shift operator $S \subseteq S$, e.g., Laplacian, adjacency matrix, weighted Laplacian, which respects the sparsity pattern of the graph. Later, in Sect. V, we will use the Laplacian as support matrix, as it is commonly used in distributed control. However, the proposed techniques do not assume the use of a specific support matrix.

Performing $k$ repeated applications of $S$ on the same signal represents the aggregation of the $k$-hop neighbourhood information. In analogy with traditional signal processing, this property can be used to define a linear graph filtering [32] that processes the multi features signal $\boldsymbol{x} \in \mathbb{R}^{N \times G}$ with $G$ features:

$$ H_S(\boldsymbol{x}) = \sum_{k=0}^{K} S^k \boldsymbol{x} \boldsymbol{H}_k \qquad (2) $$

where the weights $\boldsymbol{H}_k \in \mathbb{R}^{G \times F}$ define the output of the filter. Note that $S^k = S(S^{k-1})$, so that it can be computed locally with repeated 1-hop communications between a node and its neighbors. Hence, the computation of $H_S$ is distributed on each node.

### A. Graph Neural Network

Although $H_S$ is simple to evaluate, it can only represent a linear mapping between input and output filters. GNNs increase the expressiveness of the linear graph filters by means of pointwise nonlinearities $\rho: \mathbb{R}^{N \times F_{l-1}} \to \mathbb{R}^{N \times F_{l-1}}$ following a filter bank. Letting $H_{S,l}$ be a bank of $F_{l-1} \times F_l$ filters at layer $l$, the GNN layer is defined as

$$ \boldsymbol{x}_l = \rho(H_{S,l}(\boldsymbol{x}_{l-1})), \quad \boldsymbol{x}_{l-1} \in \mathbb{R}^{N \times F_{l-1}}. \qquad (3) $$

Starting by $l=0$ with $F_0$, the signal tensor $\boldsymbol{x}_{l_n} \in \mathbb{R}^{N \times F_{l_n}}$ is the output of a cascade of $l_n$ GNN layers. This specific type of GNN is commonly referred to as a convolution graph network because each layer utilizes a graph signal convolution (2). By the use of imitation learning or similar techniques, the GNN can learn a distributed policy by finding the optimal filter weights $\boldsymbol{H}_k$ to propagate information among the agents and generate the desired output. Notably, each agent employs an identical version of the network and exchanges intermediate quantities with the other agents in the network through the application of $S$, resulting in an overall distributed neural network. GNNs inherit some interesting properties from graph filters, such as permutational equivariance [33] and their local and distributed nature, showing superior ability to process graph signals [9], [13], [34].

### B. Gated Graph Neural Network

Recurrent models of GNNs can solve time-dependent problems. These models, similarly to recurrent neural networks (RNNs), are known as graph recurrent neural networks (GRNNs). GRNNs utilize memory to learn patterns in data sequences, where the data is spatially encoded within graphs, regardless of the team size of the agents [35]. However, traditional GRNNs encounter challenges such as vanishing gradients, which are also found in RNNs. Additionally, they face difficulties in handling long sequences in space, where certain nodes or paths within the graph might be assigned more importance than others in long-range exchanges, causing imbalances in the graph's informational encoding. Forgetting factors can be applied to mitigate this problem, reducing the influence of past or new signal on the state. A Gated Graph Neural Network (GGNN) [23] is a recurrent Graph Neural Network that uses a gating mechanism to control how the past information influences the update of the GNN states. We can add a state and an input gates, $\hat{q}, \hat{q} \in Q \subseteq [0, 1]^{N \times F}$, that are multiplied via the Hadamard product $\circ$ by the state and the inputs of the network, respectively. These two gates regulate how much the past information and the input are used to update the network's internal state. GGNNS admit the following state-space representation [36],

$$ \left\{ \begin{aligned} \hat{q} &= \sigma(\tilde{A}_S(\boldsymbol{x}) + \tilde{B}_S(\boldsymbol{u}) + \hat{b}) \\ \hat{q} &= \sigma(\tilde{A}_S(\boldsymbol{x}) + \tilde{B}_S(\boldsymbol{u}) + \tilde{b}) \\ \boldsymbol{x}^+ &= \sigma_c(\hat{q} \circ A_S(\boldsymbol{x}) + \tilde{q} \circ B_S(\boldsymbol{u}) + b) \end{aligned} \right. \qquad (4) $$

with $\sigma(x) = \frac{1}{1+e^{-x}}$ being the logistic function, and $\sigma_c(x) = \frac{e^x - e^{-x}}{e^x + e^{-x}}$ being the hyperbolic tangent. $\hat{A}_S, \hat{B}_S$ are graph filters of the forgetting gate, $\tilde{A}_S, \tilde{B}_S$ are graph filters (2) of the input gate, and $A_S$ and $B_S$ are the state graph filters (2). $\hat{b}, \tilde{b}, b \in \mathbb{R}^{N \times F}$ are respectively the biases of the gates and the state built as $\mathbf{1}_N \otimes b$ with the same bias for every agents. We consider the system under the following assumption

**Assumption 1.** The input $\boldsymbol{u}$ is unity-bounded: $\boldsymbol{u} \in \mathcal{U} \subseteq [-1, 1]^{N \times G}$, i.e. $||\boldsymbol{u}||_\infty \le 1$.

We identify the induced $\infty$-norm as $|| \cdot ||_\infty$. We used the following notation for the filters in the system (4):

$$ \begin{aligned} A_{0,K} &\triangleq [A_0, \dots, A_K]^T & B_{0,K} &\triangleq [B_0, \dots, B_K]^T \\ \tilde{A}_{0,K} &\triangleq [\tilde{A}_0, \dots, \tilde{A}_K]^T & \hat{A}_{0,K} &\triangleq [\hat{A}_0, \dots, \hat{A}_K]^T \\ \tilde{B}_{0,K} &\triangleq [\tilde{B}_0, \dots, \tilde{B}_K]^T & \hat{B}_{0,K} &\triangleq [\hat{B}_0, \dots, \hat{B}_K]^T \end{aligned} \qquad (5) $$

where $K$ is the filters length. Then, in light of assumption 1 and knowing that $||\boldsymbol{x}||_\infty \le 1$, each gate feature $q_i$ satisfies:

$$ |\hat{q}_i| \le \sigma(||S_{I,K}||_\infty(||\hat{A}_{0,K}||_\infty + ||\hat{B}_{0,K}||_\infty) + ||\hat{b}||_\infty) \triangleq \sigma_{\hat{q}}. \qquad (6) $$

**Assumption 2.** Given any two support matrices $||S_1(t)||_\infty$ and $||S_2(t)||_\infty$, $\forall t \in \mathbb{R}^+$ associated with two different graphs, they are bounded by the same $||\tilde{S}||_\infty$; moreover, they are lower bounded by $||\tilde{S}||_\infty$.

In our previous work [22], we derived the following stability condition
**Assumption 2.** Given any two support matrices $||S_1(t)||_\infty$ and $||S_2(t)||_\infty$, $\forall t \in \mathbb{R}^+$ associated with two different graphs, they are bounded by the same $||\tilde{S}||_\infty$; moreover, they are lower bounded by $||\tilde{S}||_\infty$.

In our previous work [22], we derived the following stability condition

**Theorem 1.** Under Assumptions 1 and 2, a sufficient condition for the system (4) to be $\delta$ISS is $A_\delta \le 1$; where

$$A_\delta \triangleq \sigma_{\hat{\theta}} ||\bar{S}_{0,K}||_\infty ||A_{0,K}||_\infty + \frac{1}{4} ||\bar{S}_{0,K}||_\infty^2 ||\hat{A}_{0,K}||_\infty ||A_{0,K}||_\infty \\ + \frac{1}{4} ||\bar{S}_{0,K}||_\infty^2 ||\tilde{A}_{0,K}||_\infty ||B_{0,K}||_\infty. \quad (7)$$

**Remark 1.** In practice, Assumption 2 is met by limiting the cardinality of $\#N_i < N$ for every agent i in the team. For normalized Laplacian, the assumption is met without any further restrictions on the graph topology.

### III. LIQUID-GRAPH TIME CONSTANT NETWORK

In light of the recent findings on the continuous time neural ODE, we want to propose a novel continuous-time graph neural network: Liquid-Graph Time-Constant (LGTC) network. Inspired by the LTC network proposed in [24], the LGTC is described by the following ODE:

$$\left\{ \begin{aligned} \dot{\mathbf{f}} &= \rho(\hat{A}_S(\mathbf{x}) + \mathbf{b}_\mathbf{x}) + \rho(\hat{B}_S(\mathbf{u}) + \mathbf{b}_\mathbf{u}) \\ \dot{\mathbf{x}} &= -(\mathbf{b} + \mathbf{f}) \circ \mathbf{x} - \sum_{k=1}^K \mathbf{S}^k \mathbf{x} \mathbf{A}_k + \mathbf{f} \circ \sigma_c(B_S(\mathbf{u})) \end{aligned} \right. \quad (8)$$

where $\rho(\mathbf{x}) = \text{ReLU}(\mathbf{x})$. Like in GGNN system, $\hat{A}_S, \hat{B}_S, A_S, B_S$ are graph filters (2) and $\mathbf{b}_\mathbf{x}, \mathbf{b}_\mathbf{u}, \mathbf{b} \in \mathbb{R}^{N \times F}$ are biases of the model. This model disposes of an input-state-dependent varying time-constant allowing single elements of the hidden state to identify specialized dynamical systems for input features arriving at each time step.

In the following results on the system, we used the vector operator $\mathbf{X}$: that rearranges the elements of matrix $\mathbf{X}$ in a vector. For the system in (8), the following lemma holds

**Lemma 1.** If $b_i > 0$ and $\mathbf{A}_k^T \otimes \mathbf{S}^k \ge 0$, for $k = [1, \dots, K]$, the state feature $i$ at time $t$ is bounded in the range $[-1, 1]$ if $\mathbf{x}(0) \in \mathcal{X} \subseteq [-1, 1]^{N \times F}$.

*Proof.* Defining the following Lypunov function

$$V = \frac{1}{2} \mathbf{x}|_{:}^{T} \mathbf{x}|_{:}$$

its time derivative along the system trajectories is

$$\begin{aligned} \dot{\mathbf{V}} = &-\mathbf{x}|_{:}^{T} \operatorname{diag}(\mathbf{b}|_{:}) \mathbf{x}|_{:} - \mathbf{x}|_{:}^{T} \operatorname{diag}(\mathbf{f}|_{:}) \mathbf{x}|_{:} \\ &- \mathbf{x}|_{:}^{T} \sum_{k=1}^{K} \mathbf{A}_{k}^{T} \otimes \mathbf{S}^{k} \mathbf{x}|_{:} + \mathbf{x}|_{:}^{T} \operatorname{diag}(\mathbf{f}|_{:}) \sigma_{c}(B_{S}(\mathbf{u}))|_{:} \end{aligned}$$

By imposing $b_i \ge 0$ and $\sum_{k=1}^K \mathbf{A}_k^T \otimes \mathbf{S}^k \ge 0$ for $k = [1, \dots, K]$, we have

$$\dot{\mathbf{V}} < -\mathbf{x}|_{:}^{T} \operatorname{diag}(\mathbf{f}|_{:}) \mathbf{x}|_{:} + \mathbf{x}|_{:}^{T} \operatorname{diag}(\mathbf{f}|_{:}) \sigma_c(B_S(\mathbf{u}))|_{:}$$

When $x_i > 1$, $\dot{V}_i < 0$, because the term $\mathbf{x}|_{:}^T \operatorname{diag}(\mathbf{f}|_{:}) \sigma_c(B_S(\mathbf{u}))|_{:} \le 1^T \operatorname{diag}(\mathbf{f}|_{:}) \mathbf{1}$ and reminding that $f_i > 0$. The proof is completed by noting that $\dot{V}_i < 0$ when $x_i < -1$ with a similar reasoning. $\square$

We named the vector field in eq. (8) with $F(x, u, S, t): \mathcal{X} \times \mathcal{U} \times \mathcal{S} \times \mathbb{R}_{\ge 0} \to \mathbb{R}^{N \times F}$ and its jacobian with $D_x F = \frac{\partial F}{\partial x}|_\bullet$. Giving the induced infinite log-norm definition $\mu_\infty(\mathbf{X}) = \max_i(x_{ii} + \sum_{j=1, j \ne 1}^n |x_{ij}|)$ [37], the followings results are known [15]:

**Definition III.1.** The vector field $F$ is strongly infinitesimally contractivng if:

$$\mu_\infty(D_x F) < -c \quad (9)$$

With $c > 0$ known as contraction rate.

**Definition III.2.** For a vector map $F$ satisfying III.1 and given $l_u, l_S$ the Lipschitz constants of $F$ in the input $\mathbf{u}$ and $\mathbf{S}$, any two solutions $\mathbf{x}_1$ and $\mathbf{x}_2$ with initial conditions $\mathbf{x}_1(0), \mathbf{x}_2(0)$ and inputs $\mathbf{u}_1(t), \mathbf{u}_2(t)$ satisfy the $\delta$ISS relationship

$$\begin{aligned} \|\mathbf{x}_1(t) - \mathbf{x}_2(t)\|_{\infty} &\le e^{-ct} \|\mathbf{x}_1(0) - \mathbf{x}_2(0)\|_{\infty} \\ &+ \frac{l_u}{c} (1 - e^{-ct}) \sup_{\tau \in [0, t]} \|\mathbf{u}_1(\tau) - \mathbf{u}_2(\tau)\|_{\infty} \\ &+ \frac{l_S}{c} (1 - e^{-ct}) \sup_{\tau \in [0, t]} \|S_1(\tau) - S_2(\tau)\|_{\infty} \end{aligned} \quad (10)$$

For the neural network 8, we provide the following theorem

**Theorem 2.** Under Assumption 2, with $\mathbf{x}(0) \in \mathcal{X}$, system (8) is $\delta$ISS if the following constraints are satisfied

$$c \ge 0; \quad \mathbf{b} \ge 0; \quad \mu_\infty\left(\sum_{k=1}^K \mathbf{A}_k^T \otimes \mathbf{S}^k\right) \ge 0 \quad (11)$$

$$\text{with} \quad c = \|\mathbf{b}\|_{\infty} + \|A_{1,K}\|_{\infty} \|\tilde{S}_{1,K}\|_{\infty} \\ + \|\mathbf{b}_{\mathbf{x}}\|_{\infty} - \|\hat{A}_{0,K}^T\|_{\infty} \|\tilde{S}_{0,K}\|_{\infty}. \quad (12)$$

*Proof.* Let be $\sigma_u = \sigma_c(B_S(\mathbf{u})) \subseteq [-1, 1]^{N \times F}$, the log-norm of the jacobian $\mu_\infty(D_x F)$ has its superior in:

$$\begin{aligned} \sup_x \mu_\infty(D_x F) = &\sup_x \mu_\infty(-\operatorname{diag}(\mathbf{b}|_{:}) - \operatorname{diag}(f|_{:})) \\ &- \sum_{k=1}^K \mathbf{A}_k^T \otimes \mathbf{S}^k + \operatorname{diag}(\sigma_u|_{:} - x|_{:}) \operatorname{diag}(Df) \sum_{k=0}^K \hat{\mathbf{A}}_k^T \otimes \mathbf{S}^k \\ \le &-\|\mathbf{b}\|_{\infty} + \mu_\infty(-\sum_{k=1}^K \mathbf{A}_k^T \otimes \mathbf{S}^k) + \sup_x -\|f|_{:}\|_{\infty} \\ &+ \mu_\infty(\operatorname{diag}(\sigma_u|_{:} - x|_{:}) \operatorname{diag}(Df) \sum_{k=0}^K \hat{\mathbf{A}}_k^T \otimes \mathbf{S}^k) \end{aligned}$$

where $D_x f \in [0, 1]^{NF}$ is the derivative of $f$ in $\mathbf{x}|_{:}$. The last inequality follows by $\mu_\infty(\mathbf{X}) = \|\mathbf{X}\|_\infty$ when $\mathbf{X} > 0$ and by the fact that $\mathbf{b}|_{:} \ge 0$ and $f|_{:} \ge 0$. Moreover, we can note that $-\|f|_{:}\|_\infty \le -\|\rho(\hat{A}_S(\mathbf{x}) + \mathbf{b}_\mathbf{x})\|_\infty$. Giving $\sigma_u|_{:} \in [-1, 1]^{NF}$ and $-\mu_\infty(\mathbf{A}_k^T \otimes \mathbf{S}^k) \le 0$, by lemma 1, $\mathbf{x}|_{:} \in [-1, 1]^{NF}$.
**Theorem 1.** Under Assumptions 1 and 2, a sufficient condition for the system (4) to be $\delta$ISS is $A_\delta \le 1$; where

$$A_{\delta} \triangleq \sigma_{\hat{\delta}} ||\bar{S}_{I,K}||_{\infty} ||A_{0,K}||_{\infty} + \frac{1}{4} ||\bar{S}_{I,K}||_{\infty}^{2} ||\hat{A}_{0,K}||_{\infty} ||A_{0,K}||_{\infty} \\ + \frac{1}{4} ||\bar{S}_{I,K}||_{\infty}^{2} ||\tilde{A}_{0,K}||_{\infty} ||B_{0,K}||_{\infty}. \quad (7)$$

Where $\bar{S}_{I,K}$ is the defined as in eq. (5)

**Remark 1.** In practice, Assumption 2 is met by limiting the cardinality of $\#N_i < N$ for every agent $i$ in the team. For normalized Laplacian, the assumption is met without any further restrictions on the graph topology.

### III. LIQUID-GRAPH TIME CONSTANT NETWORK

In light of the recent findings on the continuous time neural ODE, we want to propose a novel continuous-time graph neural network: Liquid-Graph Time-Constant (LGTC) network. Inspired by the LTC network proposed in [24], the LGTC is described by the following ODE:

$$\left\{ \begin{aligned} \dot{\mathbf{f}} &= \rho(\hat{A}_S(\mathbf{x}) + \mathbf{b}_\mathbf{x}) + \rho(\hat{B}_S(\mathbf{u}) + \mathbf{b}_\mathbf{u}) \\ \dot{\mathbf{x}} &= -(\mathbf{b} + \mathbf{f}) \circ \mathbf{x} - \sum_{k=1}^{K} \mathbf{S}^k \mathbf{x} \mathbf{A}_k + \mathbf{f} \circ \sigma_c(B_S(\mathbf{u})) \end{aligned} \right. \quad (8)$$

where $\rho(\mathbf{x}) = \text{ReLU}(\mathbf{x})$. Like in GGNN system, $\hat{A}_S$, $\hat{B}_S$, $A_S$, $B_S$ are graph filters (2) and $\mathbf{b}_\mathbf{x}, \mathbf{b}_\mathbf{u}, \mathbf{x} \in \mathbb{R}^{N \times F}$ are biases of the model. This model disposes of an input-state-dependent varying time-constant allowing single elements of the hidden state to identify specialized dynamical systems for input features arriving at each time step.

In the following results on the system, we used the vector operator $\mathbf{X}_|$ that rearranges the elements of matrix $\mathbf{X}$ in a vector. For the system in (8), the following lemma holds

**Lemma 1.** If $b_i > 0$ and $\mathbf{A}_k^T \otimes \mathbf{S}^k \ge 0$, for $k = [1, \dots, K]$, the state feature $i$ at time $t$ is bounded in the range $[-1, 1]$ if $\mathbf{x}(0) \in \mathcal{X} \subseteq [-1, 1]^{N \times F}$.

*Proof.* Defining the following Lypunov function

$$V = \frac{1}{2} \mathbf{x}_|^T \mathbf{x}_|$$

Its time derivative along the system trajectories is

$$\dot{V} = -\mathbf{x}_|^T \text{diag}(\mathbf{b}_|) \mathbf{x}_| - \mathbf{x}_|^T \text{diag}(\mathbf{f}_|) \mathbf{x}_| \\ - \mathbf{x}_|^T \sum_{k=1}^K \mathbf{A}_k^T \otimes \mathbf{S}^k \mathbf{x}_| + \mathbf{x}_|^T \text{diag}(\mathbf{f}_|) \sigma_c(B_S(\mathbf{u}))$$

By imposing $b_i \ge 0$ and $\sum_{k=1}^K \mathbf{A}_k^T \otimes \mathbf{S}^k \ge 0$ for $k = [1, \dots, K]$, we have

$$\dot{V} < -\mathbf{x}_|^T \text{diag}(\mathbf{f}_|) \mathbf{x}_| + \mathbf{x}_|^T \text{diag}(\mathbf{f}_|) \sigma_c(B_S(\mathbf{u}))$$

When $x_i > 1$, $\dot{V}_i < 0$, because the term $\mathbf{x}_|^T \text{diag}(\mathbf{f}_|) \sigma_c(B_S(\mathbf{u})) \le 1^T \text{diag}(\mathbf{f}_|)1$ and reminding that $f_i > 0$. The proof is completed by noting that $\dot{V}_i < 0$ when $x_i < -1$ with a similar reasoning. $\square$

We denote the vector field of the dynamic system (8) with $F(x, u, S, t): \mathcal{X} \times \mathcal{U} \times \mathcal{S} \times \mathbb{R}_{\ge 0} \to \mathbb{R}^{N \times F}$ and its jacobian with $D_x F = \frac{\partial F}{\partial x}$. Giving the induced infinite log-norm definition

$\mu_\infty(\mathbf{X}) = \max_i(x_{ii} + \sum_{j=1, j \ne 1}^n |x_{ij}|)$ [37], the followings results are known [15]:

**Definition III.1.** The vector field $F$ is strongly infinitesimally contracting if:

$$\mu_\infty(D_x F) < -c \quad (9)$$

With $c > 0$ known as contraction rate.

**Definition III.2.** For a vector map $F$ satisfying III.1 and given $l_u, l_S$ the Lipschitz constants of $F$ in the input $\mathbf{u}$ and $\mathbf{S}$, any two solutions $\mathbf{x}_1$ and $\mathbf{x}_2$ with initial conditions $\mathbf{x}_1(0), \mathbf{x}_2(0)$ and inputs $\mathbf{u}_1(t), \mathbf{u}_2(t)$ satisfy the $\delta$ISS relationship

$$\begin{aligned} \|\mathbf{x}_1(t) - \mathbf{x}_2(t)\|_{\infty} &\le e^{-ct} \|\mathbf{x}_1(0) - \mathbf{x}_2(0)\|_{\infty} \\ &\quad + \frac{l_u}{c} (1 - e^{-ct}) \sup_{\tau \in [0, t]} \|\mathbf{u}_1(\tau) - \mathbf{u}_2(\tau)\|_{\infty} \\ &\quad + \frac{l_S}{c} (1 - e^{-ct}) \sup_{\tau \in [0, t]} \|S_1(\tau) - S_2(\tau)\|_{\infty} \end{aligned} \quad (10)$$

For the neural network (8), we provide the following theorem

**Theorem 2.** Under Assumption 2, with $\mathbf{x}(0) \in \mathcal{X}$, system (8) is $\delta$ISS if the following constraints are satisfied

$$c \ge 0; \quad \mathbf{b} \ge 0; \quad \mu_\infty\left(\sum_{k=1}^K \mathbf{A}_k^T \otimes \mathbf{S}^k\right) \ge 0 \quad (11)$$

$$\text{with} \quad c = \|\mathbf{b}\|_{\infty} + \|A_{1,K}\|_{\infty} \|\tilde{S}_{1,K}\|_{\infty} \\ + \|\mathbf{b}_{\mathbf{x}}\|_{\infty} - \|\hat{A}_{0,K}^T\|_{\infty} \|\tilde{S}_{0,K}\|_{\infty}. \quad (12)$$

*Proof.* Let $\sigma_u = \sigma_c(B_S(\mathbf{u})) \subseteq [-1, 1]^{N \times F}$, the log-norm of the jacobian $\mu_\infty(D_x F)$ has its superior in:

$$\begin{aligned} \sup_x \mu_\infty(D_x F) &= \sup_x \mu_\infty(-\text{diag}(\mathbf{b}_|) - \text{diag}(\mathbf{f}_|)) \\ &\quad - \sum_{k=1}^K \mathbf{A}_k^T \otimes \mathbf{S}^k + \text{diag}(\sigma_{u|} - \mathbf{x}_|) \text{diag}(D_x \mathbf{f}) \sum_{k=0}^K \hat{\mathbf{A}}_k^T \otimes \mathbf{S}^k \\ &\le -\|\mathbf{b}\|_{\infty} + \mu_\infty(-\sum_{k=1}^K \mathbf{A}_k^T \otimes \mathbf{S}^k) + \sup_x [-\|\mathbf{f}_|\|_{\infty} \\ &\quad + \mu_\infty(\text{diag}(\sigma_{u|} - \mathbf{x}_|) \text{diag}(D_x \mathbf{f}) \sum_{k=0}^K \hat{\mathbf{A}}_k^T \otimes \mathbf{S}^k)] \end{aligned}$$

where $D_x \mathbf{f} \in [0, 1]^{NF}$ is the derivative of $\mathbf{f}$ in $\mathbf{x}_|$. The last inequality follows by $\mu_\infty(\mathbf{X}) = \|\mathbf{X}\|_\infty$ when $\mathbf{X} > 0$ and by the fact that $\mathbf{b}_| \ge 0$ and $\mathbf{f}_| \ge 0$. Moreover, we can note that $-\|\mathbf{f}_|\|_\infty \le -\|\rho(\hat{A}_S(\mathbf{x}) + \mathbf{b}_\mathbf{x})\|_\infty$. Giving $\sigma_{u|} \in [-1, 1]^{NF}$ and $-\mu_\infty(\mathbf{A}_k^T \otimes \mathbf{S}^k) \le 0$, by lemma 1, $\mathbf{x}_| \in [-1, 1]^{NF}$. Therefore, we get

$$\begin{aligned} \mu_\infty(D_x F) &\le - \|\mathbf{b}\|_\infty - \left\| \sum_{k=1}^K \mathbf{A}_k^T \otimes \mathbf{S}^k \right\|_\infty - \\ &\quad \|\rho(\hat{A}_S(\mathbf{x} + \mathbf{b}_\mathbf{x}))\|_\infty + \max_{d \in [-1, 1]^{NF}} 2\mu_\infty(\text{diag}(d) \sum_{k=0}^K \hat{\mathbf{A}}_k^T \otimes \mathbf{S}^k) \\ &\le - \|\mathbf{b}\|_\infty - \left\| \sum_{k=1}^K \mathbf{A}_k^T \otimes \mathbf{S}^k \right\|_\infty - \|\rho(\hat{A}_S(\mathbf{x} + \mathbf{b}_\mathbf{x}))\|_\infty + \\ &\quad 2 \max(\mu_\infty(-\sum_{k=0}^K \hat{\mathbf{A}}_k^T \otimes \mathbf{S}^k), \mu_\infty(\sum_{k=0}^K \hat{\mathbf{A}}_k^T \otimes \mathbf{S}^k)) \end{aligned}$$
Therefore, we get

$$
\begin{align*}
\mu_{\infty}(D_x F) &\le -\|\mathbf{b}\|_{\infty} - \left\|\sum_{k=1}^{K} \mathbf{A}_k^T \otimes \mathbf{S}^k\right\|_{\infty} - \\
&\qquad \left\|\rho(\hat{A}_S(\mathbf{x} + \mathbf{b}_x))\right\|_{\infty} + \max_{d \in [-1, 1]^{NF}} 2\mu_{\infty}(\operatorname{diag}(d) \sum_{k=0}^{K} \hat{\mathbf{A}}_k^T \otimes \mathbf{S}^k) \\
&\le -\|\mathbf{b}\|_{\infty} - \left\|\sum_{k=1}^{K} \mathbf{A}_k^T \otimes \mathbf{S}^k\right\|_{\infty} - \left\|\rho(\hat{A}_S(\mathbf{x} + \mathbf{b}_x))\right\|_{\infty} + \\
&\qquad 2\max(\mu_{\infty}(-\sum_{k=0}^{K} \hat{\mathbf{A}}_k^T \otimes \mathbf{S}^k), \mu_{\infty}(\sum_{k=0}^{K} \hat{\mathbf{A}}_k^T \otimes \mathbf{S}^k))
\end{align*}
$$

To further reduce the previous inequality, we can use the sub
multiplicity property of the infinite norm and the fact that the
last term in the previous inequality is $\le \left\|\sum_{k=1}^{K} \hat{\mathbf{A}}_k^T \otimes \mathbf{S}^k\right\|_{\infty}$
to derive

$$
\mu_{\infty}(D_x F) \leq - \|\mathbf{b}\|_{\infty} - \|A_{1,K}\|_{\infty} \|\tilde{S}_{1,K}\|_{\infty} \\ - \|\mathbf{b}_x\|_{\infty} + \|\hat{A}_{0,K}^T\|_{\infty} \|\tilde{S}_{0,K}\|_{\infty} = -c
$$

From the statement in the Theorem (2), the system is
contractive under the definition III.1 and δISS under the
definition III.2 with

$$
\begin{align}
l_u &= (2\|\hat{B}_{0,K}\|_{\infty} + (\|\hat{A}_{0,K}\|_{\infty}\|\bar{S}_{0,K}\|_{\infty} + \|\hat{B}_{0,K}\|_{\infty}\|\bar{S}_{0,K}\|_{\infty} \nonumber \\
&\quad + \|\mathbf{b}_x\|_{\infty} + \|\mathbf{b}_u\|_{\infty})\|B_{0,K}\|_{\infty}\|\bar{S}_{0,K}\|_{\infty} \nonumber \\
l_S &= \binom{k+1}{2} ((\|\hat{A}_{0,K}\|_{\infty}\|\bar{S}_{0,K}\|_{\infty} + \|\hat{B}_{0,K}\|_{\infty}\|\bar{S}_{0,K}\|_{\infty} \nonumber \\
&\quad + \|\mathbf{b}_x\|_{\infty} + \|\mathbf{b}_u\|_{\infty})\|B_{1,K-1}\|_{\infty} + 2\|\hat{A}_{1,K-1}\|_{\infty} \nonumber \\
&\quad + 2\|\hat{B}_{1,K-1}\|_{\infty})\|\bar{S}_{1,K-1}\|_{\infty} - \|\tilde{S}_{1,K-1}\|_{\infty}\|A_{1,K-1}\|_{\infty}) \tag{13}
\end{align}
$$

We omitted the derivation of $l_u$ and $l_S$ which can be found by
computing $\sup_u \|D_u F\|_{\infty}$, $\sup_S \|D_s F\|_{\infty}$, respectively. $\square$

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

The higher expressivity of the LTC model has been proved
in [24] via principal component analysis on trajectory depth.
The results of this analysis showed that LTC performs at
worst like a discrete RNN model over short trajectories, but
its expressive power does not decrease on longer trajectories.
We can derive the same conclusion for LGTC, which shares
with LTC the main components but over a network of
agents. This higher expressivity allows us to reduce the
communicated state variables. Therefore, we can select a
subset of state features $F' \ll F$ and a subset of input
features $G' \ll G$ to communicate. Denoting by $\mathbf{x}_{F'}$ the
features to communicate and by $\mathbf{x}_{F-F'}$ the features to not
communicate, we can encode the reduced communication in
the equations above by changing the filter bank in eq (2) as
$\sum_{k=0}^{K}[\mathbf{I}_N \mathbf{x}_{F-F'}, \mathbf{S}^k \mathbf{x}_{F'}] \mathbf{H}_k$ with $\mathbf{I}_N$ the identity matrix of
dimension $N$. With these reduced graph filters, Theorem 2
still holds true.

To compute the state of the system in eq. (8) at any time
point T, the neural network must implement and ODE solver,
like Runga-Kutta (RK), that simulates the system starting
from a trajectory x(0), to x(T). When used in a prediction
framework, the neural layer (8) will forward the state x(T)
to the next layers. Since the system is stiff, the authors
in [24] proposed an hybrid solver for LTC networks that has

the additional advantage of reducing the vanishing gradient
problem due to the ODE resolution. We modified this solver
to handle our network as follows:

$$
\begin{equation}
\begin{split}
\mathbf{x}(t + \mathbf{\Delta}) ={}& \frac{\mathbf{x}(t) + \Delta(-\sum_{k=1}^{K} \mathbf{S}^k \mathbf{x}(t) \mathbf{A}_k + f \circ \sigma_c(B_S(\mathbf{u})))}{1 + \Delta(b + f)} \quad (14) \\
& \text{of (14), starting from } t = 0 \text{ until reaching } t = T. \text{ The support}
\end{split}
\end{equation}
$$

With $\Delta = T/n$ being a fixed step size for $n$ evaluations of (14), starting from $t=0$ until reaching $t=T$. The support matrix $\mathbf{S}$ and the input $\mathbf{u}$ are fixed in the simulation interval $[0, T]$. Each iteration of (14) requires to communicate again the network state, e.g. $\mathbf{x}(t+\mathbf{\Delta})$ to evaluate $\mathbf{x}(t+2\mathbf{\Delta})$. Therefore, solving the ODE nullifies the reduction in the communicated variables for the high expressivity of the model. For this reason, we provide in the next section a closed-form solution for LGTC that does not require an ODE solver, thus significantly gaining in computation time and communication load.

IV. CLOSED-FORM APPROXIMATION

By fixing the input and support matrix between the interval
$t = [0, T]$, the ODE solution should satisfy the following
integral:

$$
\mathbf{x}(\mathbf{T})|_{t=T} = \mathbf{x}(0)|_{t=T} + \int_{0}^{T} F(\mathbf{x}(\tau), \mathbf{S}(0), \mathbf{u}(0))d\tau \quad (15)
$$

However, solving this integral is not straightforward due to
the non-linearity in F. We aim to find an approximate form
of (15) that keeps the same dynamic properties of the original
system. Given the following quantities

$$
\begin{align*}
f &= -D_x f \operatorname{diag}(\mathbf{x}(T)) \sum_{k=0}^{K} \hat{\mathbf{A}}_k^T \otimes \mathbf{S}(0)^k + \sum_{k=1}^{K} \mathbf{A}_k^T \otimes \mathbf{S}(0)^k \\
f_x &= \rho(\hat{A}_S(\mathbf{x}(T)) + \mathbf{b}_x) \\
f_\sigma &= \rho(\hat{B}_S(\mathbf{u}(0)) + \mathbf{b}_u) + \rho(\hat{A}_S(\mathbf{x}(0)) + \mathbf{b}_x)
\end{align*}
$$

we propose the following closed-form approximation

$$
\begin{equation}
\begin{split}
\mathbf{x}(\mathbf{T})|_{t=T} ={}& e^{-(\operatorname{diag}(\mathbf{b}_x|_t) + \operatorname{diag}(f_{x|t}^T) + f)^T} \operatorname{diag}(\sigma(2f_\sigma)|_t) \mathbf{x}(0)|_t \\
& + \operatorname{diag}(1 - \sigma(2f_\sigma)|_t) \sigma_u|_t.
\end{split}
\tag{16}
\end{equation}
$$

Reminding that $\sigma_u = \sigma_c(B_S(\mathbf{u}(0)))$, one can easily demon-
strate that $\mathbf{x}(\mathbf{T}) \in \mathcal{X}$ since its expression is a convex
combination with coefficient $\sigma(2f_\sigma)$ and $1 - \sigma(2f_\sigma)$. This
closed-form solution maintains the contraction rate as stated
by the following theorem.

**Theorem 3.** With $t = [0, T]$, the state trajectory (16) is a closed-form approximation of the dynamic system (8) with the same contraction rate $c$ under the constraints of Theorem 2.

*Proof.* We start the proof by highlighting the relationship between the superior of the induced log-norm and the superior of the time derivative of the induced norm of the matrix exponential

$$
\sup_x \mu_{\infty}(D_x F) = \frac{\partial}{\partial t} \sup_x \| \exp(D_x Ft) \|_{\infty} \Big|_{t=0^+} \quad (17)
$$

Therefore, we can show that the equivalence between the
right-hand side of the previous relation and the induced
To further reduce the previous inequality, we can use the sub multiplicity property of the infinite norm and the fact that the last term in the previous inequality is $\le \sum_{k=1}^{K} \hat{A}_k^T S^k \|_{\infty}$ to derive

$$ \mu_{\infty}(D_x F) \le -\|b\|_{\infty} - \|A_{1,K}\|_{\infty} \|\tilde{S}_{1,K}\|_{\infty} \\ - \|b_x\|_{\infty} + \|\hat{A}_{0,K}^T\|_{\infty} \|\tilde{S}_{0,K}\|_{\infty} = -c $$

From the statement in the Theorem 2, the system is contractive under the definition III.1 and δISS under the definition III.2 with

$$ l_u = (2\|\hat{B}_{0,K}\|_{\infty} + (\|\hat{A}_{0,K}\|_{\infty} \|\tilde{S}_{0,K}\|_{\infty} + \|\hat{B}_{0,K}\|_{\infty} \|\tilde{S}_{0,K}\|_{\infty} \\ + \|b_x\|_{\infty} + \|b_u\|_{\infty}) \|B_{0,K}\|_{\infty}) \|\tilde{S}_{0,K}\|_{\infty} \\ l_S = \begin{pmatrix} k+1 \\ 2 \end{pmatrix} (((\|\hat{A}_{0,K}\|_{\infty} \|\tilde{S}_{0,K}\|_{\infty} + \|\hat{B}_{0,K}\|_{\infty} \|\tilde{S}_{0,K}\|_{\infty} \\ + \|b_x\|_{\infty} + \|b_u\|_{\infty}) \|B_{1,K-1}\|_{\infty} + 2\|\hat{A}_{1,K-1}\|_{\infty} \\ + 2\|\hat{B}_{1,K-1}\|_{\infty}) \|\tilde{S}_{1,K-1}\|_{\infty} - \|\tilde{S}_{1,K-1}\|_{\infty} \|A_{1,K-1}\|_{\infty}) \quad (13) $$

We omitted the derivation of $l_u$ and $l_S$ which can be found by computing $\sup_u \|D_u F\|_{\infty}$, $\sup_S \|D_S F\|_{\infty}$, respectively. □

The higher expressivity of the LTC model has been proved in [24] via principal component analysis on trajectory depth. The results of this analysis showed that LTC performs at worst like a discrete RNN model over short trajectories, but its expressive power does not decrease on longer trajectories. We can derive the same conclusion for LGTC, which shares with LTC the main components but over a network of agents. This higher expressivity allows us to reduce the communicated state variables. Therefore, we can select a subset of state features $F' \ll F$ and a subset of input features $G' \ll G$ to communicate. Denoting by $x_{F'}$ the features to communicate and by $x_{F-F'}$ the features to not communicate, we can encode the reduced communication in the equations above by changing the filter bank in eq (2) as $\sum_{k=0}^{K} [I_N x_{F-F'}, S^k x_{F'}] H_k$ with $I_N$ the identity matrix of dimension $N$. With these reduced graph filters, Theorem 2 still holds true.

To compute the state of the system in eq. (8) at any time point $T$, the neural network must implement and ODE solver, like Runga-Kutta (RK), that simulates the system starting from a trajectory $x(0)$, to $x(T)$. When used in a prediction framework, the neural layer (8) will forward the state $x(T)$ to the next layers. Since the system is stiff, the authors in [24] proposed an hybrid solver for LTC networks that has the additional advantage of reducing the vanishing gradient problem due to the ODE resolution. We modified this solver to handle our network as follows:

$$ x(t + \Delta) = \frac{x(t) + \Delta(-\sum_{k=1}^{K} S^k x(t) A_k + f \circ \sigma_c(B_S(u)))}{1 + \Delta(b + f)} \quad (14) $$

With $\Delta = T/n$ being a fixed step size for $n$ evaluations of (14), starting from $t=0$ until reaching $t=T$. The support matrix $S$ and the input $u$ are fixed in the simulation interval $[0, T]$. Each iteration of (14) requires to communicate again the network state, e.g. $x(t + \Delta)$ to evaluate $x(t + 2\Delta)$. Therefore, solving the ODE nullifies the reduction in the

communicated variables for the high expressivity of the model. For this reason, we provide in the next section a closed-form solution for LGTC that does not require an ODE solver, thus significantly gaining in computation time and communication load.

## IV. CLOSED-FORM APPROXIMATION

By fixing the input and support matrix between the interval $t = [0, T]$, the ODE solution should satisfy the following integral:

$$ x(T) = x(0) + \int_{0}^{T} F(x(\tau), S(0), u(0)) d\tau \quad (15) $$

However, solving this integral is not straightforward due to the non-linearity in $F$. We aim to find an approximate form of eq. (15) that keeps the same dynamic properties of the original system. Given the following quantities

$$ \begin{aligned} f_i &= -D_x f \operatorname{diag}(x(T)_i) \sum_{k=0}^{K} \hat{A}_k^T \otimes S(0)^k + \sum_{k=1}^{K} A_k^T \otimes S(0)^k \\ f_x &= \rho(\hat{A}_S(x(T)) + b_x) \\ f_\sigma &= \rho(\hat{B}_S(u(0)) + b_u) + \rho(\hat{A}_S(x(0)) + b_x) \end{aligned} \quad (16) $$

we propose the following closed-form approximation

$$ x(T)_i = e^{-(\operatorname{diag}(b_i) + \operatorname{diag}(f_x)_i + f_i)^T \operatorname{diag}(\sigma(2f_\sigma))} x(0)_i + \operatorname{diag}(1 - \sigma(2f_\sigma)) \sigma_{u_i} \quad (16) $$

Reminding that $\sigma_u = \sigma_c(B_S(u(0)))$, one can easily demonstrate that $x(T) \in \mathcal{X}$ since its expression is a convex combination with coefficient $\sigma(2f_\sigma)$ and $1 - \sigma(2f_\sigma)$. This closed-form solution maintains the contraction rate as stated by the following theorem.

**Theorem 3.** With $t = [0, T]$, the state trajectory (16) is a closed-form approximation of the dynamic system (8) with the same contraction rate $c$ of the Theorem 2.

*Proof.* The theorem can be demonstrated by showing that the closed-loop system maintains the same dynamic properties. We start the proof by highlighting the relationship between the superior of the induced log-norm and the superior of the time derivative of the induced norm of the matrix exponential

$$ \sup_x \mu_\infty(D_x F) = \frac{\partial}{\partial t} \sup_x \|\exp(D_x F t)\|_{\infty} \Big|_{t=0^+} \quad (17) $$

Therefore, we can show that the equivalence between the right-hand side of the previous relation and the induced norm of the jacobian ($D_x x(T)$) of eq. (16). Let's name $F_{ex}$ the argument of the exponential in the eq. (16). $\|D_x x(t)\|_{\infty}$
norm of the jacobian ($D_x x(T)$) of (16). Let's name $F_{ex}$ the argument of the exponential in the (16). $||D_x x(t)||_\infty$ satisfies the following:

$$ \begin{aligned} ||D_x x(t)|:||_\infty = & ||e^{F_{ex}} \operatorname{diag}(\sigma(2f_\sigma))+:| + \\ & \operatorname{diag}(x:)e^{F_{ex}}(\operatorname{diag}(\sigma(2f_\sigma)):) \sum_{k=0}^{K} \hat{A}_k^T \otimes S^k t + \\ & \operatorname{diag}(x:)e^{F_{ex}}2\operatorname{diag}(D_\sigma:) \sum_{k=0}^{K} \hat{A}_k^T \otimes S^k - \\ & \operatorname{diag}(\sigma_u)2\operatorname{diag}(D_\sigma:) \sum_{k=0}^{K} \hat{A}_k^T \otimes S^k ||_\infty + \\ & \le ||e^{F_{ex}}||_\infty + ||e^{F_{ex}} \sum_{k=0}^{K} \hat{A}_k^T \otimes S^k t||_\infty + \\ & 0.5(||e^{F_{ex}}||_\infty + 1)|| \sum_{k=0}^{K} \hat{A}_k^T \otimes S^k ||_\infty \\ & \le ||e^{F_{ex}}||_\infty + ||e^{F_{ex}} \sum_{k=0}^{K} \hat{A}_k^T \otimes S^k t||_\infty + || \sum_{k=0}^{K} \hat{A}_k^T \otimes S^k ||_\infty \end{aligned} \quad (18) $$

where the last inequalities come from reasoning similar to the proof in the theorem 2 and the Lipschitz constant of $\sigma(2x)$ equal to 0.5. The following time derivative proves the equality (17)

$$ \frac{\partial}{\partial t} \sup_x ||D_x x(t)||_\infty \Big|_{t=0^+} = -||b||_\infty - ||A_{1,K}||_\infty ||\tilde{S}_{1,K}||_\infty - \\ ||b_x||_\infty + ||\hat{A}_{0,K}||_\infty ||\tilde{S}_{0,K}||_\infty \quad (19) $$

Moreover, one can easily demonstrate that the Lipschitz constant in $u$ and $S$ is upper bounded by $l_u$, $l_S$ respectively, taking into account that $te^{-(F)t}; F > 0$ can be upper bounded by 1. □

With a fixed $T$, the state trajectory (16) can serve as a discrete-time system that in training and evaluation requires one iteration/communication to be evaluated, against the $n$ iteration required by the continuous system. However, equation (16) can not be directly distributed over the agents due to the kronacker product terms appearing in $f$. Moreover, notably, the exponential term in recursion causes vanishing gradient problems in training. Therefore, we modified the trajectory (16) to tackle these two issues. First, we multiply the kronacker product terms by $x/(x + \epsilon)$, with $/$ being the division element by element and a small $\epsilon$. Second, we replace the exponential term with a sigmoidal nonlinearity which decreases much more smoothly. The resulting closed-form graph time constant (CfGC) system is:

$$ \left\{ \begin{array}{l} f_{\sigma} = \rho(\hat{B}_{S}(u) + b_{u}) + \rho(\hat{A}_{S}(x) + b_{x}) \\ f_{x} = \rho(\hat{A}_{S}(x) + b_{x}) \\ f = -D_{x}f\hat{A}_{S}(x) + \sum_{k=1}^{K} S^{k}x\mathbf{A}_{k}/(x + \epsilon) \\ \mathbf{x}^{+} = (x \circ \sigma(-\mathbf{b} + f_{x} + f)t + \pi) - \sigma_{u} \circ \sigma(2f_{\sigma}) + \sigma_{u} \end{array} \right. \quad (20) $$

This system still satisfies the theorem 3 for $\epsilon \to 0$ since $\sigma(-Ft+\pi)|_{t=0} \approx 1$. The model in (20) is different from the one in (8), even if, when they have the same weight matrices,

<figure><img src="image_3.png" /><figcaption>Fig. 1: **Flocking control:** a group of agents (yellow dots) move in order to reach the same velocity and to avoid collision. The leader (red dot) moves in order to reach the target (blue cross) and avoids the collision with the other agents.</figcaption></figure>

they have approximately the same dynamic properties. We demonstrate this good approximation in the result Section V.

## V. VALIDATION EXAMPLE

We experimentally validated the neural networks proposed in a flocking control example. Despite its simplicity, this scenario has been previously studied [38] using different approaches such as GNN, LSTM, allowing for a direct comparison with the proposed model. We also explicitly compare GraphODE [29], GGNN, LGTC and CfGC, satisfying the conditions in Theorem 1, 2 and 3 respectively. The stability condition is imposed by the following regularization in addition to the training loss:

$$ \Pi = \sum_{i} \mathrm{Softplus}(p_i) \quad (21) $$

With $p$ having only one element equal to $\delta A_i - 1$ to satisfy the Theorem 1 and equal to $[-c, -b, -\mu_\infty(\sum_{k=1}^K A_k^T \otimes S^k)]$ for Theorem 2 and 3. The Softplus is introduced to have a smooth ReLU function that can be regulated through its $\beta$ to enforce a higher contraction rate and consequently fast convergence. We used $\beta = 10$.

In the following, we show a case study involving flocking control (Fig. 1) with a leader. In this problem, the agents are initialized to follow random velocities while the goal is to have them all travelling at the same velocity while avoiding collisions with each other. Moreover, one of the agents takes the leader role, conducting the team toward a target unknown to the other agents. Flocking is a canonical problem in decentralized robotics [39], [40].

We considered $N$ agents described by the position $r(t) \in \mathbb{R}^{N \times 2}$ and the velocity $v(t) \in \mathbb{R}^{N \times 2}$ with a double integrator dynamics

$$ r(t+1) = r(t) + T v(t); \quad v(t+1) = v(t) + T u(t); $$

with the discrete acceleration $u(t) \in \mathbb{R}^{N \times 2}$ taken as system input. Note that the agent dynamics is used for building the dataset and for simulation purposes, but it is not provided to
satisfies the following :

$$
\begin{align}
\|D_x \mathbf{x}(t)\|_{\infty} &= \|e^{F_{ex}} \operatorname{diag}(\sigma(2f_{\sigma})\| + \nonumber \\
& \quad \operatorname{diag}(\mathbf{x}_l)e^{F_{ex}}(\operatorname{diag}(\sigma(2f_{\sigma})\| \sum_{k=0}^{K} \hat{\mathbf{A}}_k^T \otimes \mathbf{S}^k t + \nonumber \\
& \quad \operatorname{diag}(\mathbf{x}_l)e^{F_{ex}} 2 \operatorname{diag}(D_{\sigma_l}) \sum_{k=0}^{K} \hat{\mathbf{A}}_k^T \otimes \mathbf{S}^k - \nonumber \\
& \quad \operatorname{diag}(\sigma_u) 2 \operatorname{diag}(D_{\sigma_u}) \sum_{k=0}^{K} \hat{\mathbf{A}}_k^T \otimes \mathbf{S}^k\|_{\infty} \nonumber \\
&\le \|e^{F_{ex}}\|_{\infty} + \|e^{F_{ex}} \sum_{k=0}^{K} \hat{\mathbf{A}}_k^T \otimes \mathbf{S}^k t\|_{\infty} + \nonumber \\
& \quad 0.5(\|e^{F_{ex}}\|_{\infty} + 1)\| \sum_{k=0}^{K} \hat{\mathbf{A}}_k^T \otimes \mathbf{S}^k\|_{\infty} \nonumber \\
&\le \|e^{F_{ex}}\|_{\infty} + \|e^{F_{ex}} \sum_{k=0}^{K} \hat{\mathbf{A}}_k^T \otimes \mathbf{S}^k t\|_{\infty} + \| \sum_{k=0}^{K} \hat{\mathbf{A}}_k^T \otimes \mathbf{S}^k\|_{\infty} \tag{18}
\end{align}
$$

where the last inequalities come from similar reasons to
the proof in the theorem 2 and the Lipschitz constant of
$\sigma(2x)$ equal to 0.5. The following time derivative proves the
equality (17)

$$
\frac{\partial}{\partial t} \sup_x \|D_x \mathbf{x}(\mathbf{t})\|_{\infty} \Big|_{t=0^+} = - \|\mathbf{b}\|_{\infty} - \|A_{1,K}\|_{\infty} \|\tilde{S}_{1,K}\|_{\infty} - \\
\|\mathbf{b}_x\|_{\infty} + \|\hat{A}_{0,K}\|_{\infty} \|\tilde{S}_{0,K}\|_{\infty} \tag{19}
$$

Moreover, one can easily demonstrate that the Lipschitz
constant in **u** and **S** is upper bounded by $l_u$, $l_S$ respectively,
taking into account that $te^{-(F)t}$; $F > 0$ can be upper
bounded by 1.
□

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

With a fixed T, the state trajectory (16) can serve as a
discrete-time system that in training and evaluation requires
one iteration/communication to be evaluated, against the n
iteration required by the continuous system. However, equa-
tion (16) can not be directly distributed over the agents due
to the kronacker product terms appearing in fᵢ. Moreover,
notably, the exponential term in recursion causes vanishing
gradient problems in training. Therefore, we modified the
trajectory (16) to tackle these two issues. First, we multiply
the kronacker product terms by x/(x + ε), with / being
the division element by element and a small ε. Second, we
replace the exponential term with a sigmoidal nonlinearity
which decreases much more smoothly. The resulting closed-
form graph time constant (CfGC) system is:

$$
\begin{equation}
\begin{cases}
f_{\sigma} = \rho(\hat{B}_{S}(\mathbf{u}) + \mathbf{b}_{u}) + \rho(\hat{A}_{S}(\mathbf{x}) + \mathbf{b}_{x}) \\
f_{x} = \rho(\hat{A}_{S}(\mathbf{x}) + \mathbf{b}_{x}) \\
f_{i} = -D_{x}f\hat{A}_{S}(\mathbf{x}) + \sum_{k=1}^{K} \mathbf{S}^{k}\mathbf{x}\mathbf{A}_{k}/(\mathbf{x} + \epsilon) \\
\mathbf{x}^{+} = (\mathbf{x} \circ \sigma(-\mathbf{b} + f_{x} + f_{i})t + \pi) - \sigma_{u} \circ \sigma(2f_{\sigma}) + \sigma_{u}
\end{cases}
\tag{20}
\end{equation}
$$

This system still satisfies the theorem 3 for $\epsilon \to 0$ since
$\sigma(-Ft+\pi)|_{t=0} \approx 1$. The model in (20) is different from the
one in (8), even if, when they have the same weight matrices,
they have approximately the same dynamic properties. We
demonstrate this good approximation in the result Section V.

<figure><img src="image_3.png" /><figcaption>Fig. 1: **Flocking control:** a group of agents (yellow dots) move in order to reach the same velocity and to avoid collision. The leader (red dot) moves in order to reach the target (blue cross) and avoids the collision with the other agents.</figcaption></figure>

V. VALIDATION EXAMPLE

We experimentally validated the neural networks proposed
in a flocking control example. Despite its simplicity, this
scenario has been previously studied [38] using different
approaches such as GNN, LSTM, allowing for a direct
comparison with the proposed model. We also explicitly
compare GraphODE [29], GGNN, LGTC and CfGC, satis-
fying the conditions in Theorem 1, 2 and 3 respectively. The
stability condition is imposed by the following regularization
in addition to the training loss:

$$
\Pi = \sum_i \text{Softplus}(p_i) \quad (21)
$$

With $p_i$ being only one element equal to $\delta A_i - 1$ to satisfy
the Theorem 1 and equal to $[-c, -b, -\mu_\infty(\sum_{k=1}^K A_k^T \otimes S^k)]$
for Theorem 2 and 3. The Softplus is introduced to have
a smooth ReLU function that can be regulated through its
$\beta$ to enforce a higher contraction rate and consequently fast
convergence. We used $\beta = 10$.

In the following, we show a case study involving flocking
control (Fig. 1) with a leader. In this problem, the agents
are initialized to follow random velocities while the goal
is to have them all travelling at the same velocity while
avoiding collisions with each other. Moreover, one of the
agents takes the leader role, conducting the team toward a
target unknown to the other agents. Flocking is a canonical
problem in decentralized robotics [39], [40].

We considered *N* agents described by the position **r**(*t*) ∈ ℝ*N×2* and the velocity **v**(*t*) ∈ ℝ*N×2* with a double integrator dynamics

$$
\mathbf{r}(t+1) = \mathbf{r}(t) + T\mathbf{v}(t); \quad \mathbf{v}(t+1) = \mathbf{v}(t) + T\mathbf{u}(t);
$$

with the discrete acceleration **u**(*t*) ∈ ℝ*N×2* taken as system
input. Note that the agent dynamics is used for building the
dataset and for simulation purposes, but it is not provided to
the learning algorithm. The flocking expert controller [38]
for the follower agents is given by

$$
\mathbf{u}_f(t) = \frac{1}{N} \sum_{i=1}^{N} \mathbf{v}_i(t) - \nabla_r C A(\mathbf{r}(t), \mathbf{r}_j(t)) \Big|_{j=1,...,N} \quad (22a)
$$
the learning algorithm. The flocking expert controller [38] for the follower agents is given by

$$u_f(t) = \frac{1}{N} \sum_{i=1}^{N} v_i(t) - \nabla_r CA(r(t), r_j(t))|_{j=1...N} \quad (22a)$$

and for the leader it is

$$u_l(t) = -W_p(r_l(t) - d(t)) - \nabla_{r_l} CA(r_l(t), r_j(t))|_{j=1...N} \quad (22b)$$

where $W_p$ is a gain, $r_l \in \mathbb{R}^2$ is the leader position, $\nabla_r CA(r(t), r_j(t))/\nabla_{r_l} CA(r_l(t), r_j(t))$ are the gradient of the collision avoidance potential with respect to the position of the agents/leader $r/r_l$, evaluated at the position $r(t)/r_l(t)$ and the position of every other agent $r_j(t)$ at time $t$. The i-element of $\nabla_r CA$ for each robot $i$ with respect to robot $j$ is given by [41]

$$\nabla_{r_i} CA(r_{ij}) = \begin{cases} -\frac{r_{ij}}{\|r_{ij}\|_2^4} - \frac{r_{ij}}{\|r_{ij}\|_2^2} & \text{if } \|r_{ij}\|_2^2 \le R_{CA} \\ \mathbf{0} & \text{otherwise} \end{cases} \quad (23)$$

with $r_{ij} = r_i - r_j$ and $R_{CA} > 0$ indicating the minimum acceptable distance between agents. This potential function is a non-negative, non-smooth function that goes to infinity when the distance reduces and grows when the distance exceeds $R_{CA}$, in order to avoid the team losing connectivity [41]. $u_f(t), u_l(t)$ are a centralized controller since computing them requires agent $i$ to have instantaneous evaluation of $\frac{1}{N} \sum_{i=1}^N v_i(t)$ and $r_j(t)$ of every other agent $j$ in the team. $R_{CA}$ and $W_p$ are tunable parameters of the controllers.

## Neural Network Architecture

We assume that the agents form a communication graph when they are in a sphere of radius $R$ between each others and that exchanges occur at the sampling time $T = 0.05$ s, so that the action clock and the communication clock coincide. The input features vector $w_i \in \mathbb{R}^{10}$ of the robot $i$ for the designed neural network is

$$w_i = \begin{bmatrix} v_i, \sum_{j \in \mathcal{N}S_i} \frac{r_{ij}}{\|r_{ij}\|_2^4}, \sum_{j \in \mathcal{N}S_i} \frac{r_{ij}}{\|r_{ij}\|_2^2}, \\ \{0_2, r_l - d\}, \{[0, 1], [1, 0]\} \end{bmatrix} \quad (24)$$

where $\mathcal{N}S_i$ is the set of the sensing agents within a sphere of radius $R_{CA}$ centred in the robot $i$. Moreover, the vector $w_i$ contains the zero vector $0_2 \in \mathbb{R}^{1 \times 2}$ and the one-hot encoding $[0, 1]$, if the agent is a follower, while $(r_l - d)$ and $[1, 0]$, if the agent is a leader. We chose the one-hot encoding instead of the binary one because it allows differentiating the neural network weights between the leader and the follower. Note that we assume all the information in the vector $w_i$ to be locally available at the sampling/control time $T$.

The core of the neural network for the flocking control is a layer of GGNN, GraphODE, LGTC or CfGC with $F = 50$ features in the hidden state and filter length $K = 2$. Note that the choice of $K$ affects the complexity of the stability

condition imposed since it will constrain more parameters. For LGTC and GraphODE, we implement the hybrid solver (14) and RK4, respectively. The input features are first processed by a cascade of two fully connected layers of 128 nodes before feeding the graph neural network. A readout of two layers with 128 nodes combines the F-features GNN hidden state to get the bidimensional control $u$ saturated to the maximum admissible control. The input layers and the readout that encapsulate the graph neural network shape a more realistic setting to test the stability of the graph neural network that is usually used in combination with other kinds of neural models. To reduce the communicated variables, we used $F' = 4$ and $G = 4$ for a subset of the state and the input to communicate. Following the Remark 1, we used normalized Laplacian as a support matrix.

## Training

We collected a dataset by recording 60 trajectories, further separated into three subsets of training, validation and test set using the proportion 70% - 10% - 20%, respectively. Each trajectory is generated by randomly positioning the agents in a square such that their inter-distance is between 0.6 m and 1.0 m and that their initial velocities are picked at random from the interval $[-2, 2]$ m/s in each direction. The leader is randomly selected among the agents and the target position is randomly located within a square of length 20 m centered at the location of the leader. Regardless of the target location, the trajectories have a duration of 2.5 s and input saturation at 5 m/s$^2$. Moreover, the 60 trajectories are recorded with a random number of agents among $N = [4, 6, 10, 12, 15]$. We fixed the communication range to $R = 4$ m and the sensing to $R_{CA} = 1$ m. We trained the models for 120 epochs and executed the DAGGER algorithm [42] every 20 epochs. The algorithm evaluates the expert controller in (22) on the enrolled state trajectories applying the learned control and adding them to the training set. Note that, thanks to the use of DAGGER, we do not need a large dataset. We solve the imitation learning problem using the ADAM algorithm [43] with a learning rate $1e-3$ and forgetting factors 0.9 and 0.999. The loss function used for imitation learning is the mean squared error between the output of the model and the optimal control action.

## Results

In Fig. 2, we show a comparison between GraphODE, GGNN, LGTC and CfGC controllers for the flocking control case. We evaluate the 4 controllers on 2 sets of experiments with 20 trajectories each. In the experiments, we varied team size and communication range to test the robustness of the controllers. Figure 2 reports the leader position error evaluated after a fixed time of 2.5s with respect to the leader starting location, i.e. $e_f/e_s$ with $e_f, e_s$ respectively being the final and the initial square distance of the leader from the target. We also show the average flocking error in the interval $[0, 2.5s]$ in logarithmic scale. We assume that the communication happens within the sampling time $T$.
and for the leader it is

$$ \mathbf{u}_l(t) = -W_p(\mathbf{r}_l(t) - \mathbf{d}(t)) - \nabla_{\mathbf{r}_l} CA(\mathbf{r}_l(t), \mathbf{r}_j(t))|_{j=1...N} \quad (22b) $$

where $W_p$ is a gain, $\mathbf{r}_l \in \mathbb{R}^2$ is the leader position, $\nabla_{\mathbf{r}_l} CA(\mathbf{r}_l(t), \mathbf{r}_j(t))/\nabla_{\mathbf{r}_l} CA(\mathbf{r}_l(t), \mathbf{r}_j(t))$ are the gradient of the collision avoidance potential with respect to the position of the agents/leader $\mathbf{r}/\mathbf{r}_l$, evaluated at the position $\mathbf{r}(t)/\mathbf{r}_l(t)$ and the position of every other agent $\mathbf{r}_j(t)$ at time $t$. The i-element of $\nabla_{\mathbf{r}_l} CA$ for each robot $i$ with respect to robot $j$ is given by [41]

$$ \nabla_{\mathbf{r}_i} CA(\mathbf{r}_{ij}) = \begin{cases} -\frac{\mathbf{r}_{ij}}{|\mathbf{r}_{ij}|^2} - \frac{\mathbf{r}_{ij}}{|\mathbf{r}_{ij}|^2} & \text{if } |\mathbf{r}_{ij}|^2 \le R_{CA} \\ \mathbf{0} & \text{otherwise} \end{cases} \quad (23) $$

with $\mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j$ and $R_{CA} > 0$ indicating the minimum acceptable distance between agents. This potential function is a non-negative, non-smooth function that goes to infinity when the distance reduces and grows when the distance exceeds $R_{CA}$, in order to avoid the team losing connectivity [41]. $\mathbf{u}_f(t), \mathbf{u}_l(t)$ are a centralized controller since computing them requires agent $i$ to have instantaneous evaluation of $\frac{1}{N} \sum_{i=1}^N \mathbf{v}_i(t)$ and $\mathbf{r}_j(t)$ of every other agent $j$ in the team. $R_{CA}$ and $W_p$ are tunable parameters of the controllers.

## Neural Network Architecture

We assume that the agents form a communication graph when they are in a sphere of radius $R$ between each others and that exchanges occur at the sampling time $T = 0.05$ s, so that the action clock and the communication clock coincide. The input features vector $\mathbf{w}_i \in \mathbb{R}^{10}$ of the robot $i$ for the designed neural network is

$$ \mathbf{w}_i = \left[ \mathbf{v}_i, \sum_{j \in \mathcal{N}S_i} \frac{\mathbf{r}_{ij}}{|\mathbf{r}_{ij}|^2}, \sum_{j \in \mathcal{N}S_i} \frac{\mathbf{r}_{ij}}{|\mathbf{r}_{ij}|^2}, \{\mathbf{0}_2, \mathbf{r}_l - \mathbf{d}\}, \{[0, 1], [1, 0]\} \right] \quad (24) $$

where $\mathcal{N}S_i$ is the set of the sensing agents within a sphere of radius $R_{CA}$ centred in the robot $i$. Moreover, the vector $\mathbf{w}_i$ contains the zero vector $\mathbf{0}_2 \in \mathbb{R}^{1 \times 2}$ and the one-hot encoding $[0, 1]$, if the agent is a follower, while $(\mathbf{r}_l - \mathbf{d})$ and $[1, 0]$, if the agent is a leader. We chose the one-hot encoding instead of the binary one because it allows differentiating the neural network weights between the leader and the follower. Note that we assume all the information in the vector $\mathbf{w}_i$ to be locally available at the sampling/control time $T$.

The core of the neural network for the flocking control is a layer of GGNN, GraphODE, LGTC or CfGC with $F = 50$ features in the hidden state and filter length $K = 2$. Note that the choice of $K$ affects the complexity of the stability condition imposed since it will constrain more parameters. For LGTC and GraphODE, we implement the hybrid solver (14) and RK4, respectively. The input features are first processed by a cascade of two fully connected layers of 128 nodes before feeding the graph neural network. A readout of two layers with 128 nodes combines the F-features GNN

hidden state to get the bidimensional control $\mathbf{u}$ saturated to the maximum admissible control. The input layers and the readout that encapsulate the graph neural network shape a more realistic setting to test the stability of the graph neural network that is usually used in combination with other kinds of neural models. To reduce the communicated variables, we used $F' = 4$ and $G = 4$ for a subset of the state and the input to communicate. Following the Remark 1, we used normalized Laplacian as a support matrix.

## Training

We collected a dataset by recording 60 trajectories, further separated into three subsets of training, validation and test set using the proportion 70% - 10% - 20%, respectively. Each trajectory is generated by randomly positioning the agents in a square such that their inter-distance is between 0.6 m and 1.0 m and that their initial velocities are picked at random from the interval $[-2, 2]$ m/s in each direction. The leader is randomly selected among the agents and the target position is randomly located within a square of length 20 m centered at the location of the leader. Regardless of the target location, the trajectories have a duration of 2.5 s and input saturation at 5 m/s$^2$. Moreover, the 60 trajectories are recorded with a random number of agents among $N = [4, 6, 10, 12, 15]$. We fixed the communication range to $R = 4$ m and the sensing to $R_{CA} = 1$ m. We trained the models for 120 epochs and executed the DAGGER algorithm [42] every 20 epochs. The algorithm evaluates the expert controller in (22) on the enrolled state trajectories applying the learned control and adding them to the training set. Note that, thanks to the use of DAGGER, we do not need a large dataset. We solve the imitation learning problem using the ADAM algorithm [43] with a learning rate $1e-3$ and forgetting factors 0.9 and 0.999. The loss function used for imitation learning is the mean squared error between the output of the model and the optimal control action.

## Results

In Fig. 2, we show a comparison between GraphODE, GGNN, LGTC and CfGC controllers for the flocking control case. We evaluate the 4 controllers on 2 sets of experiments with 20 trajectories each. In the experiments, we varied team size and communication range to test the robustness of the controllers. Figure 2 reports the leader position error evaluated after a fixed time of 2.5s with respect to the leader starting location, i.e. $e_f/e_s$ with $e_f, e_s$ respectively being the final and the initial square distance of the leader from the target. We also show the average flocking error in the interval $[0, 2.5s]$ in logarithmic scale. We assume that the communication happens within the sampling time $T$.

In the first experiment, the controller scalability is evaluated for team sizes $N = [4, 10, 25, 50]$ with a communication range at 4 m. Figure 2 shows empirically the higher prediction abilities of LGTC and CfGC compared to GGNN, as they are up to 40% for the flocking error and 10% for the leader error closer to the expert controller. This is despite LGTC and CfGC communicating fewer state and input
<figure><img src="image_4.png" /><figcaption>Fig. 2: Flocking and Leader Error for GGNN, LGCT, CfGC and GraphODE, varying the team size N with a communication range of 4m and a variable communication range with N = 25.</figcaption></figure>

In the first experiment, the controller scalability is evaluated for team sizes $N = [4, 10, 25, 50]$ with a communication range at 4 m. Figure 2 shows empirically the higher prediction abilities of LGTC and CfGC compared to GGNN, as they are up to 40% for the flocking error and 10% for the leader error closer to the expert controller. This is despite LGTC and CfGC communicating fewer state and input variables as explained in the previous section. Moreover, the performances of LGTC and CfGC are close to each other, confirming the good approximation of CfGC. Sometimes, for example, for the flocking error with $N = 10$, we can see CfGC performing better than LGTC of $1-2\%$. This can be justified by the fact that CfGC is simpler to train, as it is a discrete model, and so it reduces the vanishing gradient phenomena. In general, GraphODE is the one that performs the worst both for leader and flocking errors. This can be explained easily by no direct use of the input features, as they are only used to initialize the internal state, resulting in an auto-regressive model.

In the second set of experiments, we evaluate the performances of the networks when the communication range differs from the training range $R = 4$ m. For 2 m all the models show a drop in the performances with 0.18 m/s in the flocking error for GGNN, LGTC and CfGC, and 0.3 m/s for GraphODE. This drop appears due to the lower ability to broadcast the information among the agents when they space out over 2 m. As expected, when the communication range increases, the flocking errors decrease for all controllers, since they can communicate with more agents at the same time. However, in this case, the leader error increases, since the leader is often “encapsulated” by the other team agents and it is thus forced to follow them to not collide, causing a slower convergence to the target. Note that this behaviour also affects the expert controller. Notably, also for this experiment, LGTC and CfGC have better performances compared to the others and are similar to each other.

## VI. CONCLUSIONS

In this work, we model a new continuous-time graph neural network, Liquid-Graph Time-constant (LGTC). Its enhanced prediction capabilities allow us to reduce communication load, one of the main drawbacks of GNNs. We analyze its dynamic properties leveraging contraction analysis and find a closed-form approximation (CfGC) that

preserves the dynamic properties of the system. We validate the proposed model in a flocking control case and compare it with $\delta$ISS GGNN and the GraphODE. Results suggest that, despite the reduced communicated vector, LGTC and CfGC have performances closer to the centralized expert and more robust to parametric changes in a deployment scenario, such as communication radius and team size. Future works will consider how to further reduce the communication load, still high compared to model-based approaches, as well as an experimental validation on quadrotor UAVs.

## REFERENCES

[1] A. Dorri, S. S. Kanhere, and R. Jurdak, "Multi-agent systems: A survey," *Ieee Access*, vol. 6, pp. 28 573-28 593, 2018.

[2] J. Cortés and M. Egerstedt, "Coordinated control of multi-robot systems: A survey," *SICE Journal of Control, Measurement, and System Integration*, vol. 10, no. 6, pp. 495-503, 2017.

[3] G.-B. Huang, L. Chen, C. K. Siew, et al., "Universal approximation using incremental constructive feedforward networks with random hidden nodes," *IEEE Trans. Neural Networks*, vol. 17, no. 4, pp. 879-892, 2006.

[4] S. Batra, Z. Huang, A. Petrenko, T. Kumar, A. Molchanov, and G. S. Sukhatme, "Decentralized control of quadrotor swarms with end-to-end deep reinforcement learning," in *Conference on Robot Learning. PMLR*, 2022, pp. 576-586.

[5] C. He, Y. Wan, Y. Gu, and F. L. Lewis, "Integral reinforcement learning-based multi-robot minimum time-energy path planning subject to collision avoidance and unknown environmental disturbances," *IEEE Control Systems Letters*, vol. 5, no. 3, pp. 983-988, 2020.

[6] N. Majcherczyk, N. Srishankar, and C. Pinciroli, "Flow-fl: Data-driven federated learning for spatio-temporal predictions in multi-robot systems," in *2021 IEEE International Conference on Robotics and Automation (ICRA)*. IEEE, 2021, pp. 8836-8842.

[7] F. Gama, E. Tolstaya, and A. Ribeiro, "Graph neural networks for decentralized controllers," in *ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*. IEEE, 2021, pp. 5260-5264.

[8] Q. Li, W. Lin, Z. Liu, and A. Prorok, "Message-aware graph attention networks for large-scale multi-robot path planning," *IEEE Robotics and Automation Letters*, vol. 6, no. 3, pp. 5533-5540, 2021.

[9] E. Tolstaya, J. Paulos, V. Kumar, and A. Ribeiro, "Multi-robot coverage and exploration using spatial graph neural networks," in *2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)*. IEEE, 2021, pp. 8944-8950.

[10] A. Khan, A. Ribeiro, V. Kumar, and A. G. Francis, "Graph neural networks for motion planning," *arXiv preprint arXiv:2006.06248*, 2020.

[11] X. Ji, H. Li, Z. Pan, X. Gao, and C. Tu, "Decentralized, unlabeled multi-agent navigation in obstacle-rich environments using graph neural networks," in *2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)*. IEEE, 2021, pp. 8936-8943.

[12] Y. Zhou, J. Xiao, Y. Zhou, and G. Loianno, "Multi-robot collaborative perception with graph neural networks," *IEEE Robotics and Automation Letters*, vol. 7, no. 2, pp. 2289-2296, 2022.
<figure><img src="image_4.png" /><figcaption>Fig. 2: Flocking and Leader Error for GGNN, LGCT, CfGC and GraphODE, varying the team size N with a communication range of 4m and a variable communication range with N = 25.</figcaption></figure>

variables as explained in the previous section. Moreover, the performances of LGTC and CfGC are close to each other, confirming the good approximation of CfGC. Sometimes, for example, for the flocking error with $N = 10$, we can see CfGC performing better than LGTC of 1–2%. This can be justified by the fact that CfGC is simpler to train, as it is a discrete model, and so it reduces the vanishing gradient phenomena. In general, GraphODE is the one that performs the worst both for leader and flocking errors. This can be explained easily by no direct use of the input features, as they are only used to initialize the internal state, resulting in an auto-regressive model.

In the second set of experiments, we evaluate the performances of the networks when the communication range differs from the training range $R = 4$ m. For 2 m all the models show a drop in the performances with 0.18 m/s in the flocking error for GGNN, LGTC and CfGC, and 0.3 m/s for GraphODE. This drop appears due to the lower ability to broadcast the information among the agents when they space out over 2 m. As expected, when the communication range increases, the flocking errors decrease for all controllers, since they can communicate with more agents at the same time. However, in this case, the leader error increases, since the leader is often “encapsulated” by the other team agents and it is thus forced to follow them to not collide, causing a slower convergence to the target. Note that this behaviour also affects the expert controller. Notably, also for this experiment, LGTC and CfGC have better performances compared to the others and are similar to each other.

## VI. CONCLUSIONS

In this work, we model a new continuous-time graph neural network, Liquid-Graph Time-constant (LGTC). Its enhanced prediction capabilities allow us to reduce communication load, one of the main drawbacks of GNNs. We analyze its dynamic properties leveraging contraction analysis and find a closed-form approximation (CfGC) that preserves the dynamic properties of the system. We validate the proposed model in a flocking control case and compare it with $\delta$ISS GGNN and the GraphODE. Results suggest that, despite the reduced communicated vector, LGTC and CfGC have performances closer to the centralized expert and more robust to parametric changes in a deployment scenario, such as communication radius and team size. Future works will

consider how to further reduce the communication load, still high compared to model-based approaches, as well as an experimental validation on quadrotor UAVs.

## REFERENCES

[1] A. Dorri, S. S. Kanhere, and R. Jurdak, “Multi-agent systems: A survey,” *Ieee Access*, vol. 6, pp. 28 573–28 593, 2018.

[2] J. Cortés and M. Egerstedt, “Coordinated control of multi-robot systems: A survey,” *SICE Journal of Control, Measurement, and System Integration*, vol. 10, no. 6, pp. 495–503, 2017.

[3] G.-B. Huang, L. Chen, C. K. Siew, et al., “Universal approximation using incremental constructive feedforward networks with random hidden nodes,” *IEEE Trans. Neural Networks*, vol. 17, no. 4, pp. 879–892, 2006.

[4] S. Batra, Z. Huang, A. Petrenko, T. Kumar, A. Molchanov, and G. S. Sukhatme, “Decentralized control of quadrotor swarms with end-to-end deep reinforcement learning,” in *Conference on Robot Learning. PMLR*, 2022, pp. 576–586.

[5] C. He, Y. Wan, Y. Gu, and F. L. Lewis, “Integral reinforcement learning-based multi-robot minimum time-energy path planning subject to collision avoidance and unknown environmental disturbances,” *IEEE Control Systems Letters*, vol. 5, no. 3, pp. 983–988, 2020.

[6] N. Majcherczyk, N. Srishankar, and C. Pinciroli, “Flow-fl: Data-driven federated learning for spatio-temporal predictions in multi-robot systems,” in *2021 IEEE International Conference on Robotics and Automation (ICRA)*. IEEE, 2021, pp. 8836–8842.

[7] F. Gama, E. Tolstaya, and A. Ribeiro, “Graph neural networks for decentralized controllers,” in *ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)*. IEEE, 2021, pp. 5260–5264.

[8] Q. Li, W. Lin, Z. Liu, and A. Prorok, “Message-aware graph attention networks for large-scale multi-robot path planning,” *IEEE Robotics and Automation Letters*, vol. 6, no. 3, pp. 5533–5540, 2021.

[9] E. Tolstaya, J. Paulos, V. Kumar, and A. Ribeiro, “Multi-robot coverage and exploration using spatial graph neural networks,” in *2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)*. IEEE, 2021, pp. 8944–8950.

[10] A. Khan, A. Ribeiro, V. Kumar, and A. G. Francis, “Graph neural networks for motion planning,” *arXiv preprint arXiv:2006.06248*, 2020.

[11] X. Ji, H. Li, Z. Pan, X. Gao, and C. Tu, “Decentralized, unlabeled multi-agent navigation in obstacle-rich environments using graph neural networks,” in *2021 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)*. IEEE, 2021, pp. 8936–8943.

[12] Y. Zhou, J. Xiao, Y. Zhou, and G. Loianno, “Multi-robot collaborative perception with graph neural networks,” *IEEE Robotics and Automation Letters*, vol. 7, no. 2, pp. 2289–2296, 2022.

[13] M. Tzes, N. Bousias, E. Chatzipantazis, and G. J. Pappas, “Graph neural networks for multi-robot active information acquisition,” *arXiv preprint arXiv:2209.12091*, 2022.

[14] L. Hewing, K. P. Wabersich, M. Menner, and M. N. Zeilinger, “Learning-based model predictive control: Toward safe learning in control,” *Annual Review of Control, Robotics, and Autonomous Systems*, vol. 3, pp. 269–296, 2020.

[15] A. Davydov, S. Jafarpour, and F. Bullo, “Non-euclidean contraction theory for robust nonlinear stability,” *IEEE Transactions on Automatic Control*, vol. 67, no. 12, pp. 6667–6681, 2022.
[13] M. Tzes, N. Bousias, E. Chatzipantazis, and G. J. Pappas, "Graph neural networks for multi-robot active information acquisition," arXiv preprint arXiv:2209.12091, 2022.

[14] L. Hewing, K. P. Wabersich, M. Menner, and M. N. Zeilinger, "Learning-based model predictive control: Toward safe learning in control," *Annual Review of Control, Robotics, and Autonomous Systems*, vol. 3, pp. 269–296, 2020.

[15] A. Davydov, S. Jafarpour, and F. Bullo, "Non-euclidean contraction theory for robust nonlinear stability," *IEEE Transactions on Automatic Control*, vol. 67, no. 12, pp. 6667–6681, 2022.

[16] B. Song, J.-J. Slotine, and Q.-C. Pham, "Stability guarantees for continuous r1 control," *arXiv preprint arXiv:2209.07324*, 2022.

[17] H. Tsukamoto, S.-J. Chung, and J.-J. Slotine, "Learning-based adaptive control using contraction theory," in *2021 60th IEEE Conference on Decision and Control (CDC)*. IEEE, 2021, pp. 2533–2538.

[18] F. Bonassi, M. Farina, and R. Scattolini, "On the stability properties of gated recurrent units neural networks," *Systems & Control Letters*, vol. 157, p. 105049, 2021.

[19] E. Terzi, F. Bonassi, M. Farina, and R. Scattolini, "Learning model predictive control with long short-term memory networks," *International Journal of Robust and Nonlinear Control*, vol. 31, no. 18, pp. 8877–8896, 2021.

[20] F. Bayer, M. Bürger, and F. Allgöwer, "Discrete-time incremental iss: A framework for robust nmpc," in *2013 European Control Conference (ECC)*. IEEE, 2013, pp. 2068–2073.

[21] W. D'Amico, A. La Bella, and M. Farina, "An incremental input-to-state stability condition for a generic class of recurrent neural networks," *arXiv preprint arXiv:2210.09721*, 2022.

[22] A. Marino, C. Pacchierotti, and P. R. Giordano, "Input state stability of gated graph neural networks," *IEEE Transactions on Control of Network Systems*, pp. 1–12, 2024.

[23] L. Ruiz, F. Gama, and A. Ribeiro, "Gated graph recurrent neural networks," *IEEE Transactions on Signal Processing*, vol. 68, pp. 6303–6318, 2020.

[24] R. Hasani, M. Lechner, A. Amini, D. Rus, and R. Grosu, "Liquid time-constant networks," in *Proceedings of the AAAI Conference on Artificial Intelligence*, vol. 35, no. 9, 2021, pp. 7657–7666.

[25] R. Hasani, M. Lechner, A. Amini, L. Liebenwein, A. Ray, M. Tschaikowski, G. Teschl, and D. Rus, "Closed-form continuous-time neural networks," *Nature Machine Intelligence*, vol. 4, no. 11, pp. 992–1003, 2022.

[26] E. Haber and L. Ruthotto, "Stable architectures for deep neural networks," *Inverse problems*, vol. 34, no. 1, p. 014004, 2017.

[27] R. T. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud, "Neural ordinary differential equations," *Advances in neural information processing systems*, vol. 31, 2018.

[28] M. Lechner, R. Hasani, M. Zimmer, T. A. Henzinger, and R. Grosu, "Designing worm-inspired neural networks for interpretable robotic control," in *2019 international conference on robotics and automation (ICRA)*. IEEE, 2019, pp. 87–94.

[29] M. Poli, S. Massaroli, J. Park, A. Yamashita, H. Asama, and J. Park, "Graph neural ordinary differential equations," *arXiv preprint arXiv:1911.07532*, 2019.

[30] Y. Qin, W. Ju, H. Wu, X. Luo, and M. Zhang, "Learning graph ode for continuous-time sequential recommendation," *IEEE Transactions on Knowledge and Data Engineering*, 2024.

[31] Y. Su, B. Ren, and K. Zhang, "Graph ode recurrent neural networks for traffic flow forecasting," in *2022 IEEE 5th International Conference on Electronics and Communication Engineering (ICECE)*. IEEE, 2022, pp. 178–182.

[32] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vanderheynst, "The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains," *IEEE Signal Processing Magazine*, vol. 30, no. 3, pp. 83–98, 2013.

[33] F. Gama, J. Bruna, and A. Ribeiro, "Stability properties of graph neural networks," *IEEE Transactions on Signal Processing*, vol. 68, pp. 5680–5695, 2020.

[34] Q. Li, F. Gama, A. Ribeiro, and A. Prorok, "Graph neural networks for decentralized path planning," in *Proc. International Conference on Autonomous Agents and Multiagent Systems*, 2020, pp. 1901–1903.

[35] F. Gu, H. Chang, W. Zhu, S. Sojoudi, and L. El Ghaoui, "Implicit graph neural networks," *Advances in Neural Information Processing Systems*, vol. 33, pp. 11984–11995, 2020.

[36] D. Lukovnikov, J. Lehmann, and A. Fischer, "Improving the long-range performance of gated graph neural networks," *arXiv preprint arXiv:2007.09668*, 2020.

[37] T. Ström, "On logarithmic norms," *SIAM Journal on Numerical Analysis*, vol. 12, no. 5, pp. 741–753, 1975.

[38] E. Tolstaya, F. Gama, J. Paulos, G. Pappas, V. Kumar, and A. Ribeiro, "Learning decentralized controllers for robot swarms with graph neural networks," in *Conference on robot learning*. PMLR, 2020, pp. 671–682.

[39] W. Yu, G. Chen, and M. Cao, "Distributed leader-follower flocking control for multi-agent dynamical systems with time-varying velocities," *Systems & Control Letters*, vol. 59, no. 9, pp. 543–552, 2010.

[40] L. E. Beaver and A. A. Malikopoulos, "An overview on optimal flocking," *Annual Reviews in Control*, vol. 51, pp. 88–99, 2021.

[41] H. G. Tanner, A. Jadbabaie, and G. J. Pappas, "Stable flocking of mobile agents part i: dynamic topology," in *42nd IEEE International Conference on Decision and Control (IEEE Cat. No. 03CH37475)*, vol. 2. IEEE, 2003, pp. 2016–2021.

[42] S. Ross, G. Gordon, and D. Bagnell, "A reduction of imitation learning and structured prediction to no-regret online learning," in *Proceedings of the fourteenth international conference on artificial intelligence and statistics*. JMLR Workshop and Conference Proceedings, 2011, pp. 627–635.

[43] D. P. Kingma and J. Ba, "Adam: A method for stochastic optimization," *CoRR*, vol. abs/1412.6980, 2014.
[16] B. Song, J.-J. Slotine, and Q.-C. Pham, "Stability guarantees for continuous rl control," arXiv preprint arXiv:2209.07324, 2022.

[17] H. Tsukamoto, S.-J. Chung, and J.-J. Slotine, "Learning-based adaptive control using contraction theory," in 2021 60th IEEE Conference on Decision and Control (CDC). IEEE, 2021, pp. 2533-2538.

[18] F. Bonassi, M. Farina, and R. Scattolini, "On the stability properties of gated recurrent units neural networks," Systems & Control Letters, vol. 157, p. 105049, 2021.

[19] E. Terzi, F. Bonassi, M. Farina, and R. Scattolini, "Learning model predictive control with long short-term memory networks," *International Journal of Robust and Nonlinear Control*, vol. 31, no. 18, pp. 8877-8896, 2021.

[20] F. Bayer, M. Bürger, and F. Allgöwer, "Discrete-time incremental is: A framework for robust nmpc," in *2013 European Control Conference (ECC)*. IEEE, 2013, pp. 2068-2073.

[21] W. D'Amico, A. La Bella, and M. Farina, "An incremental input-to-state stability condition for a generic class of recurrent neural networks," *arXiv preprint arXiv:2210.09721*, 2022.

[22] A. Marino, C. Pacchierotti, and P. R. Giordano, "Input state stability of gated graph neural networks," *IEEE Transactions on Control of Network Systems*, pp. 1-12, 2024.

[23] L. Ruiz, F. Gama, and A. Ribeiro, "Gated graph recurrent neural networks," *IEEE Transactions on Signal Processing*, vol. 68, pp. 6303-6318, 2020.

[24] R. Hasani, M. Lechner, A. Amini, D. Rus, and R. Grosu, "Liquid time-constant networks," in *Proceedings of the AAAI Conference on Artificial Intelligence*, vol. 35, no. 9, 2021, pp. 7657-7666.

[25] R. Hasani, M. Lechner, A. Amini, L. Liebenwein, A. Ray, M. Tschaikowski, G. Teschl, and D. Rus, "Closed-form continuous-time neural networks," *Nature Machine Intelligence*, vol. 4, no. 11, pp. 992-1003, 2022.

[26] E. Haber and L. Ruthotto, "Stable architectures for deep neural networks," *Inverse problems*, vol. 34, no. 1, p. 014004, 2017.

[27] R. T. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud, "Neural ordinary differential equations," *Advances in neural information processing systems*, vol. 31, 2018.

[28] M. Lechner, R. Hasani, M. Zimmer, T. A. Henzinger, and R. Grosu, "Designing worm-inspired neural networks for interpretable robotic control," in *2019 international conference on robotics and automation (ICRA)*. IEEE, 2019, pp. 87-94.

[29] M. Poli, S. Massaroli, J. Park, A. Yamashita, H. Asama, and J. Park, "Graph neural ordinary differential equations," *arXiv preprint arXiv:1911.07532*, 2019.

[30] Y. Qin, W. Ju, H. Wu, X. Luo, and M. Zhang, "Learning graph ode for continuous-time sequential recommendation," *IEEE Transactions on Knowledge and Data Engineering*, 2024.

[31] Y. Su, B. Ren, and K. Zhang, "Graph ode recurrent neural networks for traffic flow forecasting," in *2022 IEEE 5th International Conference on Electronics and Communication Engineering (ICECE)*. IEEE, 2022, pp. 178-182.

[32] D. I. Shuman, S. K. Narang, P. Frossard, A. Ortega, and P. Vanderheynst, "The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains," *IEEE Signal Processing Magazine*, vol. 30, no. 3, pp. 83-98, 2013.

[33] F. Gama, J. Bruna, and A. Ribeiro, "Stability properties of graph neural networks," *IEEE Transactions on Signal Processing*, vol. 68, pp. 5680-5695, 2020.

[34] Q. Li, F. Gama, A. Ribeiro, and A. Prorok, "Graph neural networks for decentralized path planning," in *Proc. International Conference on Autonomous Agents and Multiagent Systems*, 2020, pp. 1901-1903.

[35] F. Gu, H. Chang, W. Zhu, S. Sojoudi, and L. El Ghaoui, "Implicit graph neural networks," *Advances in Neural Information Processing Systems*, vol. 33, pp. 11984-11995, 2020.

[36] D. Lukovnikov, J. Lehmann, and A. Fischer, "Improving the long-range performance of gated graph neural networks," *arXiv preprint arXiv:2007.09668*, 2020.

[37] T. Ström, "On logarithmic norms," *SIAM Journal on Numerical Analysis*, vol. 12, no. 5, pp. 741-753, 1975.

[38] E. Tolstaya, F. Gama, J. Paulos, G. Pappas, V. Kumar, and A. Ribeiro, "Learning decentralized controllers for robot swarms with graph neural networks," in *Conference on robot learning*. PMLR, 2020, pp. 671-682.

[39] W. Yu, G. Chen, and M. Cao, "Distributed leader-follower flocking control for multi-agent dynamical systems with time-varying velocities," *Systems & Control Letters*, vol. 59, no. 9, pp. 543-552, 2010.

[40] L. E. Beaver and A. A. Malikopoulos, "An overview on optimal flocking," *Annual Reviews in Control*, vol. 51, pp. 88-99, 2021.

[41] H. G. Tanner, A. Jadbabaie, and G. J. Pappas, "Stable flocking of mobile agents part i: dynamic topology," in *42nd IEEE International Conference on Decision and Control (IEEE Cat. No. 03CH37475)*, vol. 2. IEEE, 2003, pp. 2016-2021.

[42] S. Ross, G. Gordon, and D. Bagnell, "A reduction of imitation learning and structured prediction to no-regret online learning," in *Proceedings of the fourteenth international conference on artificial intelligence and statistics*. JMLR Workshop and Conference Proceedings, 2011, pp. 627-635.

[43] D. P. Kingma and J. Ba, "Adam: A method for stochastic optimization," *CoRR*, vol. abs/1412.6980, 2014.