P S ¸ekercio Glu
email: pelin.sekercioglu@onera.fr.

I Sarras

J Marzat

Distributed Bipartite Containment Tracking over Signed Networks with Multiple Leaders

Pelin S ¸ekercio glu Elena Panteley Ioannis Sarras Antonio Loría Julien Marzat

Abstract-We address the distributed bipartite containment tracking-control problem for second-order systems steered by multiple cooperative and antagonistic leaders. The influence of the former is represented via positive links, while negative interconnections are used for the latter, thereby generating a signed network. Because of the presence of diverse leaders, having positive or negative influence on the followers, complete consensus is impossible. Instead, the followers' respective states may only converge to a residual compact set, not predefined, but resulting from the agents' initial conditions and the network's topology. For this set, which is called bipartite containment set, we establish global exponential stability and we compute the specific equilibria to which all agents converge inside the containment set (this is called multiconsensus). In addition, we provide strict Lyapunov functions for signed-graph networks and establish robustness of the bipartite containment tracking control. Finally, we illustrate our theoretical findings via numerical simulations. Index Terms-Multiconsensus, bipartite containment tracking, containment, autonomous vehicles, signed graphs

I. INTRODUCTION

In the study of coordination of multiagent networks, a large number of problems have been extensively addressed, such as consensus for first-order [START_REF] Olfati-Saber | Consensus problems in networks of agents with switching topology and time-delays[END_REF], second-order [START_REF] Ren | Distributed multi-vehicle coordinated control via local information exchange[END_REF] and linear high-order dynamics [START_REF] Wieland | On consensus in multi-agent systems with linear high-order agents[END_REF]. In particular, when there is a leader in the network, all followers converge to the leader's state [START_REF] Ni | Leader-following consensus of multi-agent systems under fixed and switching topologies[END_REF]. However, this so-called leader-follower consensus does not occur when the network contains more than one leader. For this case, in which multiple consensus equilibria may appear, it appears more appropriate to speak of containment control [START_REF] Cao | Distributed Coordination of Multi-agent Networks: Emergent Problems, Models, and Issues[END_REF]. This problem consists in making all the followers' states converge to a convex hull that is determined by the leaders' initial conditions and is solvable if for each follower there exists at least one leader that has a directed path to that follower.

There are multiple studies on distributed containment control, e.g., for social networks [START_REF] Kan | Containment control for a social network with state-dependent connectivity[END_REF] or for networks of singleintegrators [START_REF] Cao | Distributed containment control with multiple stationary or dynamic leaders in fixed and switching directed networks[END_REF], [START_REF] Ji | Containment control in mobile networks[END_REF], double-integrators [START_REF] Cao | Distributed containment control for multiple autonomous vehicles with double-integrator dynamics: algorithms and experiments[END_REF], [START_REF] Liu | Necessary and sufficient conditions for containment control of networked multi-agent systems[END_REF], and general linear autonomous systems [START_REF] Li | Distributed containment control of multi-agent systems with general linear dynamics in the presence of multiple leaders[END_REF]. However, most of the current research on the consensus or containment problems for multiagent systems is on cooperative networks, i.e., the coordination of nodes is achieved only by cooperative interactions. Yet, there are many scenarii in which agents may be competitive or even disruptive. These may appear, e.g., in robotics applications, in the context of herding control [START_REF] Sebastián | Adaptive multirobot implicit control of heterogeneous herds[END_REF], in aerospace applications involving control of multiple satellites that must avoid face debris represented as non-cooperative leaders, and in the context of social networks that include deceiving influencers injecting disinformation [START_REF] Altafini | Consensus problems on networks with antagonistic interactions[END_REF], to mention a few. In this case, the interconnections may be either positive or negative and at least two consensus equilibria appear-we speak of bipartite consensus [START_REF] Altafini | Consensus problems on networks with antagonistic interactions[END_REF], [START_REF] Valcher | On the consensus and bipartite consensus in high-order multi-agent dynamical systems with antagonistic interactions[END_REF]. In the case that multiple leaders (cooperative and/or antagonistic) appear, moreover, the overall behavior is even more complex [START_REF] Meng | Extended structural balance theory and method for cooperative-antagonistic networks[END_REF] and more than two consensus equilibria may appear [START_REF] Monaco | On multi-consensus and almost equitable graph partitions[END_REF]- [START_REF] Hu | Emergent collective behaviors on coopetition networks[END_REF]. Thus, signed interconnections and leader multiplicity are distinct aspects that may appear in a number of usual scenarios involving networked systems.

To analyze directed signed networks containing multiple leaders, in [START_REF] Meng | Bipartite containment tracking of signed networks[END_REF] the notion of containment control is extended to bipartite containment tracking-control. Succeeding [START_REF] Meng | Bipartite containment tracking of signed networks[END_REF], bipartite containment has been also studied in [START_REF] Zhang | Cooperative bipartite containment control for multiagent systems based on adaptive distributed observer[END_REF], [START_REF] Meng | High-order bipartite containment control in multiagent systems over time-varying cooperation-competition networks[END_REF] and [START_REF] Wang | Fixed-time bipartite containment of multi-agent systems subject to disturbance[END_REF]. In general, it consists in having all the followers' states converge to a convex hull determined by the network's topology and the initial conditions of the leaders and the mirror leaders. Moreover, in [START_REF] Meng | Bipartite containment tracking of signed networks[END_REF] and [START_REF] Wang | Fixed-time bipartite containment of multi-agent systems subject to disturbance[END_REF] limit points for the followers states are given explicitly. It is also important to mention that in [START_REF] Zhang | Cooperative bipartite containment control for multiagent systems based on adaptive distributed observer[END_REF]- [START_REF] Wang | Fixed-time bipartite containment of multi-agent systems subject to disturbance[END_REF] only structurally balanced networks are considered. Roughly speaking, these are networks whose representing graph contains positive-only cycles or none at all-see Section III and [START_REF] Zhang | Bipartite consensus of multi-agent systems over signed graphs: state feedback and output feedback control approaches[END_REF] for details. Such networks exclude important cases, such as certain agents capable of being friends with the enemies of their friends. Furthermore, in [START_REF] Zhang | Cooperative bipartite containment control for multiagent systems based on adaptive distributed observer[END_REF], only cooperative leaders are considered and in [START_REF] Wang | Fixed-time bipartite containment of multi-agent systems subject to disturbance[END_REF] interconnections between the followers are assumed to be only undirected.

As in [START_REF] Meng | Bipartite containment tracking of signed networks[END_REF], in this paper we study the bipartite containment tracking problem over structurally balanced and unbalanced signed networks with multiple cooperative or competitive leaders. Also, as in [START_REF] Meng | Bipartite containment tracking of signed networks[END_REF], we allow leaders to have neighbors. In contrast to the latter reference, however, our results apply to second-order systems, which better describe mechanical systems and a variety of (feedback linearizable) autonomous vehicles [START_REF] Tzafestas | Introduction to mobile robot control[END_REF]. In contrast to [START_REF] Zhang | Cooperative bipartite containment control for multiagent systems based on adaptive distributed observer[END_REF]- [START_REF] Wang | Fixed-time bipartite containment of multi-agent systems subject to disturbance[END_REF] our main statements hold for the general case of structurally unbalanced networks and, beyond bipartite containment tracking, we provide explicit estimates of the limit points of the followers. Furthermore, we give sufficient conditions for exponential stability of the containment set and, in contrast to all references mentioned previously, our proofs are constructive; we provide a strict Lyapunov function regardless of whether the network is structurally balanced or unbalanced. Exponential stability is important because it covers the much weaker property of mere convergence to the limit points (or to the interior of a convex hull); providing strict Lyapunov functions is significant because they are a basis to establish, in addition, input-to-state stability. Thus, our main results guarantee robustness of the containment set with respect to additive perturbations.

From a technical viewpoint, our main results are based on the framework introduced in [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF]. We recast the bipartite containment problem into one of stability of a set of the appropriately defined errors. Then, we generalize a statement in [START_REF] Panteley | Strict Lyapunov functions for consensus under directed connected graphs[END_REF] on the Lyapunov characterization of exponential stability of sets, for linear systems with one pole at zero and others having negative real part. The technical result we provide here applies to systems with several poles at the origin. Then, in addition to exponential stability of the containment set, we provide the explicit limit values of the followers by a matrix determined by all eigenvectors associated to the zero eigenvalues. In particular, when the considered graph is a traditional cooperative network, the bipartite containment tracking results stay valid for containment tracking.

The remainder of the paper is organized as follows. In the next section, we describe the bipartite containment-tracking control problem and how it may be recast as one of stability analysis. In Section III we present the technical statements described above. In Section IV we present our main results. Some numerical examples are presented in Section V, and we wrap up the paper with some closing remarks in Section VI.

II. PROBLEM FORMULATION

Consider a group of n second-order dynamical systems modeled by ẋ1i

= x 2i , x 1i , x 2i ∈ R (1a) ẋ2i = u i , u i ∈ R, i ∈ I N (1b)
where I N := {1, 2, . . . , n}.

Remark 1: The system (1) is a basic representation of (feedback-linearizable) mechanical systems; it is for notational simplicity, and without loss of generality, that we assume that x 1i ∈ R, but the contents of this paper apply to systems of higher dimension-see Section V.

• The consensus problem for [START_REF] Olfati-Saber | Consensus problems in networks of agents with switching topology and time-delays[END_REF], that is, to ensure that

lim t→∞ [x j (t) -x i (t)] = 0, ∀ i, j ≤ n, (2)
where x i = x 1i x 2i , via distributed control, is completely solved under various conditions on the interconnections and the resulting network's topology. For instance, for static directed networks, it is well-known that under the consensuscontrol law,

u i = -k 1 n j=1 a ij (x 1i -x 1j ) -k 2 n j=1 a ij (x 2i -x 2j ), (3)
where a ij ∈ R is the adjacency weight between nodes i and j, with k 1 , k 2 > 0, the expressions in (2) hold if and only if the underlying graph contains a directed spanning tree [START_REF] Ren | Distributed consensus in multi-vehicle cooperative control[END_REF]. Moreover, the consensus equilibrium may be computed exactly. Indeed, if there exists a directed spanning tree, the resulting Laplacian matrix,

L := [ ij ] ∈ R n×n , where ij = k∈I N a ik i = j -a ij i = j, (4)
has exactly one eigenvalue equal to zero. Therefore, the consensus equilibrium is uniquely determined by the left eigenvector v l corresponding to that eigenvalue, x 1m := v l x 1 (0) [START_REF] Ren | Distributed consensus in multi-vehicle cooperative control[END_REF]. Furthermore, one can construct a strict Lyapunov function to establish exponential stability of the origin in the space of the synchronization errors e := x 1 -v r x 1m [START_REF] Panteley | Strict Lyapunov functions for consensus under directed connected graphs[END_REF], where v r is the right eigenvector corresponding to zero eigenvalue. Such is the case, at least, when all nodes in the network are cooperative, that is, if a ij ≥ 0 for all i, j ≤ n. For networks in which some of the nodes are competitive, we have a ij < 0 for some i, j ≤ n. In such scenario, the distributed consensus control law (3) becomes

u i = -k 1 n j=1 |a ij |(x 1i -sgn(a ij )x 1j ) -k 2 n j=1 |a ij |(x 2i -sgn(a ij )x 2j ) -k 3 x 2i , (5)
where k 1 , k 2 > 0, and k 3 ≥ 0. If k 3 = 0, bipartite containment is reached but the system's trajectories may grow unbounded.

In the case that k 3 > 0, in addition to bipartite containment, the velocities tend to zero. Under the distributed control law [START_REF] Cao | Distributed Coordination of Multi-agent Networks: Emergent Problems, Models, and Issues[END_REF], agents on a directed signed network containing a directed spanning tree or a leader achieve bipartite consensus if and only if the underlying graph is structurally balanced [START_REF] Altafini | Consensus problems on networks with antagonistic interactions[END_REF]. If the signed network is structurally balanced, the associated Laplacian matrix has a simple zero eigenvalue [START_REF] Altafini | Consensus problems on networks with antagonistic interactions[END_REF] and its associated right eigenvector has all entries equal to ±1, which results in agents converging to the same state in modulus but different in signs. In such networks, moreover, it appears natural to encounter several competitive leader nodes that inject disinformation into the network. In this particular case, bipartite consensus can no longer be achieved due to the existence of multiple leaders and antagonistic interactions in the network. Agents can then achieve bipartite containment [START_REF] Meng | Bipartite containment tracking of signed networks[END_REF], where followers converge to the convex hull spanned by all leaders' trajectories and their mirrored trajectories, in order to avoid competitive leaders in the system.

In this paper we analyze the behavior of the networked systems (1) in closed loop with the distributed control law [START_REF] Cao | Distributed Coordination of Multi-agent Networks: Emergent Problems, Models, and Issues[END_REF] and under the assumption that several cooperative or competitive leaders interfere. Since consensus is unattainable, the control goal is set to achieving bipartite containment tracking, that is,

lim t→∞ [|x j (t)| -max i∈L |x i (t)|] ≤ 0, j ∈ F, (6)
where L and F are the sets of leaders and followers respectively-see the Standing Assumption below. This is the problem solved, e.g., in [START_REF] Meng | Bipartite containment tracking of signed networks[END_REF]- [START_REF] Wang | Fixed-time bipartite containment of multi-agent systems subject to disturbance[END_REF]. Beyond the inequality in (6), we give the explicit limit value of the followers' states and establish exponential stability of the containment set. This is significant because it allows to give a tight bound depending only on the initial conditions of the leaders. This is done under the following Standing Assumption:

1) The network contains m leaders, which can be organized into k groups of p i leaders included in a strongly connected subgraph, where 1 ≤ k ≤ m ≤ n, i ∈ I K , and

k i=1 p i = m -cf. [19, Definition 1]. 2)
Given each follower ν j , i.e., with j ∈ F with F := {m + 1, m + 2, . . . , n}, there exists at least one leader ν i i.e., with i ∈ L := {1, 2, . . . , m}, such that there exists at least one path from

ν i to ν j -cf. [19, Condition 1].
The Standing Assumption is mild; in the case that there exists no more than one leader group, it boils down to the necessary condition for consensus that imposes the existence of a spanning tree. As the networks considered here contain, a priori, multiple leaders, the resulting Laplacian matrix has as many zero eigenvalues and associated eigenvectors as the number of groups of leaders [START_REF] Caughman | Kernels of directed graph Laplacians[END_REF]. This also results in multiple convergence points for the agents. Therefore, in contrast to the consensus equilibrium x 1m = v l x 1 (0) for single-leader networks, the final states of the agents are determined by all eigenvectors associated to the zero eigenvalues. One of this paper's contributions is to show that under the control law [START_REF] Cao | Distributed Coordination of Multi-agent Networks: Emergent Problems, Models, and Issues[END_REF] and the Standing Assumption, the limit-values of the agents are given by

x m := (V ⊗ I 2 )x, (7)
where V is a matrix determined by all the eigenvectors associated to the k null eigenvalues of the Laplacian. More precisely, the matrix V is given by,

V := k i=1      v ri,1 v ri,2 . . . v ri,n      v li,1 v li,2 . . . v li,n , (8)
where v ri and v li denote, respectively, the right and left eigenvectors of the matrix L corresponding to the ith 0 eigenvalue, and, for each j ≤ n, v ri,j and v li,j denote, respectively, the jth element of the ith right and left eigenvectors. In particular, as we demonstrate further below, a possible form of the right and left eigenvectors associated with the zero eigenvalues is

v r1 =        1 p1 0 p2 . . . 0 p k ξ 1        , v r2 =        0 p1 1 p2 . . . 0 p k ξ 2        , . . . , v r k =        0 p1 0 p2 . . . 1 p k ξ k        , (9a)
v l1 =        ρ p1 0 p2 . . . 0 p k 0 n-m        , v l2 =        0 p1 ρ p2 . . . 0 p k 0 n-m        , . . . , v l k =        0 p1 0 p2 . . . ρ p k 0 n-m        , (9b)
where

ξ i ∈ R n-m is defined in Lemma 1 farther below, 1 pi ∈ R pi is a vector of ones, 0 pi ∈ R pi and 0 n-m ∈ R n-m
are vectors of zeros, and

ρ pi ∈ R pi . If p i = 1
, where i ≤ k, we have ρ pi = 1, so the corresponding left eigenvector v li has a unique non-zero element equal to 1.

If p i > 1, v li has p i elements belonging to (-1, 1), where pi l=1 |ρ l | = 1.
The definition of the eigenvectors given in [START_REF] Cao | Distributed containment control for multiple autonomous vehicles with double-integrator dynamics: algorithms and experiments[END_REF] covers that of [START_REF] Caughman | Kernels of directed graph Laplacians[END_REF]Corollary 4.2] and [16, Proposition 3], which are restricted to unsigned networks. With this under consideration, we establish bipartite containment of the system (1) and, more significantly, that x → x m -as defined in (7)-exponentially.

III. ANALYSIS APPROACH

Our main results are based on original statements for networks with an associated Laplacian having multiple zero eigenvalues. First, following the framework laid in [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF], we show how to construct the matrix V in ( 7)-( 8), to define the average states of the agents. Then, we extend the method of [START_REF] Panteley | Strict Lyapunov functions for consensus under directed connected graphs[END_REF], to construct strict Lyapunov functions for linear systems with one zero eigenvalue, to the case of multiple null eigenvalues.

To that end, we first recall certain notations and definitions that will be used later. A signed graph is said to be structurally balanced if it may be split into two disjoint sets of vertices

V 1 and V 2 , where V 1 ∪ V 2 = V, V 1 ∩ V 2 = 0 such that for every i, j ∈ V p , p ∈ [1, 2] if [a ij ] ≥ 0 while for every i ∈ V p , j ∈ V q , p, q ∈ [1, 2], p = q if [a ij ] ≤ 0. It is structurally unbalanced, otherwise.
The following definition is also recalled from [START_REF] Caughman | Kernels of directed graph Laplacians[END_REF] to introduce some useful sets of vertices in a graph. Definition 1: A set R of vertices is called a reach if it is a maximal reachable set that consists of a leader group and its followers. For each reach R i of a graph, we define the exclusive part of R i as the set H i = R i \ j =i R j , i.e., the set of followers influenced only by the leader group i, and the common part of R i as the set C i = R i \H i , i.e., the set of followers influenced by leaders other than the ith leader group.

The following statement, which leads to the construction of the matrix V, is an original contribution of this paper and extends Corollary 4.2 of [START_REF] Caughman | Kernels of directed graph Laplacians[END_REF] to the case of signed networks.

Lemma 1: Let G denote a directed signed graph and let L denote the associated Laplacian matrix. Suppose G has n vertices and k reaches. Then, the algebraic and geometric multiplicity of the eigenvalue 0 is equal to k.

Furthermore, if G is structurally unbalanced, the eigenspace generated by the eigenvectors associated to the null eigenvalues, has a basis defined by the vectors γ i ∈ R n , with i ≤ k, whose elements satisfy the following:

1) γ i,j = 0 for j / ∈ R i , 2) |γ i,j | = 1 for j ∈ H i , 3) |γ i,j | ≤ 1 for j ∈ C i , 4) j |γ i | ≤ 1 n . On the other hand, in the more restrictive case that G is structurally balanced, 1') γ i,j = 0 for j / ∈ R i 2') γ i,j = 1, if (ν j , ν i ) ∈ V 1 -1, if ν j ∈ V 1 , ν i ∈ V 2 for j ∈ H i 3') γ i,j ∈ (0, 1), if (ν j , ν i ) ∈ V 1 (-1, 0), if ν j ∈ V 1 , ν i ∈ V 2 for j ∈ C i 4') j |γ i | = 1 n
, where V 1 and V 2 are the two disjoint sets of vertices, i ∈ I K , j ∈ I N , and γ i,j denotes the jth element of γ i .

Proof: Let L + := D -|A| denote the associated Laplacian matrix of G + , where D is the in-degree matrix and A is the adjacency matrix of G. G + and G have identical reaches. Then, the fact that the algebraic and geometric multiplicity of 0 equals k follows from [START_REF] Caughman | Kernels of directed graph Laplacians[END_REF]Theorem 3.2]. Now, for structurally-balanced graphs, Items 1) and 2), follow from Definition 1, under the Standing Assumption. Item 3) follows from [START_REF] Meng | Extended structural balance theory and method for cooperative-antagonistic networks[END_REF], while Item 4) results from computing the null space of L, which is generated by γ i ∈ R n such that Lγ i = 0 for each i ≤ k. This gives L k i=1 γ i = 0. Then, under the Standing Assumption, the Laplacian matrix being equal to

L = 0 0 -A lf L f + ∆ |A lf | ,
where [START_REF] Cao | Distributed containment control for multiple autonomous vehicles with double-integrator dynamics: algorithms and experiments[END_REF]. The sum of the remaining rows of γ i s gives the following, [START_REF] Meng | Bipartite containment tracking of signed networks[END_REF]. Thus, we conclude that

A lf ∈ R (n-m)×m and L f + ∆ |A lf | ∈ R (n-m)×(n-m) , we have -A lf k i=1 γ im + (L f + ∆ |A lf | ) k i=1 γ in-m = 0. This gives us (L f + ∆ |A lf | ) k i=1 γ in-m = A lf k i=1 γ im , where k i=1 γ im = 1 m from
k i=1 γ in-m = (L nr + ∆ |Arnr| ) -1 A rnr 1 m and it follows that j |γ in-m | ≤ 1 n-m from Lemma 6 of
γ i,j = 0 for j / ∈ R i , |γ i,j | = 1 for j ∈ H i , |γ i,j | ≤ 1 for j ∈ C i , and j |γ i | ≤ 1 n .
On the other hand, if the graph is structurally balanced, we may apply the gauge transformation, which consists in a change of coordinates performed by the matrix D = diag(σ), where σ = [σ 1 , ..., σ n ], σ j ∈ {1, -1}, j ∈ I N [START_REF] Altafini | Consensus problems on networks with antagonistic interactions[END_REF]. Let , we obtain the following for a structurally balanced signed network: γ i,j = σ j γ Di,j = 0 for j / ∈ R i , γ i,j = σ j γ Di,j = σ j for i ∈ H i , γ i,j = σ j γ Di,j ∈ σ j (0, 1) for j ∈ C i and j γ i = j Dγ Di = D1 n , which gives j |γ i | = 1 n . Items 1') to 4') follow. From Lemma 1, setting v ri,j = γ i,j , we obtain the form given in (9a) for the k right eigenvectors associated to the zero eigenvalues, such that for each j ≤ m, v ri,j is either equal to 1 or to 0, because the jth leader can only be in the exclusive part of its corresponding reach. The remaining rows ξ i,j of v ri belong, either to {-1, 1} or to (-1, 1), depending on network's topology and signs of the interconnections. Moreover, under the Standing Assumption and the given form in (9a) for right eigenvectors, as the Laplacian matrix has all entries equal to 0 for its first m rows, we obtain (9b) for k left eigenvectors associated to zero eigenvalues. Notice that, because of the form of the left eigenvectors, each column of the matrix V has the same properties as the basis defined in Lemma 1. Therefore, we may split V in four blocks, as follows:

V = k i=1 v ri v li = V l 0 m×(n-m) V f 0 (n-m)×(n-m) , (10)
where V l ∈ R m×m represents leaders' interactions and V f ∈ R (n-m)×m represents leader-follower interactions.

It is important to remark that the elements of V f have the properties stated in Lemma 1. This is significant because, owing to the fact that x → x m where x m is defined by [START_REF] Cao | Distributed containment control with multiple stationary or dynamic leaders in fixed and switching directed networks[END_REF], it is clear that V f is the matrix that defines the limit point for the followers. Notice that even if during their trajectories, depending on network's topology, followers' states may be influenced by other followers' states, it follows from (10) that the final states of the followers are defined only by leaders' states. Now, akin to the case of networks containing one leader, in which case the error is defined as e = x 1 -v r x 1m , with x 1m := v l x 1 , for multi-leader networks we define the consensus errors as

e := ([I -V] ⊗ I 2 )x, (11)
where e = e 1 e 2 . Then, to establish beyond the convergence statements that x 1 → x 1m and x 2 → x 2m and, consequently, the bipartite containment objectives defined by (6), we will prove the stronger property of global exponential stability of the set {(e 1 , e 2 ) = (0, 0)}. For that, we shall show how to construct strict-in the space of (e 1 , e 2 )-Lyapunov functions, based on the following proposition, which is another original contribution of this paper and extends Proposition 1 of [START_REF] Panteley | Strict Lyapunov functions for consensus under directed connected graphs[END_REF] to the case of signed networks with multiple leaders. Proposition 1: Let G be a directed signed network with multiple leaders. Then the following are equivalent:

(i) the graph has k groups of leaders and m leaders, and given each follower ν j , ∀j ∈ F, there exists at least one leader ν i , ∀i ∈ L, such that there exists at least one path from ν i to ν j , (ii) for any Q ∈ R n×n , Q = Q > 0 and for any α i > 0, there exists a matrix P (α i ) ∈ R n×n , P = P > 0 such that

P L + L P = Q - k i=1 α i (P v ri v li + v li v ri P ), (12)
where v ri , v li ∈ R are the right and left eigenvectors of L associated with the ith 0 eigenvalue. Remark 2: Note that the Proposition provides a Lyapunov characterization of the second part of the Standing Assumption.

• Proof: (i) ⇒ (ii): By assumption, the graph G has k leader groups and is connected. Then, from Lemma 1, it follows that L has k zero eigenvalues, and the rest of its eigenvalues have positive real parts:

0 = λ 1 = • • • = λ k < e(λ k+1 ) ≤ • • • ≤ e(λ n ).
Following the lines of proof as for Lemma 2 in [START_REF] Panteley | Strict Lyapunov functions for consensus under directed connected graphs[END_REF], we can write the Jordan decomposition of L as

L = U ΛU -1 = k i=1 λ i (L)v ri v li + U 1 Λ 1 U † 1 with Λ 1 ∈ C n-k×n-k , U = v r1 . . . v r k U 1 ∈ C n×n and U -1 = v l1 . . . v l k U † 1 ∈ C n×n . For any α i > 0 define R(α i ) = L + k i=1 α i v ri v li .
From this decomposition and the properties of Λ 1 , e{λ j (R)} > 0 for all j ≤ n.

-R(α i ) is Hurwitz, therefore for any Q = Q > 0 and α i > 0, i ≤ m, there exists P = P > 0 such that

-P R(α i ) -R(α i ) P = -Q, -P (L + k i=1 α i v ri v li ) -(L + k i=1 α i v ri v li ) P = -Q, P L + L P = Q - k i=1 α i (P v ri v li + v li v ri P ).
(ii) ⇒ (i): Let statement (ii) hold and assume that the Laplacian matrix has k + 1 zero eigenvalues and the rest of its eigenvalues have positive real parts. In view of Lemma 1, the assumption that the system has k groups of leaders does not hold. Now, the Jordan decomposition of L has the form

L = U ΛU -1 = k+1 i=1 λ i (L)v ri v li + U 1 Λ 1 U † 1 with U = v r1 . . . v r k+1 U 1 and U -1 = v l1 . . . v l k+1 U † 1 . Next let us consider R(α i ) = L + k i=1 α i v ri v li which admits the Jordan decomposition R := U Λ R U -1 , where Λ R :=        α 1 . . . α k 0 Λ 1        .
Clearly, R is not positive definite because one of its eigenvalues is equal to zero. Then, there exists a matrix Q = Q for which there does not exist a matrix P = P such that -P R -R P = -Q, which contradicts statement (ii).

IV. BIPARTITE CONTAINMENT OF SECOND-ORDER SYSTEMS

In this section, we present our main results. We will consider bipartite containment protocols for second-order systems, with and without absolute velocity damping. Furthermore, we establish robustness of the bipartite containment tracking in the sense of input-to-state stability with respect to external bounded perturbations.

A. Exponential stability

Consider the system (1), interconnected with the bipartite containment control law [START_REF] Cao | Distributed Coordination of Multi-agent Networks: Emergent Problems, Models, and Issues[END_REF]. We analyze the dynamics of the errors in [START_REF] Li | Distributed containment control of multi-agent systems with general linear dynamics in the presence of multiple leaders[END_REF]. Differentiating the latter on both sides, to obtain

ė = ([I -V] ⊗ I 2 ) ẋ
and using ( 1) and ( 5), we obtain the closed-loop dynamical equations

ė1 = e 2 (13a) ė2 = -k 1 Le 1 -k 2 Le 2 -k 3 e 2 . ( 13b
)
Guaranteeing the bipartite containment problem is now recast as a problem of stability analysis of the dynamical system [START_REF] Altafini | Consensus problems on networks with antagonistic interactions[END_REF]. Thus, relying on Proposition 1, our next statement provides sufficient conditions on the controller gains to achieve global exponential stability of the set {(e 1 , e 2 ) = (0, 0)}, which covers the bipartite containment tracking objective (6). Proposition 2: Consider the system (13) and let P be generated by [START_REF] Sebastián | Adaptive multirobot implicit control of heterogeneous herds[END_REF] with Q = I n . Then, under the Standing Assumption, the set {(e 1 , e 2 ) = (0, 0)} is exponentially stable if

k 1 > 0, ( 14a
)
k 2 > 2 k 1 λ P , (14b) k 2 k 3 ≥ λ P , k 3 ≥ 0, (14c)
where λ P ≥ |P | is the largest eigenvalue of P . Remark 3: For conciseness, we consider in one statement two interesting and distinct cases, in which k 3 > 0 and in which k 3 = 0. Even though the value of k 3 ≥ 0 is inconsequential for exponential stability statement to hold, it plays a role on performance and the solutions' behavior. Indeed, k 3 > 0 implies the injection of velocity damping, hence, it guarantees that the agents states x 1i converge to constant values, while x 2i → 0. In the case that k 3 = 0, containment is achieved in both variables x 1i and x 2i but the solutions grow unbounded. In the particular case of vehicles moving on a plane, these cases are tantamount to the vehicles converging to a rendezvous zone or permanently moving in formation and in containment, steered by the leaders. See Proposition 3 further below.

• Proof of Proposition 2: Let Q = Q > 0 and α > 0 be arbitrarily fixed. Since by the Standing Assumption and Proposition 1, ∃P = P > 0 such that (12) holds. Then, consider the following Lyapunov function candidate

V (e) = 1 2 |e 1 | 2 + e 1 P e 2 + µe 2 P e 2 , ( 15
)
which is positive definite (for all e as in ( 11)), under the condition ≤ 2µ |P | , µ > 0. The total time derivative of V along the trajectories yields

V (e) = -k 1 e 1 P Le 1 -e 2 [k 2 µ(P L + L P ) + 2k 3 µP -P ]e 2 + e 1 (I -k 2 P L -2k 1 µL P -k 3 P )e 2 .
Let µ > 0 be such that 2k 1 µ = k 2 = 1. Using [START_REF] Sebastián | Adaptive multirobot implicit control of heterogeneous herds[END_REF] with Q = I n and the identity P Hence, in compact form, we have

V (e) = - 1 2 e 1 2 k 1 I k 3 P k 3 P 4k 3 µP e - 1 4 k 1 |e 1 | 2 -e 2 [k 2 µI -P ]e 2 , ( 16
)
so V is negative definite if the matrix in ( 16) is positive semidefinite and [k 2 µI -P ] is positive definite, which, using 2k 1 µ = k 2 = 1 we see that it is so if and only if

k 2 I -2 k 1 k 2 λ P I > 0,
which, in turn, holds if k 2 > 2k 1 λ P , that is, under condition (14b). Then, computing the Schur complement, the condition for the matrix in ( 16) to be positive semi-definite, is 4k 3 µP -(k 3 P ) 2 k1 k 3 P ≥ 0. If k 3 = 0 the latter holds trivially. Otherwise, considering k 3 > 0, and using 2 k1 k2 = µ , we see that the latter inequality is equivalent to

2 µ µ I - k 3 k 1 µ P = 2 µ µ I - k 3 k 1 2k 1 k 2 P = 1 - k 3 k 2 λ P ≥ 0,
which is satisfied under (14c). We conclude that

V (e) ≤ - 1 4 k 1 |e 1 | 2 + [k 2 µ -λ P ]|e 2 | 2 , ( 17
)
so the statement of the proposition follows.

The following statement provides explicit expressions for the limit values of the followers' states and emphasizes the role of the gain k 3 .

Proposition 3: Consider the system (1) and the bipartite containment control law [START_REF] Cao | Distributed Coordination of Multi-agent Networks: Emergent Problems, Models, and Issues[END_REF]. Under the Standing Assumptions, the bipartite containment objective is achieved, that is, the inequalities (6) hold. Furthermore, if k 3 > 0, the final states of the followers satisfy

lim t→∞ x 1 f (t) = V f x 1 l (0) + 1 k 3 V f x 2 l (0), ( 18a
)
lim t→∞ x 2 f (t) = 0. ( 18b
)
On the other hand, if

k 3 = 0, lim t→∞ x 1 f (t) = V f x 1 l (0) + tV f x 2 l , ( 19a
)
lim t→∞ x 2 f (t) = V f x 2 l , (19b)
where x l and x f are leaders' and followers' states respectively and V f ∈ R (n-m)×m is given in [START_REF] Liu | Necessary and sufficient conditions for containment control of networked multi-agent systems[END_REF]. Proof: Differentiating the weighted average of the system (7), we obtain the dynamical equations below

ẋ1m = V ẋ1 = Vx 2 = x 2m , (20a) ẋ2m = V ẋ2 = V(-k 1 Lx 1 -k 2 Lx 2 -k 3 x 2 ) = -k 3 x 2m , (20b)
with v li L = 0 for i ≤ k. Their solutions give the following

x 2m (t) = x 2m (0)e -k3t , (21a)
x 1m (t) = x 1m (0) + x 2m (0) t 0 e -k3s ds. (21b)
From Proposition 2, we have lim t→∞ e(t) = 0, which results in lim t→∞ x(t) = x m (t). From ( 21), we obtain lim t→∞ x 2 (t) = 0 and lim t→∞ x 1 (t) = x 1m (0) + 1 k3 x 2m (0). Then, using [START_REF] Liu | Necessary and sufficient conditions for containment control of networked multi-agent systems[END_REF], we obtain the relations in [START_REF] Hu | Emergent collective behaviors on coopetition networks[END_REF]. Under the Standing Assumption and from items 4) and 4') of Lemma 1, we can write

lim t→∞ |x 1 f j (t)| = m i=1 |v i,m+j x 1 li (0) + 1 k 3 v i,m+j x 2 li (0)| ≤ m i=1 |v i,m+j ||x 1 li (0) + 1 k 3 x 2 li (0)| ≤ m i=1 |v i,m+j ||x 1 li (t)| ≤ max 1≤i≤m |x 1 li (t)|,
where v i,j are the elements of the matrix V in [START_REF] Liu | Necessary and sufficient conditions for containment control of networked multi-agent systems[END_REF]. In the case when k 3 = 0, (20) becomes

ẋ1m = x 2m , ẋ2m = 0, so x 2m (t) = x 2m (0), (22a) x 1m (t) = x 1m (0) + tx 2m (0). (22b)
From [START_REF] Wang | Fixed-time bipartite containment of multi-agent systems subject to disturbance[END_REF], we obtain lim t→∞ x 2 (t) = x 2m (0) and lim t→∞ x 1 (t) = x 1m (0) + x 2m (0)t. Then using [START_REF] Liu | Necessary and sufficient conditions for containment control of networked multi-agent systems[END_REF], we obtain the relations in [START_REF] Meng | Bipartite containment tracking of signed networks[END_REF]. Under the Standing Assumption and from items 4) and 4') of Lemma 1, we have

lim t→∞ |x 2 f j (t)| = m i=1 |v i,m+j x 2 li | ≤ m i=1 |v i,m+j | max 1≤i≤m |x 2 li | ≤ max 1≤i≤m |x 2 li |, lim t→∞ |x 1 f j (t)| = m i=1 |v i,m+j x 1 li (0) + tv i,m+j x 2 li | ≤ max 1≤i≤m |x 1 li (0)| + t max 1≤i≤m |x 2 li | ≤ max 1≤i≤m |x 1 li (t)|, so (6) follows.
It is worth noting that for unsigned networks, the achievable objective is containment [START_REF] Cao | Distributed Coordination of Multi-agent Networks: Emergent Problems, Models, and Issues[END_REF], that is,

lim t→∞ [x j (t) -max i∈L x i (t)][x j (t) -min i∈L x i (t)] ≤ 0. (23)
Then, for such networks, we recover the following statement from Proposition 3. Corollary 1: Consider the system (1) with the containment control law [START_REF] Cao | Distributed Coordination of Multi-agent Networks: Emergent Problems, Models, and Issues[END_REF] and let the associated digraph be unsigned. Under the Standing Assumption, the containment objective is achieved, that is the inequality (23) holds.

For structurally balanced signed networks, the bipartite containment set is refined so that the achievable objective is

lim t→∞ [x j (t) -max i∈L (σ i σ j x i (t))][x j (t) -min i∈L (σ i σ j x i (t))] ≤ 0, (24)
where j ∈ F and σ = [σ 1 , σ 2 , . . . , σ n ] with

σ i = 1 if (i, j) ∈ V 1 or σ i = -1 if i ∈ V 1 , j ∈ V 2 . See Section V- A for an example.
In the next section, we use the strict Lyapunov functions provided in this section in order to conduct a robustness analysis of the control law [START_REF] Cao | Distributed Coordination of Multi-agent Networks: Emergent Problems, Models, and Issues[END_REF], in the sense of input-to-state stability of the bipartite containment tracking.

B. Robustness Analysis

Consider the perturbed second-order systems

ẋ1i = x 2i (25a) ẋ2i = u i + d i (t), (25b)
where the disturbances d i : R ≥0 → R n are assumed to be essentially bounded locally integrable functions. Under the action of the control law (5), the system (25) becomes

ẋ1 = x 2 (26a) ẋ2 = -k 1 Lx 1 -k 2 Lx 2 -k 3 x 2 + d(t), (26b)
where

d := [d i ] ∈ R n .
Differentiating the errors in (11) on both sides and using [START_REF] Panteley | Strict Lyapunov functions for consensus under directed connected graphs[END_REF] we obtain

ė1 = e 2 (27a) ė2 = -k 1 Le 1 -k 2 Le 2 -k 3 e 2 + [I -V]d(t). (27b)
Then, we have the following. Proposition 4: The closed-loop system [START_REF] Ren | Distributed consensus in multi-vehicle cooperative control[END_REF], under the Standing Assumption, is input-to-state stable with respect to essentially bounded, locally integrable external disturbances if the conditions in ( 14) hold.

Proof: Consider the Lyapunov function [START_REF] Meng | Extended structural balance theory and method for cooperative-antagonistic networks[END_REF] which is positive definite, under the condition ≤ 2µ |P | , µ > 0. The total time derivative of V along the trajectories of (27) yields

V (e) = ∂V ∂e 1 e 2 + ∂V ∂e 2 (-k 1 Le 1 -k 2 Le 2 -k 3 e 2 ) + ∂V ∂e 2 [I -V]d. (28)
Then, from [START_REF] Meng | Interval bipartite consensus of networked agents associated with signed digraphs[END_REF], we obtain

V (e) ≤ - 1 4 [ k 1 |e 1 | 2 + (k 2 µ -λ P )|e 2 | 2 ] + ∂V ∂e 2 |[I -V]||d| ≤ - 1 4 k 1 |e 1 | 2 - 1 4 (k 2 µ -λ P )|e 2 | 2 + λ P e 1 |[I -V]|d + 2µλ P e 2 |[I -V]||d|.
We know that 0 ≤ |[I -V]| ≤ |I| + |V| ≤ 2, because all eigenvalues of I are equal to 1 and all eigenvalues of |V| are either 1 or 0. Let δ > 0 such that c 1 := 1 4 ( k 1 -4 δ λ P ) > 0 and c 2 := 1 4 (k 2 µ -λ P -8 δ µλ P ) > 0. Then,

V (e) ≤ -c 1 |e 1 | 2 -c 2 |e 2 | 2 + c 3 |d| 2 ,
with c 3 = λ P δ( + 2µ) > 0. Thus, the system ( 27) is inputto-state stable with respect to the bounded external input.

Corollary 2: The system [START_REF] Ren | Distributed consensus in multi-vehicle cooperative control[END_REF], under the Standing Assumption and over an unsigned network, is input-to-state stable with respect to an essentially bounded, locally integrable external disturbance if the conditions in ( 14) hold.

V. SIMULATION RESULTS

We provide some numerical examples on a system of nonholonomic unicycle mobile robots over both structurally balanced and unbalanced networks. Let r i = [r xi r yi ] be the inertial position, θ i the orientation, s i the linear speed, ω i the angular speed, m i the mass, J i the moment of inertia, F i the applied force and τ i the applied torque. Each robot has the dynamical equation

      ṙxi ṙyi θ ṡi ω       =         s i cos θ i s i sin θ i ω i 0 0 0         +       0 0 0 0 0 0 1 mi 0 0 1 Ji       η i , (29)
where η i := F i τ i . We choose a reference point p i = r i + δ i cos θ i sin θ i located at a distance δ i = 0.1m along the line that is perpendicular to the wheels' axis and we define

ζ :=       r xi + δ i cos θ i r yi + δ i sin θ i s i cos θ i -δ i ω i sin θ i s i sin θ i + δ i ω i cos θ i θ i       . (30)
In transformed coordinates, with

p i = [ζ 1i ζ 2i ] , we have ζ1i ζ2i = ζ 3i ζ 4i , ( 31a
) ζ3i ζ4i = -s i ω i sin θ i -L i ω 2 i cos θ i s i ω i cos θ i -L i ω 2 i sin θ i + 1 mi cos θ i -Li Ji sin θ i 1 mi sin θ i Li Ji cos(θ i ) η i , ( 31b
) ζ5i = - 1 2δ i ζ 3i sin ζ 5i + 1 2δ i ζ 4 I cos ζ 5i . (31c)
The feedback linearizing control η i is given by

η i = 1 mi cos θ i -Li Ji sin θ i 1 mi sin θ i Li Ji cos(θ i ) -1 × u i - -s i ω i sin θ i -L i ω 2 i cos θ i s i ω i cos θ i -L i ω 2 i sin θ i , (32)
which gives

[ ζ1i ζ2i ] = [ζ 3i ζ 4i
] and [ ζ3i ζ4i ] = u i . Thus, we implemented (32) with u i as in [START_REF] Cao | Distributed Coordination of Multi-agent Networks: Emergent Problems, Models, and Issues[END_REF],

x 1i = [ζ 1i ζ 2i ] , x 2i = [ζ 3i ζ 4i ]
, m i = 8kg and J i = 0.12kg/m 2 for each i ∈ I N .

A. Structurally Balanced Networks

Consider a network with five leaders ν i , i ≤ 5, organized in three leader groups {ν 1 , ν 2 , ν 3 }, {ν 4 }, {ν 5 }, and four followers ν j , 6 ≤ j ≤ 9, communicating over a directed graph as the one depicted in Figure 1, below. The Laplacian matrix corresponding to the graph is

ν 1 ν 3 ν 6 ν 8 ν 2 ν 4 ν 9 ν 5 ν 7 1 1 1 3 1 -1 1 -1 5 -3
L =               1 0 -1 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -3 0 0 0 0 4 0 1 0 0 0 0 0 0 -1 2 0 1 0 0 0 -5 0 1 0 6 0 0 0 0 0 3 0 -1 0 4              
and its eigenvalues are λ L = 0, 0, 0, 1.5 ± 0.86i, 2, 3.59, 4, 6.41. The network may be bipartitioned into two subgroups, V 1 = {ν 1 , ν 2 , ν 3 , ν 5 , ν 6 , ν 7 }, V 2 = {ν 4 , ν 8 , ν 9 }, so is structurally balanced, and

V l =      
0.33 0.33 0.33 0 0 0.33 0.33 0.33 0 0 0.33 0.33 0.33 0 0 0 0 0 1 0 0 0 0 0 1 We notice that V l and V f have the properties stated on Items 1 ) -4 ) of Lemma 1, as the network is structurally balanced. For V l , we have two elements equal to one, corresponding to isolated leaders ν 4 and ν 5 , as they are in the exclusive part of the reaches. Leaders ν 1 , ν 2 and ν 3 are in the exclusive part of a reach so the corresponding elements of the right eigenvector are equal to one. Thus, as they are interconnected within a strongly-connected graph, the absolute value of the elements of the left eigenvector are less than one. V f does not have an element equal to ±1 since all the followers are influenced by more than one leader. Followers ν 6 and ν 8 are not influenced by the leader ν 5 , so the corresponding elements on the fifth column are equal to zero. We also remark that the sum of the absolute value of the terms on each row is equal to 1. Now, let P be generated by [START_REF] Sebastián | Adaptive multirobot implicit control of heterogeneous herds[END_REF] with Q = I n and α = 20, then we obtain λ P = 0.5193. Consider the system (29) and the bipartite containment law (32) with k 1 = 0.8, k 2 = 1.5, and k 3 = 1.2, which satisfy the conditions in [START_REF] Valcher | On the consensus and bipartite consensus in high-order multi-agent dynamical systems with antagonistic interactions[END_REF]. The respective agents' initial states are r x (0) = [3.5 1 -0.5 -1.6 1 -6. The simulation results are showed in Figures 2 and3. The followers converge to the bipartite containment set spanned by cooperative leaders' final states and competitive leaders' mirrored final states and all agents' velocities converge to zero. Using the relations in [START_REF] Hu | Emergent collective behaviors on coopetition networks[END_REF] and the coordinate transformation, we obtain the following limit values for followers' states, lim t→∞ r x (t) = [1.99 In a second run of simulations, we tested the bipartitecontainment control law (32) on the system (29), using the same initial conditions and controller gains as above. Therefore, in the simulations we consider that the system is subject to a disturbance d i = σ i (t) 1 1 , where σ i (t) is given by

      , (33)
V f =     0
σ i (t) =    tanh(t -10) -1 + 1 (t+10) i ∈ {7, 9} -tanh(t -10) + 1 -1 (t+10) i = 6 0 i ∈ {1, 2, 3, 4, 5, 8}.
The simulation results can be appreciated in Figures 4 and5. During the first 10s, the perturbation d(t) prevents the achievement of bipartite containment tracking but, after the perturbation vanishes, the trajectories of the followers towards the bipartite containment set spanned by cooperative leaders final states and competitive leaders' mirrored final states and all agents' velocities converge to zero. We obtain the same limit values as before for the followers.

B. Structurally Unbalanced Networks

Consider a network containing five leaders ν i , i ≤ 5, three leader groups {ν 1 , ν 2 , ν 3 }, {ν 4 }, {ν 5 }, and four followers ν j , 6 ≤ j ≤ 9, provided in Figure 6. The Laplacian matrix corresponding to the graph is

L =               1
0 -1 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -3 0 0 0 0 4 0 -1 0 0 0 0 0 0 -1 2 0 1 0 0 0 -5 0 -1 0 6 0 0 0 0 0 3 0 -1 0 4

             
and its eigenvalues are λ L = 0, 0, 0, 1.5 ± 0.86i, 2, 3.59, 4, 6.41. The network may not be bipartitioned into two subgroups, so is structurally unbalanced. The matrix V l is calculated as in (33) and V f is calculated as We notice that V f has the properties stated on Items 1) -4) of Lemma 1, since the network is structurally unbalanced. Since all followers are influenced by more than one leader, there are no agents corresponding to the exclusive part of a reach and as a result V f does not have an element equal to ±1 but the absolute value of each element is less than one. We also remark that the sum of the absolute value of the terms on each row is less than 1.

V f =     0
Now, let P be generated by [START_REF] Sebastián | Adaptive multirobot implicit control of heterogeneous herds[END_REF] with Q = I n and α = 20, then we obtain λ P = 0.4839. Consider the system (29) and the bipartite containment law (32) with the same controller gains as before. The respective agents' inertial positions and linear speeds are the same as before, while the orientations are θ(0) = [0.25 0.61 -1.27 0.32 1.42 0.46 0.98 1. 19 1.03] . The simulation results are depicted in Figures 7 and8.

Note that all followers converge to the convex hull spanned by leaders' final states and mirrored final states and all agents' velocities converge to zero. Using relations in [START_REF] Hu | Emergent collective behaviors on coopetition networks[END_REF] and the coordinate transformation, we obtain the following limit values for followers' states, lim t→∞ r x (t) = [1.5 1.1 0.01 -0.94] , lim t→∞ r y (t) = [1.79 -0.01 0.89 1.54] , lim t→∞ s(t) = [0 0 0 0] .

We remark that for structurally balanced signed networks, the followers' final states converge to the bipartite containment set spanned by cooperative leaders' final states and antagonistic leaders' mirrored states. However, in the case of a structurally unbalanced network, we cannot observe this as we cannot bipartition the agents into two disjoint subsets.

VI. CONCLUSIONS

We presented a Lyapunov approach to analyze the exponential stability of the bipartite containment tracking problem of double-integrators over multiple-leader signed networks. Via a change of coordinates, we have shown the existence and characterized a bound for the convergence of the followers. Moreover, we have generalized the Lyapunov equation characterization of the Hurwitz property of a matrix to matrices having more than one zero eigenvalue, which allowed us to construct strict Lyapunov functions. Disposing of strict Lyapunov functions allowed us to establish the robustness of the system with a bounded disturbance. Even though, paper, we only address bipartite containment tracking of second-order systems, the stability and robustness properties may serve for more complex systems.

  L D = DLD denote the unsigned Laplacian matrix of the transformed network. Then, we may express the Laplacian L in Jordan canonical form, as L = DL D D = DP D ΛP -1 D D, where P D = v D1 . . . v D k P D1 ∈ C n×n , P -1 D = w D1 . . . w D k P † D1 ∈ C n×n , and Λ ∈ C n×n , with v Di and w Di , i ≤ m are the right and left eigenvectors associated to k zero eigenvalues of L D . From the Jordan decomposition, we can see that the basis of the null space of L is given by LDγ Di = LDv Di = 0, i ≤ k and has a basis defined by the columns of γ = Dγ D , where γ D = [γ D1 • • • γ D k ] and the columns {γ Di } constitute the basis of the associated eigenspace of L D . Then, using [28, Corollary 4.2]

k

  i=1 v ri v li e 1 = P k i=1 v ri v li e 2 = 0 we obtain the following:-k 1 e 1 P Le 1 = -k 1 e 1 (P L + L P )e 1 = -1 2 k 1 |e 1 | 2 , -e 2 µk 2 (P L + L P )e 2 = -µk 2 e 2 e 2 , e 1 [I -(L P + P L)]e 2 = -e 1 [I -I]e 2 = 0.

Fig. 1 .

 1 Fig. 1. Network 1: A network of nine mobile robots

  5 5.5 -3.5 6] , r y (0) = [3 2.5 -1.5 1.3 -3 -1 -3 -3 -2.5] , θ(0) = [0.25 0.61 -1.27 0.32 1.42 0.46 0.98 1.19 1.03] and s(0) = [1.2 1.2 -1.4 -1.6 1.3 -2.2 0.4 -1.1 1.2] .

Fig. 2 .

 2 Fig. 2. Bipartite containment tracking of (29) with the control (5) on the plane. The filled dots are the final states of the agents. The diamonds represent the mirrored final states of the leaders.

Fig. 3 .

 3 Fig.3. Bipartite containment tracking of (29) with (5) on velocity.

Fig. 4 .

 4 Fig. 4. Bipartite containment tracking of (25) on the plane.

Fig. 5 .

 5 Fig. 5. Bipartite containment tracking of (25) on velocity.

Fig. 6 .

 6 Fig. 6. Network 2: A network of nine mobile robots

Fig. 7 .

 7 Fig. 7. Bipartite containment tracking of (1) on position. The filled dots are the final states of the agents. The diamonds represent the mirrored final states of the leaders.

Fig. 8 .

 8 Fig. 8. Bipartite containment tracking of (1) on velocity.


Rinske Alkemade

Emanuele Boattini

Laura Filion

Frank Smallenburg

Comparing machine learning techniques for predicting glassy dynamics

Comparing machine learning techniques for predicting glassy dynamics Rinske M. Alkemade, [START_REF] Berthier | Theoretical perspective on the glass transition and amorphous materials[END_REF] Emanuele Boattini, 1 Laura Filion, [START_REF] Berthier | Theoretical perspective on the glass transition and amorphous materials[END_REF] and Frank Smallenburg 2 1) Soft Condensed Matter, Debye Institute of Nanomaterials Science, Utrecht University, Utrecht, Netherlands 2) Université Paris-Saclay, CNRS, Laboratoire de Physique des Solides, 91405 Orsay, France

In the quest to understand how structure and dynamics are connected in glasses, a number of machine learning based methods have been developed that predict dynamics in supercooled liquids. These methods include both increasingly complex machine learning techniques, and increasingly sophisticated descriptors used to describe the environment around particles. In many cases, both the chosen machine learning technique and choice of structural descriptors are varied simultaneously, making it hard to quantitatively compare the performance of different machine learning approaches. Here, we use three different machine learning algorithms -linear regression, neural networks, and GNNs -to predict the dynamic propensity of a glassy binary hard-sphere mixture using as structural input a recursive set of order parameters recently introduced by Boattini et al. [Phys. Rev. Lett. 127, 088007 (2021)]. As we show, when these advanced descriptors are used, all three methods predict the dynamics with nearly equal accuracy. However, the linear regression is orders of magnitude faster to train making it by far the method of choice.

The relationship between local structure and dynamics in glassy systems has been a heavily debated question in condensed matter for several decades [START_REF] Berthier | Theoretical perspective on the glass transition and amorphous materials[END_REF][START_REF] Royall | The role of local structure in dynamical arrest[END_REF][START_REF] Tanaka | Revealing key structural features hidden in liquids and glasses[END_REF] . Over the last few years, one avenue for exploring this relationship has been an effort to predict dynamical behavior based on local structural features using various machine learning (ML) methods. Pioneered by Cubuk et al. [START_REF] Cubuk | Identifying structural flow defects in disordered solids using machine-learning methods[END_REF] using support vector machines (SVMs) to predict rearrangement probabilities in glassy mixtures, this area of research has now embraced a wide variety of ML techniques, including e.g. linear regression, convolutional neural networks (CNNs), graph neural networks (GNNs), autoencoders, and community inference, see e.g. Refs. 5-10. This raises the question of what ML technique one should choose when predicting the dynamics of a glassy system.

This question is far from straightforward, since in addition to choosing a machine learning technique, one also has to make a choice with respect to the encoding of the local structure in terms of data that can be interpreted by an ML algorithm. For most ML approaches, the structure around a particle is encoded into a set of structural order parameters that capture e.g. the local density and symmetry of the distribution of neighbors around that particle. However, some sophisticated ML approaches can work from much more restricted data. For example, in Ref. 5, it was shown that graph neural networks are capable of predicting the dynamic propensity of a glassy Lennard-Jones mixture based solely on encoding the structure into a graph of nearest neighbors and the pair distances between them. Such data would not be sufficient for e.g. a simple linear regression approach, but using a GNN it was enough to drastically outperform both SVMs and CNNs trained with more sophisticated input. More recently, however, it was shown that even a simple linear regression approach can rival the predictive power of GNNs if supplied with sufficiently intelligent input data [START_REF] Boattini | Averaging local structure to predict the dynamic propensity in supercooled liquids[END_REF] . In particular, Boattini et al. proposed a method to iteratively construct generations of structural order parameters that successively take into account locally averaged order in expanding shells. These descriptors, in combination with linear regression preformed essentially as well as a GNN which was fed just the particle coordinates.

This raises the intriguing question of whether more sophisticated ML techniques supplied with more intelligently cho-sen structural parameters can result in even better predictions. Here, we use three different ML algorithms -linear regression, neural networks, and GNNs -to predict the dynamic propensity of a glassy binary hard-sphere mixture, based on the hierarchical set of order parameters from Ref. 8, and compare and contrast the results. As we show, out of the three methods, linear regression provides the best compromise between accuracy and efficiency when combined with these advanced structural descriptors.

I. MODEL AND DESCRIPTORS

A. Model

The glassy system that we use here to compare the three different ML methods is a binary hard-sphere mixture at packing fraction η = 0.58. It consists of hard spheres of two sizes, with a size ratio of σ S /σ L = 0.85, where σ L(S) is the diameter of a large (small) particle. The composition x L = N L /(N L + N S ) = 0.3, where N L(S) is the number of large (small) spheres. Note that this is the same glassy mixture as was studied in Refs. 6, 8, and 11.

We simulate the evolution of our system using event-driven molecular dynamics (EDMD) [START_REF] Rapaport | The Event-Driven Approach to N-Body Simulation[END_REF] . The simulations are performed in the microcanonical ensemble (constant number of large and small particles N L and N S , volume V , and kinetic energy U). The time unit of our simulation is defined as

τ = mσ 2 L /k B T
where k B is Boltzmann's constant and m is the particle mass. Note that we set the masses of all particles to be equal. All simulated systems contained 2000 particles in total (600 large, 1400 small).

To generate the initial configurations, we use a separate EDMD simulation in which the particles grow over time until the desired packing fraction η = 0.58 is reached. After this, the system is equilibrated for at least 10 5 τ. Note that from previous work [START_REF] Boattini | Autonomously revealing hidden local structures in supercooled liquids[END_REF] , we know that the relaxation time of this system is on the order of τ α = 10 4 τ.

B. Dynamic propensity

To characterize the dynamical heterogeneity in our glassy system, we use the dynamic propensity [START_REF] Widmer-Cooper | How reproducible are dynamic heterogeneities in a supercooled liquid?[END_REF][START_REF] Widmer-Cooper | On the study of collective dynamics in supercooled liquids through the statistics of the isoconfigurational ensemble[END_REF] . This quantity is closely related to the mean squared displacement, and captures how far individual particles on average move over time. To measure the dynamic propensity, the evolution of a glassy system is measured multiple times, each time starting from the same initial configuration, while assigning each particle a random velocity drawn from a Maxwell-Boltzmann distribution at the desired temperature. This ensemble is called the isoconfigurational ensemble [START_REF] Widmer-Cooper | How reproducible are dynamic heterogeneities in a supercooled liquid?[END_REF] . To obtain the propensity ∆r i (t) of particle i, we average the absolute distance it traveled over the time interval t over all trajectories

∆r i (t) = |r i (t) -r i (0)| iso , (1)
where • iso indicates the average taken over all trajectories in the isoconfigurational ensemble.

To measure the propensity, we average over simulations starting from 100 different initial snapshots, and for each initial snapshot we average over 50 trajectories with different initial velocities. The dynamic propensity is measured at a logarithmically spaced set of time intervals t, between t/τ = 10 -2 and t/τ = 10 5 .

To illustrate the behavior of the dynamic propensity, we plot in Fig. 1 the globally averaged dynamic propensity as a function of time.

C. Structural descriptors

To describe the local environments of particles, we use the structural order parameters used in Ref. 8, which consist of a combination of radial densities and angular functions measured in different shells around a particle.

For the radial functions, we use essentially the same descriptors as used in Refs. 5, 8, and 9. These descriptors measure a weighted particle density inside a spherical shell with a thickness of approximately 2δ at distance r with respect to a reference particle i. The functions are defined as

G (0) i (r, δ , s) = ∑ j =i,s j =s e - (r i j -r) 2 2δ 2 .
(

) 2
Here i is the reference particle, s is the particle type and r i j is the distance between particles i and j. The summation is carried out over all particles with particle type s j = s, which means that the radial density that is measured is type-specific.

The angular descriptors that we use are based on bond order parameters [START_REF] Lechner | Accurate determination of crystal structures based on averaged local bond order parameters[END_REF][START_REF] Steinhardt | Bond-orientational order in liquids and glasses[END_REF] . These bond order parameters expand the local environment in terms of spherical harmonics. To obtain the angular descriptors for a particle i, we first calculate the complex coefficients

q (0) i (l, m, r, δ ) = 1 Z ∑ i = j e - (r i j -r) 2 2δ 2 Y m l (r i j ). (3)
Here Y m l (r i j ) is the spherical harmonic of order l, with m an integer that runs from -l to l, and Z is a normalization factor given by

Z = ∑ i = j e - (r i j -r) 2 2δ 2 . ( 4
)
Note that although the summation runs over all particles, again the exponent makes sure that mainly particles within a spherical shell at distance r and thickness 2δ will contribute to q (0) i (l, m, r, δ ). Finally we sum over m to obtain the rotationally invariant angular descriptors

q (0) i (l, r, δ ) = 4π 2l + 1 m=l ∑ m=-l |q (0) i (l, m, r, δ )| 2 .
(

) 5
Due to the symmetries of the spherical harmonics, q (0) (l, r, δ ) for a certain l is expected to detect l-fold symmetry in the environment at the chosen distance r. Boattini et al. [START_REF] Boattini | Averaging local structure to predict the dynamic propensity in supercooled liquids[END_REF] showed that the propensity prediction of a particle improves significantly, when the prediction is based not only on the structural parameters associated with the particle itself, but also on averaged structural information of neighbouring particles. Inspired by the architecture of graph neural networks, this was done by recursively constructing higherorder averaged structural parameters, which are are defined as

x (n) i = ∑ j:r i j <r c x (n-1) j e -r i j /r c ∑ j:r i j <r c e -r i j /r c , (6)
where x i can be any of the radial or angular order parameters of particle i. Additionally, x

(n) i represents the n th generation of parameter x i , and the sum runs over all neighboring particles within a cutoff distance r c , including i itself. The cutoff value r c is chosen to be r c /σ L = 2.1, which approximately corresponds to the second minimum of the radial distribution function [START_REF] Boattini | Averaging local structure to predict the dynamic propensity in supercooled liquids[END_REF] . However, as already shown in Ref. 8, the exact value does not have a significant influence on our results.

In total, we consider 354 0th-generation structural descriptors: 162 radial descriptors and 192 angular descriptors. For the radial descriptors we use 46 equally spaced spherical shells in the interval r/σ L = [0.86, 2.0], 20 equally spaced spherical shells in the interval r/σ L = (2.0, 3.0] and 15 equally spaced spherical shells in the interval r/σ L = (3.0, 4.5]. For the angular descriptors we consider l = 1 to 12 in 16 equally spaced spherical shells in the interval [1, 2.5]. The full environment of each particle is then described with up to 3 generations of these 354 parameters each, leading to a total of 1062 parameters.

Before using the structural parameters as an input for the machine learning algorithms, they are standardized by evaluating

x st i = x i - x σ x , (7)
where x i is the vector containing all parameters associated with particle i, x st i is the standardized parameter vector, and where x and σ x are respectively the mean and standard deviation of the parameter vector considering all particles of the same species as i. The standardization ensures that all descriptors have zero mean and unit variance, which can be helpful when using regularization in machine learning techniques. This will be discussed in more detail below.

II. MACHINE LEARNING METHODS

In this paper, we compare three different machine learning approaches for predicting the dynamic propensity based on the structural parameters introduced above. In particular, we compare linear regression (LR), neural networks (NN), and graph neural networks (GNN). Unless otherwise specified, we train separate models for large and small particles, and separate models for each time interval at which we are trying to predict the dynamic propensity.

Note that each of these approaches has a number of hyperparameters that tune the model fitted by the method to the supplied training data. For example, this can be the number of layers inside the neural network, parameters controlling regularization techniques that reduce overfitting, or the learning rate of the optimization algorithm for NNs and GNNs.

A. Linear regression

Linear regression is the simplest of the three methods, and simply finds the best linear combination of all structural de-FIG. 2. A Neural network consisting of an input layer, hidden layers and an output layer. The orange and blue circles represent input and output parameters respectively. The red and green arrows represent the weights that connect each layer in the neural network.

scriptors to predict the dynamic propensity. Here, we make use of L2-regularization (also known as Ridge regression) to reduce overfitting [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF] . This approach penalizes large weights in the linear fit. Note that this is the reason we standardized our structural parameters: since the different parameters have the same mean and variance, the effect of the regularization on each parameter is the same. For linear regression the only hyperparameter that can be tuned is α, which sets the strength of the large-weight penalty in Ridge regression.

B. Neural network

Neural networks are loosely based on the biological neural networks that make up our brains. A neural network consists of multiple layers of connected nodes, see Fig. 2, which mimic the neurons and synapses in the brain. The first and last layer are respectively called the input and output layer, which in our case take the structural parameters of each particle as an input, and give the predicted propensity for a certain time as the output. All the layers that lie between the input and output layer are called hidden layers. In a fully connected feed-forward neural network, as we use here, each node in a specific layer is connected to all the nodes in the following layer.

Due to the connections between nodes, information can be passed through the neural network. Each connection between nodes is associated with a so called 'weight'. The values of the nodes in the hidden layer and the output layer are found by multiplying the values of the nodes in the previous layer with the associated weights and adding a bias. The result is then passed to a non-linear function, in our case a Rectified Linear Unit (ReLU) [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF] , to yield the value of the node. Hence, the value of a node a m in layer l is calculated as

a (l) m = f ∑ n w (l) mn a (l-1) n + b (l) m . ( 8
)
Here f is the ReLU function, w

mn is the weight associated with the connection between nodes n and m, a We train the neural networks using the Python package Pytorch [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF] , and in particular use an Adam optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] . This optimizer is an extension to the stochastic gradient descent procedure, and is used to find an efficient path to a locally optimal set of weights and biases via backpropagation [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF] .

For neural networks there are many more hyperparameters that can be tuned, including the number of layers, and the number of nodes in each layer, as well as parameters associated with the learning process, such as the learning rate and the batch size. As the neural network turns out to be very sensitive to overfitting, we also run the neural network with ridge regression and drop-out [START_REF] Srivastava | Dropout: A simple way to prevent neural networks from overfitting[END_REF] .

C. Graph neural networks

GNNs [START_REF] Battaglia | Relational inductive biases, deep learning, and graph networks[END_REF][START_REF] Battaglia | Interaction networks for learning about objects, relations and physics[END_REF][START_REF] Scarselli | The graph neural network model[END_REF] are a relatively new class of machine learning techniques that combine neural networks with a graph-like data structure. As there are many variations of GNNs, our description is necessarily somewhat specific to the GNN we use in this paper.

In contrast to LR and NNs, which consider each particle as a separate training example, the GNN takes in an entire configuration at once. As a result, the GNN does not make a propensity prediction on a single-particle basis, but instead simultaneously makes a prediction for all particles (of a chosen species) in a configuration of the system.

In a GNN, input data is first mapped to a data structure which consists of a graph that holds numerical data at its nodes and edges (see Figure 3). In our case, each node corresponds to a particle in the configuration for which we are trying to predict propensities. When two particles are closer to each other than a certain distance r c = 2.1σ L , the corresponding nodes in the graph are connected by an edge.

As an input, every node holds the structural parameters of the corresponding particle as well as its species, and each edge holds information about the distance between the two connected particles. Like a neural network, a GNN corresponds to a non-linear function of its input parameters, that gets calculated in multiple consecutive layers, called graph layers. Each layer itself is a non-linear function, which takes as an input a graph structure containing information on all nodes and edges, and outputs a new graph structure with the same graph topology but updated numerical data at the nodes and edges. For each node and edge, the update that takes place within this graph layer takes into account not only the information already at that node or edge, but also data from neighboring edges or nodes, as illustrated in Fig. 3. The internal functions that perform this update consist of fully connected feed-forward neural networks as described above. In addition to these graph layers, the full GNN incorporates an encoder layer, which maps the input node data to the data structure used in the graph layers, and a decoder, which makes a propensity prediction for each particle of a chosen species, based on the updated node data. Both of these layers are standard feed-forward neural networks. Note that analogous to the GNN in Ref. 5, our GNN additionally provides each graph layer with information about the graph data output by the encoder. In other words, a node update takes into account i) the node's current values, ii) the aggregated current values of all edges connected to the node, and iii) the node's values as they were just after the encoder layer.

The core benefit of the GNN is that the prediction for the propensity of a given particle can include information about the structure of multiple shells of neighboring particles. The distance over which information is included can be controlled by the number of included graph layers. As such, the GNN inherently includes the feature of recursively considering the average local structure of shells of neighbors.

III. RESULTS

In this section we first look at each method independently, and explore the influence of the different hyperparameter settings, and the inclusion of different generations of structural order parameters. In all cases, to compare the predictions to the measured propensity we use the Pearson correlation coefficient.

A. Linear regression

The only hyperparameter for linear regression is the regularization parameter α. As optimizing the fit for a single parameter is trivial, here we only present results that correspond to the best choice of α, optimized between 10 -5 and 10 4 for the large particles. In Figure 4 we show the linear regression performance for different generations of order parameters. In all cases, when we refer to a generation, we include all lower generations as well. Note that these results are consistent with Ref. 8. We clearly see that the predictions from the zeroth generation of descriptors are significantly worse than the ones including higher-generation data, at least for longer times. In particular, we see that the information of the higher-order generations only starts to influence the performance when the system enters the caging regime. This is what we expected: before entering the caging regime not enough time has passed for particles to be influenced by particles from further away, meaning that higher-order generations will not add relevant information about the expected trajectories. Although adding the second generation of order parameters still improves the predictions for the propensity, the effect is small in comparison to the improvement of adding the first generation. Adding even higher generations (not shown here) does not significantly improve the performance beyond this point [START_REF] Boattini | Averaging local structure to predict the dynamic propensity in supercooled liquids[END_REF] . I. Hyperparameters for different neural networks used to predict the propensity. Each row corresponds to a line in Fig. 5. The Hidden layers entry contains both information about the number of hidden layers in the network, as well as the number of nodes in each hidden layer (shown in brackets). Drop-out shows the fraction of nodes in the second hidden layer that is set to zero during drop-out and α represents the value of the parameter associated with L2-regularization.

B. Neural networks

For the neural networks, in addition to the number of generations and regularization parameter, there are many hyperparameters to optimize. Since trying out all possible combination of settings would not be feasible, we limit ourselves to a small number of different combinations, shown in Table I. The network is trained on large and small particles separately in around 500 epochs. In order to limit overfitting, after each epoch we evaluate the performance of the NN on the test data set, and eventually use the network that had the highest correlation for this test data. In total for each time interval and for each generation of structural parameters, we trained five networks. Their performance is shown in Fig. 5. Here, all networks are trained using three generations of structural parameters.

Although the overall behavior of the neural network accuracy over time is similar to that found for linear regression, we see significant variation in performance between the different hyperparameter choices. In particular, we see that only using a single NN layer (blue line) or no regularization (blue and yellow lines) lead to worse performance, especially at shorter times. However, once we include at least two NN layers and regularization, the results are relatively robust: the other three lines essentially coincide except for noise.

In order to examine the influence of different generations of order parameters on the performance, in Fig. 6 we show the performance of the NN for different generations, with the hyperparameters optimized for each time interval. Contrary to what we saw in the case of linear regression, providing the NN with more generations of order parameters does not always improve the performance, something that is especially clear for short times. This is likely the result of the higher number of weights and biases that the NN training needs to optimize when higher-order descriptors are included.

C. Graph neural networks

In Ref. 5, Bapst et al. demonstrated that hyperparameters did not play a large role in the accuracy of their GNNs for predicting dynamic propensity. Since the training of a GNN is considerably more expensive than a normal neural network, here we do not focus on optimizing the hyperparameters, but instead consider a few different variations of supplying data to the GNN. As a baseline, we consider a GNN with four graph layers, which considers the zero'th generation of structural parameters as the node inputs, and the absolute distance between particles as the edge inputs. The GNN is trained separately to predict the propensities of the large and small particles, but takes the information of all particles into account for both trainings. As a variation on this baseline, we also consider i) a GNN that predicts both the propensity of both the large and small species simultaneously, ii) a GNN that incorporates all three generations as node data, and iii) a GNN that uses the x, y, and z components of the vectors between neighboring particles as edge data. A complete summary of all relevant (hyper)parameters is shown in Table II.

In Fig. 7, we show the performance of the different vari-FIG. 7. Correlation between the dynamic propensity as predicted by graph neural networks and measured in simulations as a function of time, analyzed for different variations of the GNN model, which can be found in Table II.

ations of GNN. Clearly, none of the changes made to the GNN inputs have a significant effect on the overall performance. The most significant difference is associated with the networks that trained both small and large particles simultaneously -this adaptation performed worse for our system. In contrast to linear regression and neural networks, where taking into account additional generations clearly led to some improvement in the performance, evidently the inherent local averaging of the GNN is already sufficient to incorporate this type of structural information. Additionally, we observe that compared to the NN, the GNN is less sensitive to overfitting, resulting in relatively smooth lines. This is actually quite remarkable, since the number of weights and biases that need to be optimized in a GNN is significantly larger than in the case of a NN. Moreover, none of the training runs in Fig. 7 used any regularization methods. The fact that GNNs, compared to NNs, have less trouble converging and are less sensitive to overfitting, might be due to the fact that they are trained and evaluated on entire snapshots at once. Realistically, in a snapshot we know that there are strong correlations in mobility between neighboring particles. The fact that the GNN can take into account the mobility of neighbouring particles, will likely lead to smoother variation of the predicted propensity in space than the NN or LR, which results in fewer outliers.

D. Comparing the three methods

Finally, in Fig. 8, we compare the performance of all three methods. Again, we take for each time interval the bestperforming result from either LR, NN, or GNN. The results of the three methods are remarkably similar, suggesting that all three approaches are capable of extracting essentially the same information from the input data. Overall, the NN approach performs the worst, likely due to an inability to find the globally optimal solution to its training problem. Intuitively, the fit from linear regression could be reproduced essentially exactly by the NN, with a sufficiently good optimiza- show the number of hidden layers for the respectively the node and edge hidden layers together with the number of nodes in each layers (the number of output nodes for each of these networks is equal to 10), Node gen represents up to how many generations of structural order parameters we provide the GNN with, Edge par indicates whether the edge input is given by the absolute distance between particles (r), or the vector distance (x, y, z). Finally Tog or Sep represents whether we train the network together (Tog) for small and large particles, or train two separate networks (Sep).

tion. Hence, at any point where the NN performs less well than the LR solution, there is at least some failure to fully optimize the network. As we saw earlier, GNNs are less sensitive to their hyperparameters than NNs, and converge more easily. Moreover, during intermediate times in the caging regime and the beginning of the diffusive regime GNNs slightly outperform LR. This implies that the averaging that takes place in a GNN provides the network with slightly different information than the averaged parameters of the first and second generation. However, the improvement is extremely limited, and comes at the cost of a significantly more computationally costly training process. To illustrate this: for one choice of hyperparameters and the full set of time intervals, a typical training process takes approximately 3 minutes for LR, 24 hours for NNs, and 6.5 hours for GNNs. Note that these times are achieved by a laptop CPU for the LR, while the NN and GNN trainings made use of an Nvidia GeForce RTX 2080Ti GPU. We conclude that, given the discussion above, linear regression is the preferred method; it is fast, robust, and provides accurate predictions.

IV. CONCLUSION

In this paper, we compared three different ML methods for predicting the dynamic propensity of a glassy system at different times, namely linear regression, neural networks, and graph neural networks. We find surprisingly little difference in their performance over the full range of time intervals considered. The intuitive conclusion one can draw from the similar results of the three methods is that our prediction, at this point, is limited not by our fitting approach, but rather by the infor-mation contained in the set of structural order parameters. On the bright side, this means that given this set of descriptors, a simple and efficient ML method like linear regression is sufficient for essentially optimal predictions. On the other hand, it suggests that more advanced ML techniques are not likely to provide a solution to the question of how to further improve the prediction of dynamics in these systems. The main question that remains is: what information are we currently missing in order to improve our ability to predict the dynamic propensity, in particular in the caging regime, where the correlation between prediction and reality is currently minimal? While answering this question will require further research, our results here suggest that linear regression is likely a sufficient method for evaluating the predictive capabilities of new sets of structural order parameters.

FIG. 1 .

 1 FIG. 1. Globally averaged dynamical propensity ∆r i (t) as a function of time. The insets at the bottom illustrate three different dynamical regimes: the ballistic regime, in which particles are not yet interacting with their neighbors, the caging regime, in which particles are trapped by surrounding particles, and the diffusive regime, where particles have escaped their cages.

  (l-1) n is the information of node n in layer l -1 and b (l) m is the bias associated

FIG. 3 .

 3 FIG. 3. (left)Mapping between the initial configuration data and the graph structure. From the configuration, we determine nearest neighbors in order to set up a graph structure, and then calculate structural parameters to fill in the data at each node and edge. (right) Node and edge update in a graph layer. To update the node data, the added and averaged parameters of the connected edges are given to the node neural network, together with the information of the node itself. For the edge update, the information of both neighbouring nodes, together with the edge information itself is given to the edge neural network. The output of the node/edge network is an updated value of the node/edge.

FIG. 4 .

 4 FIG. 4. Pearson's correlation coefficient ρ between the dynamic propensity as predicted by linear regression and measured in simulations as a function of time, analyzed for three different generations of order parameters. Results are shown for the large particles.

FIG. 5 .

 5 FIG. 5. Correlation between the dynamic propensity as predicted by neural networks and measured in simulations as a function of time, analyzed for different hyperparameters as specified in Table I. Results are shown for the large particles.

FIG. 6 .

 6 FIG. 6. Correlation between the dynamic propensity as predicted by neural networks and measured in simulations as a function of time, analyzed for three different generations of order parameters. Results are shown for the large particles.

FIG. 8 .

 8 FIG. 8. Comparison of the accuracy of the dynamic propensity prediction over time as obtained with LR, NN and GNN. Results are shown for both the (a) large and (b) small particles. For each point in time we use for each method the set of hyperparameters that resulted in the best performance.

TABLE II .

 II Hyperparameters for different graph neural networks used to predict the propensity. Each row corresponds to a line in Fig.7. In the table the following abbreviations are used: GL is the number of graph layers, Node enc H.L. and Edge enc H.L. represent respectively the number of hidden layers in the node and edge encoder (for both, the number of input parameters is equal to the number of parameters associated with a node or an edge, while the number of output parameters is equal to 10. Node H.L and Edge H.L.

	Nr.	Batch size	Learning rate	GL	Node enc H.L.	Edge enc H.L.	Node H.L	Edge H.L.	Node gen	Edge par	Tog or Sep
	1	5	10 -4	4 2 (50, 50)	1 (5)	2 (30, 16) 2 (16, 16)	2 (16, 16) 2 (16, 16)	1	r	Sep
	2	5	10 -4	4 2 (50, 50)	1 (5)	2 (30, 16) 2 (16, 16)	2 (16, 16) 2 (16, 16)	3	r	Sep
	3	5	10 -4	4 2 (50, 50)	1 (5)	2 (30, 16) 2 (16, 16)	2 (16, 16) 2 (16, 16)	1	x, y, z Sep
	4	5	10 -3	4 2 (50, 50)	1 (5)	2 (30, 16) 2 (16, 16)	2 (16, 16) 2 (16, 16)	1	x, y, z Sep
	5	5	10 -4	4 2 (50, 50)	1 (5)	2 (30, 16) 2 (16, 16)	2 (16, 16) 2 (16, 16)	1	r	Tog
	6	3	10 -4	3 2 (50, 50)	1 (5)	2 (30, 16) 2 (16, 16)	2 (16, 16) 2 (16, 16)	1	r	Tog

ACKNOWLEDGEMENTS

The authors would like to thank Marjolein de Jager for many discussions. L.F. and E.B. acknowledge funding from The Netherlands Organisation for Scientific Research (NWO) (Grant No. 16DDS004), and L.F. acknowledges funding from NWO for a Vidi grant (Grant No. VI.VIDI.192.102).