Local convergence of multi-agent systems towards triangular patterns

Geometric pattern formation is an important emergent behavior in many applications involving large-scale multi-agent systems, such as sensor networks deployment and collective transportation. Attraction/repulsion virtual forces are the most common control approach to achieve such behavior in a distributed and scalable manner. Nevertheless, for most existing solutions only numerical and/or experimental evidence of their convergence is available. Here, we revisit the problem of achieving pattern formation giving sufficient conditions to prove analytically that under the influence of appropriate virtual forces, a large-scale multi-agent swarming system locally converges towards a stable and robust triangular lattice configuration. Specifically, the proof is carried out using LaSalle's invariance principle and geometry-based arguments. Our theoretical results are complemented by exhaustive numerical simulations confirming their effectiveness and estimating the region of asymptotic stability of the triangular configuration.


I. INTRODUCTION
Many natural and artificial systems consist of multiple interacting agents; their behavior being determined by both the individual agent dynamics and their interaction.In some applications the number of agents can be extremely large (large-scale multi-agent systems) and the role played by their interconnections becomes predominant over their individual dynamics [1].Examples include cell populations [2], swarming multi-robot systems [3], social networks [4] among many others.Some of the most relevant emerging behavior exhibited by these systems involve their spatial organization, coordination, and cooperation [5].A notable case is geometric pattern formation [6] where the agents are required to self-organize into some desired pattern, such as, for example, triangular lattices consisting of repeating adjacent triangles.Applications of pattern formation include sensor networks deployment [7], collective transportation and construction [8], [9], and exploration and mapping [10].
Most of the existing distributed control algorithms for geometric pattern formation rely on the use of virtual forces (or virtual potentials), [7], [11]- [18].Within this framework, agents move under the effect of forces generated by the presence of their neighboring agents and the environment, causing attraction, repulsion, alignment, etc.
Interestingly, most strategies are validated only numerically or experimentally [7], [11]- [13].Among the exceptions, in [19], a geometric control approach based on trigonometric functions is proposed to build triangular lattices, and its global convergence is proved.The extension to 3D spaces is validated analytically in [20].Moreover, harmonic approximation [21] provides necessary conditions for the local stability of a lattice.These conditions are used in [14] to numerically design a virtual force that locally stabilizes an hexagonal lattice.A general analysis of the effects of attraction/repulsion virtual forces is carried out in [22], where the authors prove that the agents converge inside a bounded region, even though the specific equilibrium configuration is not characterized.We wish to remark here that formation control [15]- [17] differs from geometric pattern formation because of a typically smaller number of agents (order of tens) with, possibly, unique identifiers, numerous roles for the agents and often some coordinated motion of the agents.Similarly, when solving flocking control problems, the emergence of coordinated motion is the crucial concern [18], [23], [24].
In this paper, we revisit the problem of geometric pattern formation using attraction/repulsion virtual forces with the aim of bridging a gap in the existing literature and deriving a general proof of convergence when considering the formation of triangular lattice configurations.When compared to previous work, our stability results (i) can be applied to most control laws based on virtual forces (or potentials), rather than only holding for specific algorithms, e.g.[19], (ii) are sufficient rather than being only necessary [21], (iii) characterize the asymptotic configuration of the agents, rather than just proving its boundedness [22], and (iv) guarantee the emergence of triangular lattices rather than less regular ones, e.g.α-lattices studied in [18].

II. MATHEMATICAL PRELIMINARIES
Given a vector v ∈ R d , we denote by [v] i its i-th element, by v its Euclidean norm, and by v := v v its direction.0 denotes a column vector of appropriate dimension with all elements equal to 0. Given a matrix A, [A] i j is its (i, j)-th element.
Given a continuous-time, autonomous dynamical system with state vector x(t) ∈ R d , and x 0 ∈ R d , we term as φ (t, x 0 ) its trajectory starting from x(0) = x 0 .Definition 1 (Equilibrium set): Definition 2 (Local asymptotic stability [25, Definition 1.8]):An equilibrium set Ξ for system (1) is locally asymptotically stable if ∀ε > 0, ∃δ > 0 such that if min y∈Ξ x 0 − y < δ , then 1) min y∈Ξ φ (t, x 0 ) − y < ε, ∀t > 0, and 2) lim t→+∞ φ (t, x 0 ) ∈ Ξ. Definition 3 (Incidence matrix): Given a digraph with n vertices and m edges, its incidence matrix B ∈ R n×m has elements defined as Moreover, the length of an edge, say (i, j) ∈ E, is p i − p j .Definition 5 (Congruent frameworks [26, p. 3]): Given a graph G = (V, E) and two frameworks (G, p) and (G, q), these are congruent if p i − p j = q i − q j ∀i, j ∈ V. Definition 6 (Rigidity matrix [26, p. 5]): Given a ddimensional framework with n ≥ 2 vertices and m edges, its rigidity matrix M ∈ R m×dn has elements defined as i and ends in vertex j, [p i − p j ] k , if edge e starts from vertex j and ends in vertex i, 0, otherwise.
(2) with k = 1, . . ., d. Definition 7 (Infinitesimal rigidity [16, p. 122]): A framework with rigidity matrix M is infinitesimally rigid if, for any infinitesimal motion, say u, 1 of its vertices, such that the length of the edges is preserved, it holds that Mu = 0.
To give a geometrical intuition of the concept of infinitesimal rigidity, we note that an infinitesimally rigid framework is also rigid [16, p. 122], according to the definition below. 2  Definition 8 (Rigidity [26, p. 3]): A framework is rigid if every continuous motion of the vertices, that preserves the length of the edges, also preserves the distances between all pairs of vertices.
As a consequence, in a rigid framework, any continuous motion that does not preserve the distance between any two pairs of vertices also does not preserve the length of at least one edge.Theorem 1 ( [16, p. 122]): A 2-dimensional framework with n ≥ 2 vertices and rigidity matrix M is infinitesimally rigid if and only if rank(M) = 2n − 3. Definition 9 (Swarm): A (planar) swarm S := {1, 2, . . ., n} is a set of n ∈ N >0 identical agents that can move on the 1 u can be interpreted as either a velocity or a small displacement. 2In rare cases, a rigid graph is not infinitesimally rigid; e.g.[26, p. 7].plane.For each agent i ∈ S, x i (t) ∈ R 2 denotes its position in the plane at time t ∈ R ≥0 .
Moreover, we call x(t as its center, and denote by r i j (t) := x i (t) − x j (t) ∈ R 2 the relative position of agent i with respect to agent j.Definition 10 (Adjacency set): Given a swarm S, the adjacency set of agent i at time t is A i (t) := { j ∈ S \ {i} : ,where R a ∈ R >0 is the maximum link length.Definition 11 (Links): A link is a pair (i, j) ∈ S × S such that j ∈ A i (t); r i j (t) is its length.The set of all links existing in a certain configuration x is denoted by E(x).
Notice that (i, j) ∈ E(x) ⇐⇒ ( j, i) ∈ E(x).Here, we assume that so that, when the swarm is in a triangular configuration, the adjacency set (Definition 10) of any agent includes only the agents in its immediate surroundings, and all the links (Definition 11) have length R (see Fig. 1).We denote by T ⊂ R 2n the set of all triangular lattice configurations; it is immediate to verify that T is unbounded and disconnected.Definition 14 (Congruent configurations): Given a configuration x , we define the set of its congruent configurations Γ(x ) as the set of configurations with congruent associated frameworks (see Definition 5), that is These configurations are obtained by translations and rotations of the framework F(x ); thus, it is immediate to verify that Γ(x ) is connected and unbounded for any x (see Fig. 2a).Also, note that x * ∈ T ⇐⇒ Γ(x * ) ⊂ T , and In the following, we omit the dependence on time when clear from the context.

III. PROBLEM STATEMENT
Consider a swarm S of n agents, with agents' dynamics described by where u i (t) ∈ R 2 is a control law to be designed.Let R s ∈ R >0 be a sensing radius and define the interaction set of agent i at time t as For the term u i (t) in ( 5), we consider the distributed virtual forces control law, given by where f : R ≥0 → R is the interaction function.Note that in general there is no specific relation between I i and A i (see Definition 10); however, we reasonably assume that R s ≥ R a , so that The following result slightly extends the one reported in [22, Lemma 1].Lemma 1: The position of the center of the swarm (see Definition 9), say x c , under the control law Proof.Exploiting ( 5) and (7), the dynamics of the center of the swarm is given by ẋc Since in a swarm the existence of any link (i, j) implies the existence of link ( j, i) (see Definition 10), in (9), for any term f ( r i j ) ri j there exists a term f ( r ji ) r ji = − f ( r i j ) ri j (because r i j = r ji and ri j = −r ji ).Therefore, the sum of the two is zero, yielding the thesis.

IV. CONVERGENCE TO A TRIANGULAR CONFIGURATION
We can now state the main result of this work, showing that, given an interaction function f (in (7)) that generates short range repulsion and long range attraction, the set of triangular configurations of the swarm is a locally asymptotically stable equilibrium set (see Definitions 1 and 2).
An exemplary interaction function fulfilling the assumption above is portrayed in Fig. 3.
Without loss of generality, we further assume that, under Assumption 1, in a sufficiently small neighborhood of a triangular configuration, all other equilibria are also triangular (supporting evidence showing that this assumption is not restrictive is reported in the Appendix).Theorem 2: Let Assumption 1 hold.Then, for any triangular configuration x * , Γ(x * ) is a locally asymptotically stable equilibrium set.Consequently, T is also a locally asymptotically stable equilibrium set.
Proof.Let us consider any triangular configuration x * ∈ T , with center x * c := 1 n ∑ n i=1 x * i and relative positions r * i j , and the set Γ(x * ) of its congruent configurations.Recalling Definition 13.(B) and (a1), we have that x * is an equilibrium point of ( 5)- (7); thus, Γ(x * ) and T are equilibrium sets.Next, we will prove local asymptotic stability of Γ(x * ) ⊂ T , which implies local asymptotic stability of T through (4).
Step 1 (Lyapunov function): Given a configuration x ∈ R 2n with center x c and inducing the links in E(x) according to Definition 11, let m := |E(x)| and order the links in E(x) arbitrarily, so that r 1 , . . ., r m refer to the relative positions r i j for (i, j) ∈ E(x).Recalling (a3), we can define the potential function P : [0, R a ] → R given by P(z) = − z R f (y) dy (see Fig. 3).Note that P(R) = 0, dP dz (z) = − f (z), and, from (a2), Then, let us consider the candidate Lyapunov function By (10), it holds that V (x) ≥ 0 ∀x ∈ R 2n , and V = 0 if and only if both x c = x * c and Definition 13.(B) holds.
Step 2 (Properties of V ): V (x) is discontinuous over R 2n (because E(x) changes when links (dis-)appear).However, V (x) is continuous and differentiable in any subset of R 2n where the set E(x) of links is constant.To find such a set, we seek conditions on x such that E(x) = E(x * ) (see Definitions 10 and 11), i.e., (12a) means that all links in E(x * ) are preserved in E(x), while (12b) means that no new links are created in E(x) with respect to E(x * ).With simple algebraic manipulations it is possible to show that (12a) and (12b) hold if x ∈ B, where and β < min i, j∈S R a − r * i j .Note that B can be interpreted as a "neighborhood" of Γ(x * ) with "width" β (see Fig. 2b).Hence, E(x) = E(x * ) in B, and thus V is continuously differentiable in B.
Step 3 (Analysis of V ): At this point, we can restrict our analysis to the set B to study the attractivity of Γ(x * ).Let us start by studying the dynamics of the agents.From ( 5)-( 7), we have Hypothesis (a4) and (8) imply that in (14) we have Then, exploiting (15) and the incidence matrix B (Definition 3) of the swarm graph, ( 14) can be rewritten as Moreover we can write the dynamics of the relative positions along a link k as ṙk = ∑ n i=1 [B] ik ẋi .Thus, exploiting Lemma 1 and ( 16), we get We can hence conclude that V (x) = 0 if and only if ẋ = 0, i.e., in correspondence of equilibrium configurations.Now, choosing β small enough, we can exclude the presence of equilibrium configurations not belonging to Γ(x * ), and therefore Step 4 (Applying LaSalle's invariance principle): To complete the proof, we define a forward invariant neighborhood of x * and then apply LaSalle's invariance principle.Given some ω ∈ R >0 , let Ω be the largest connected set containing x * such that V (x) ≤ ω ∀x ∈ Ω (see Fig. 2b).In particular, we select ω small enough that Ω ⊆ B. 3 Since V (x) ≤ ω and V (x) ≤ 0 for all x ∈ Ω, then Ω is forward invariant.Moreover, Ω is closed, because V is continuous in Ω, and Ω is the inverse image of the closed set [0, ω].Ω is also bounded because (i) translations too far from x * cause V to increase beyond ω (see (11)), and (ii) Ω ⊆ B implies that the deformations of the framework are bounded (see (13)).
Since Ω is closed and bounded, it is also compact.
As Ω is compact and forward invariant, we can apply LaSalle's invariance principle [27,Theorem 4.4], and noting that, in Ω, V (x) = 0 if and only if x ∈ Γ(x * ) (see ( 17)), we get that all the trajectories starting in Ω converge to Γ(x * ) ∩ Ω.This and the forward invariance of Ω imply that Γ(x * ) is locally asymptotically stable, and so is T because of (4).

V. NUMERICAL VALIDATION
In this section, we validate numerically the result presented in Section IV and estimate the basin of attraction of T .

A. Simulation setup
We set the desired link length to R = 1, the maximum link length to R a = (1 + √ 3)/2 ≈ 1.37, the sensing radius to R s = 3, and the number of agents to n = 100.
The interaction function f is chosen as the Physicsinspired Lennard-Jones function [5], [11], given by where we select a = b = 0.5 and c = 12; see Fig. 4. In (18), f is saturated to 1 to avoid divergence of f for z → 0.
Concerning Assumption 1, the interaction function f satisfies (a1), (a2), and (a3).Also, as shown in Fig. 4, it quickly tends to zero so that we can assume it practically satisfies (a4).The choice of not setting f (z) exactly equal to zero for z ≥ R a is intentional as it allows to account for long range attraction between the agents, which is frequently required in swarm robotics applications [22].
To assess if the swarm is in a triangular configuration, we check the conditions in Definition 13.To evaluate whether a configuration is infinitesimally rigid, we use Theorem 1.Moreover, we define the error e(t) := max k∈E(t) | r k (t) − R|, which is zero when the configuration is triangular.Also, as long as e(t) is lower than R a − R, links in the configuration of interest are neither created nor destroyed.
For each simulation, the initial positions of the agents are obtained by picking a random triangular configuration and then applying, to each agent, a random displacement drawn from a uniform distribution over a disk of radius δ ∈ R ≥0 .All simulation are run in MATLAB 4 and last 20 s; the agents' dynamics ( 5)-( 7) are integrated using the forward Euler method with a fixed time step equal to 0.01 s.Remark 1: Theorem 2, together with Definitions 9 and 13 allow a straightforward extension of the analysis to the threedimensional case (d = 3).The only cumbersome step is to assess the infinitesimal rigidity of the 3D framework of interest as Theorem 1 can no longer be applied.

B. Numerical results
To validate Theorem 2 and estimate the basin of attraction of the set of triangular configurations, we performed extensive simulations for various values of δ , and observed the steady state configurations.The results are reported in Fig. 5. Namely, we see that for δ ≤ δ thres := 0.25 all simulations converge to a triangular configuration, with a rigid framework and a negligible value of e.Then, as δ increases beyond δ thres , the average number of simulations converging to triangular configurations decreases, until for δ > 0.45 no simulation converges to a triangular configuration.Notice that e(0) ≤ 2δ , therefore δ = 0.25 corresponds to a perturbation of up to 50% of the initial length of the links, providing an estimation of the basin of attraction (region of asymptotic stability) of T .Moreover, we analysed the time evolution of e(t) in the case δ = 0.2.The results of 10 simulations are shown in Fig. 6.We find that the rigidity is preserved during all simulations, and at steady state e reaches zero, meaning that the swarm, when locally perturbed, quickly converges back to a triangular configuration, as expected from Theorem 2.

VI. CONCLUSIONS
We proved analytically local asymptotic stability of triangular lattice configurations for planar swarms under the action of a distributed control action based on virtual attraction/repulsion forces.The theoretical derivations were supported by exhaustive numerical simulations validating the theoretical results and providing an estimate of the basin of attraction.The mild hypotheses required on the interaction function that were used to prove convergence allow for wide applicability of the theoretical results.
Future work will focus on the formalization of the threedimensional case and the extension of the results to other geometric lattices, such as squares and hexagons. 4Code available at https://github.com/diBernardoGroup/SwarmSimPublic.

APPENDIX
To confirm the effectiveness of our theoretical results, we provide below further semi-analytical evidence that the set of triangular configurations T , is locally asymptotically stable, which also excludes the presence of other equilibria in an arbitrarily small neighborhood of it.To do so, we linearize system (5)-( 7) around a triangular configuration, say x * , obtaining ẋ ≈ J(x * ) (x − x * ), with J(x * ) ∈ R 2n×2n derived as follows.
Jacobian of (5)- (7): System ( 5)-( 7) can be recast as ẋ where F, G, H ∈ R m×m are diagonal matrices; [F] ii := f ( r i ), [G] ii := r i , and H := FG −1 .The Jacobian of ( 19) is where Numerical analysis: We set R = 1 and generated 760 random triangular configurations (10 per each number of agents n between 25 and 100).For each of these configurations, assuming f (in (7)) is in the form (18), we computed J using ( 20)- (21) and found that in all cases J has 3 zero eigenvalues with eigenvectors {w 0 i } i , and 2n − 3 negative eigenvalues with eigenvectors {w ± j } j .Moreover, Mw 0 i = 0 and Mw ± j = 0; thus, from Definition 7, the span of {w 0 i } corresponds to roto-translations and is a hyperplane locally tangent to Γ(x * ) (see Definition 14), while {w ± j } correspond to other motions.Therefore, the center manifold theorem [25, Theorem 5.1] yields that Γ(x * ) is a center manifold of system ( 5)- (7).Moreover, as expected from Theorem 2, the reduction principle [25, Theorem 5.2] confirms that the dynamics locally converge onto the equilibrium set Γ(x * ), and excludes the presence of other equilibria in an arbitrarily small neighborhood of it.

Fig. 1 .
Fig. 1.Triangular configurations.(a) Schematic representation of a triangular lattice; red agents belong to the adjacency set of the black agent.(b) Example of a triangular configuration with n = 100 agents.

Fig. 3 .
Fig. 3. Example of an interaction function f (top panel) and its corresponding potential P (bottom panel).

Fig. 4 .
Fig.4.Plot of the interaction function defined by(18).The zero of the function is highlighted by a red dot.

Fig. 5 .
Fig. 5. Simulations for different values of δ .(a): Terminal values of e and ρ. ρ is the fraction of trials converging to an infinitesimally rigid configuration.For e, the solid line is the mean; the shaded area is the minimum and maximum.20 simulations with random initial conditions are performed for each value of δ .(b), (c): Initial and final configurations of representative simulations for specific values of δ .

Fig. 6 .
Fig.6.Time evolution of the error e in 10 simulations with random initial conditions and δ = 0.2; the solid line is the mean, while the shaded area is the maximum and minimum.