Distributed Control for Spatial Self-Organization of Multi-Agent Swarms

TL;DR

This paper presents a distributed control approach for multi-agent swarms to self-organize spatially without position data, achieving desired density distributions in 1D and 2D domains using Laplacian-based algorithms.

Contribution

It introduces a novel pseudo-localization method and control laws enabling swarms to self-organize spatially without position information, applicable in both 1D and 2D.

Findings

Successful design of distributed control laws for density shaping.

Effective pseudo-localization algorithm for agents without position data.

Validation of methods in both 1D and 2D spatial domains.

Abstract

In this work, we design distributed control laws for spatial self-organization of multi-agent swarms in 1D and 2D spatial domains. The objective is to achieve a desired density distribution over a simply-connected spatial domain. Since individual agents in a swarm are not themselves of interest and we are concerned only with the macroscopic objective, we view the network of agents in the swarm as a discrete approximation of a continuous medium and design control laws to shape the density distribution of the continuous medium. The key feature of this work is that the agents in the swarm do not have access to position information. Each individual agent is capable of measuring the current local density of agents and can communicate with its spatial neighbors. The network of agents implement a Laplacian-based distributed algorithm, which we call pseudo-localization, to localize themselves in a new coordinate frame, and a distributed control law to converge to the desired spatial density distribution. We start by studying self-organization in one-dimension, which is then followed by the two-dimensional case.

Figures

Equations

$$d_H(M_1,M_2)=\max \{\underset{x\in M_1}{\sup }\underset{y\in M_2}{\inf }d(x,y),\underset{y\in M_2}{\sup }\underset{x\in M_1}{\inf }d(x,y)\}.$$

$${\int }_{\Omega }(\nabla \cdot \mathbf{F})d\mu ={\int }_{\partial \Omega }\mathbf{F}\cdot \mathbf{n}dS,$$

$${\int }_{\Omega }(\mathbf{F}\cdot \nabla U)d\mu ={\int }_{\partial \Omega }U(\mathbf{F}\cdot \mathbf{n})dS-{\int }_{\Omega }U(\nabla \cdot \mathbf{F})d\mu .$$

$$\frac{d}{dt}\left({\int }_{\Omega (t)}f(t,\mathbf{r})d\mu \right)={\int }_{\Omega (t)}{\partial }_t(f(t,\mathbf{r}))d\mu +{\int }_{\partial \Omega (t)}f(t,\mathbf{r})\mathbf{v}\cdot \mathbf{n}dS.$$

$$U=\frac{1}{2}{\int }_{\Omega }|f{|}^{2}d\mu ,$$

$$\dot{U}={\int }_{\Omega }f\cdot \left(\frac{df}{dt}\right)d\mu +\frac{1}{2}{\int }_{\Omega }|f{|}^2\nabla \cdot \mathbf{v}d\mu .$$

$$\frac{\partial U}{\partial t}$$

$$={\int }_{\Omega }f\cdot f_t+\frac{1}{2}{\int }_{\Omega }\nabla \cdot (|f{|}^2\mathbf{v})$$

$$={\int }_{\Omega }f\cdot f_t+{\int }_{\Omega }f\cdot (\mathbf{v}\cdot \nabla )f+\frac{1}{2}{\int }_{\Omega }|f{|}^2\nabla \cdot \mathbf{v}$$

$$={\int }_{\Omega }f\cdot (f_t+(\mathbf{v}\cdot \nabla )f)+\frac{1}{2}{\int }_{\Omega }|f{|}^2\nabla \cdot \mathbf{v}$$

$$={\int }_{\Omega }f\cdot \left(\frac{df}{dt}\right)+\frac{1}{2}{\int }_{\Omega }|f{|}^2\nabla \cdot \mathbf{v}.$$

$$\Vert u-u_{\Omega }{\Vert }_{L^p(\Omega )}\le C\Vert \nabla u{\Vert }_{L^p(\Omega )},$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{v})=0,\forall \mathbf{r}\in \mathring{M}(t),$$

$$E(\varphi )={\int }_M|\nabla \varphi {|}^{2}d{v}_g,$$

$$\Theta (x)={\int }_0^x\rho (\bar{x})d\bar{x},$$

$$\begin{array}{rl}& {\partial }_{t}X=\frac{1}{\rho }{\partial }_x\left(\frac{{\partial }_{x}X}{\rho }\right),\\ & X(t,0)=\alpha (t),X(t,L)=\beta (t),X(0,x)=X_0(x),\\ & \dot{\alpha }(t)=-\alpha (t),\dot{\beta }(t)=1-\beta (t),\end{array}$$

$$\begin{array}{rl}& {\partial }_{t}w=\frac{1}{\rho }{\partial }_x\left(\frac{{\partial }_{x}w}{\rho }\right),\\ & \frac{d}{dt}w(t,0)=-w(t,0),\frac{d}{dt}w(t,L)=-w(t,L),w(0,x)=X_0(x)-\Theta (x).\end{array}$$

$$V=\frac{1}{2}{\int }_M\rho |w{|}^2+\frac{1}{2}{\int }_M\frac{1}{\rho }|{\partial }_{x}w{|}^2.$$

$$\dot{V}={\int }_M\rho w({\partial }_{t}w)+{\int }_M\frac{1}{\rho }({\partial }_{x}w)({\partial }_t{\partial }_{x}w).$$

$$\dot{V}$$

$$\displaystyle=-\int_{M}\frac{1}{\rho}\left|\partial_{x}w\right|^{2}-\int_{M}\frac{1}{\rho}\left|\partial_{x}\left(\frac{\partial_{x}w}{\rho}\right)\right|^{2}+\frac{w+\partial_{t}w}{\rho}\partial_{x}w\bigg{|}_{L}-\frac{w+\partial_{t}w}{\rho}\partial_{x}w\bigg{|}_{0}.$$

$$\dot{V}=-{\int }_M\frac{1}{\rho }{\left|{\partial }_{x}w\right|}^2-{\int }_M\frac{1}{\rho }{\left|{\partial }_x\left(\frac{{\partial }_{x}w}{\rho }\right)\right|}^2.$$

$$\begin{array}{rl}& {\partial }_{t}X=\frac{1}{\rho }{\partial }_x\left(\frac{{\partial }_{x}X}{\rho }\right)-v{\partial }_{x}X,\\ & X(t,0)=0,X(t,L(t))=\beta (t),X(0,x)=X_0(x).\end{array}$$

$$\begin{array}{r}\frac{dX}{dt}={\partial }_{t}X+v{\partial }_{x}X=\frac{1}{\rho }{\partial }_x\left(\frac{{\partial }_{x}X}{\rho }\right)={\partial }_{\theta }\left({\partial }_{\theta }X\right)=\frac{{\partial }^{2}X}{\partial {\theta }^2}.\end{array}$$

$$\begin{array}{rl}& X_i(t+1)=\frac{1}{3}\left(X_{i-1}(t)+X_i(t)+X_{i+1}(t)\right),\\ & X_l(t)=0,X_r(t)=\beta (t),X_i(0)={X_0}_i.\end{array}$$

$$\begin{array}{rl}& {\partial }_t\rho =-{\partial }_x(\rho v),\\ & {\partial }_{t}X=\frac{1}{\rho }{\partial }_x\left(\frac{{\partial }_{x}X}{\rho }\right)-v{\partial }_{x}X,\\ & X(t,0)=0,X(t,L(t))=\beta (t),X(0,x)=X_0(x).\end{array}$$

$$\begin{array}{rl}& v(t,0)=0,{\partial }_{x}v=(\rho -p^{\ast }\circ X)-\frac{{\partial }_X{p}^{\ast }}{\rho (\rho +p^{\ast }\circ X)}{\partial }_x\left(\frac{{\partial }_{x}X}{\rho }\right),\end{array}$$

$$\displaystyle\begin{aligned} &X(t,0)=0,\hskip 14.45377pt\beta_{t}=k\left(2-\beta(t)-\frac{X_{x}}{\rho}\bigg{|}_{L(t)}\right).\end{aligned}$$

$$V=\frac{1}{2}{\int }_0^{L(t)}|\rho -p^{\ast }\circ X{|}^{2}dx+\frac{1}{2}{\int }_0^{L(t)}\rho |w{|}^{2}dx+\frac{1}{2}|w(L(t)){|}^2.$$

$$\dot{V}=$$

$$\displaystyle+\int_{0}^{L(t)}\rho w\partial_{t}w~{}dx+\frac{1}{2}\int_{0}^{L(t)}(\partial_{t}\rho)|w|^{2}~{}dx+\frac{1}{2}\rho|w|^{2}v\bigg{|}_{0}^{L(t)}+w(L)\frac{dw(L(t))}{dt}.$$

$$\dot{V}=$$

$$\displaystyle+\int_{0}^{L(t)}w\partial_{x}\left(\frac{\partial_{x}w}{\rho}\right)~{}dx-\int_{0}^{L(t)}\rho vw\partial_{x}w~{}dx-\frac{1}{2}\int_{0}^{L(t)}\partial_{x}(\rho v)|w|^{2}~{}dx+\frac{1}{2}\rho|w|^{2}v\bigg{|}_{0}^{L(t)}$$

$$+w(L)\frac{dw(L(t))}{dt}.$$

Full text

{AMS}

34B45, 35B40, 35B35, 58J32, 58J35, 58E20

1 Introduction

Self-organization in swarms refers broadly to the emergence of patterns of long-range order in large groups of dynamic agents which interact locally with each other. It is a pervasive phenomenon in nature, observed in biological [7] and other natural systems [36]. In the context of robotic systems, problems of deployment and formation control of groups of robots have been extensively studied [6, 27, 11, 32, 21]. More recently, research efforts have been undertaken to massively increase the scale of these robotic systems [30]. This transition does not merely involve an increase in the size of robotic networks, but it also introduces new theoretical challenges for their analysis and control design. In particular, large groups of agents have some essential characteristics that distinguish them from other smaller-scale counterparts. In a swarm, individual agents have no significance and only the macroscopic objectives are relevant. A swarm largely remains unaffected by the removal of a large, but discrete, number of agents. Moreover, it is difficult (and needlessly complicated) to specify the global configuration of the swarm using the states of individual agents; instead, employing macroscopic quantities such as the swarm spatial density distribution to specify its configuration is more appropriate. From an analysis and control-theoretic viewpoint, the dynamic modeling of swarms is less explored, which e.g. can be established by means of PDEs, for which control theoretic tools are less well developed in comparison to ODEs. These theoretical challenges motivate the investigation of self-organization in large-scale swarms.

In the literature, Markov-chain based methods have been widely used in addressing some of the key theoretical problems pertaining to swarm self-organization. By means of it, the swarm configuration is described through the partitioning the spatial domain in a finite number of larger size disjoint subregions, on which a probability distribution is defined. Then, the self-organization problem is reduced to the design of the transition matrix governing the evolution of this probability density function to ensure its convergence to a desired profile. A recent approach to density control using Markov chains is presented in [12], which includes additional conflict-avoidance constraints. In this setting every agent is able to determine the bin to which it belongs at every instant of time, which essentially means that individual agents have self-localization capabilities. Also, the dimensional transition matrix is synthesized in a central way at every instant of time by solving a convex optimization problem. In [3], the authors make use of inhomogeneous Markov chains to minimize the number of transitions to achieve a swarm formation. In this approach, the algorithm necessitates the estimation of the current swarm distribution, and computes the transition Markov matrices for each agent, at each instant of time. The fact that every agent needs to have an estimate of the global state (swarm distribution) at every time may not be desirable or feasible. The localization of each agent still remains to be a main assumption. Under similar conditions, one can find the manuscripts [1] and [8], which describe probabilistic swarm guidance algorithms. In [5], the authors present an approach to task allocation for a homogeneous swarm of robots. This is a Markov-chain based approach, where the goal is to converge to the desired population distribution over the set of tasks.

In the context of robotic swarms, programmable self-assembly of two-dimensional shapes with a thousand-robot swarm is demonstrated in [31]. These robots are capable of measuring distances to nearby neighbors which they use to localize themselves relative to other localized robots. Each robot then uses its position to implement an edge-following algorithm.

Another approach uses partial differential equations to model swarm behaviour, and control action is applied along the boundary of the swarm. Previous works on PDE-based methods with boundary control include [18], where the authors present an algorithm for the deployment of agents onto families of planar curves. Here, the swarm collective dynamics are modeled by the reaction-advection-diffusion PDE and the particular family of curves to which the swarm is controlled to is parametrized by the continuous agent identity in the interval of unit length. An extension of this work to deployment on a family of $2$D surfaces in $3$D space can be found in [29]. The problem of planning and task allocation is addressed in the framework of advection-diffusion-reaction PDEs in [14]. In [17] and [16], the authors present an optimal control problem formulation for swarm systems, where microscopic control laws are derived from the optimal macroscopic description using a potential function approach.

The problem of position-free extremum-seeking of an external scalar signal using a swarm of autonomous vehicles, inspired by bacterial chemotaxis, has been studied in [28].

In this work, we adopt a viewpoint outlined in [2], wherein we make an amorphous medium abstraction of the swarm, which is essentially a manifold with an agent located at each point. We then model the system using PDEs and design distributed control laws for them. An important component of this paper is the Laplacian-based distributed algorithm which we call pseudo-localization algorithm, which the agents implement to localize themselves in a new coordinate frame. The convergence properties of the graph Laplacian to the manifold Laplacian have been studied in [4], which find useful applications in this paper.

The main contribution of this paper is the development of distributed control laws for the index- and position-free density control of swarms to achieve general 1D and a large class of 2D density profiles. In very large swarms with thousands of agents, particularly those deployed indoors or at smaller scales, presupposing the availability of position information or pre-assignment of indices to individual agents would be a strong assumption. In this paper, in addition to not making the above assumptions, the agents are only capable of measuring the local density, and in the $2$D case, the density gradient and the normal direction to the boundary.

Under these assumptions, we present distributed pseudo-localization algorithms for one and two dimensions that agents implement to compute their position identifiers. Since every agent occupies a unique spatial position, we are able to rigorously characterize the resulting position assignment as a one-to-one correspondence between the set of spatial coordinates and the set of position identifiers, which corresponds to a diffeomorphism of the continuum domain. Based on this assignment, we then design control strategies for self-organization in one and two dimensions under the assumption that the motion control of agents is noiseless. The extension to the $2$D case leads to new difficulties related to the control of the swarm boundaries. To address these, we implement a variant of the $1$D pseudo-localization algorithm at the boundary during an initialization phase. A preliminary version of this work appeared in [23] where we presented an outline of the algorithms and stated some of the results. We develop them here rigorously, providing detailed proofs for our claims.

The paper is organized as follows. In Section 2, we introduce the basic notation and preliminary concepts used in the manuscript. We present the analysis of self-organization in one dimension in Section 4, where we introduce the pseudo-localization algorithm in Section 4.1 and the distributed control law in Section 4.2. After this, we generalize and extend the analysis for self-organization in two dimensions in Section 5. Section 6 contains numerical simulations of the results in the paper, and in Section 7, we present our conclusions.

2 Preliminaries

Let denote the set of all real numbers, ${\mathbb{R}}_{\geq 0}$ the set of non-negative real numbers, and n the $n$-dimensional Euclidean space. We use boldface letters to denote vectors in n. The norm $|\mathbf{x}|$ of a vector $\mathbf{x}\in^{n}$ is the standard Euclidean $2$-norm, unless otherwise specified. Let $\nabla=\left(\frac{\partial}{\partial x_{1}},\ldots\frac{\partial}{\partial x_{n}}\right)$ denote the gradient operator in n when acting on real-valued functions and the Jacobian in the context of vector-valued functions. As a shorthand, we let $\frac{\partial}{\partial z}(\cdot)=\partial_{z}(\cdot)$ for a variable $z$. Let $\Delta=\sum_{i=1}^{n}\frac{\partial^{2}}{\partial x_{i}^{2}}$ be the Laplace operator in n. We denote by either $\dot{S}$ or $\frac{dS}{dt}$ the total time derivative of $S(t)$. Given functions $f,g:\rightarrow$, we write $f=\mathcal{O}(g)$ if there exist positive constants $C$ and $c$ such that $|f(h)|\leq C|g(h)|$, for all $|h|\leq c$. Let $\mathcal{S}$ denote the set of agents in the swarm, and $N$ its cardinality. For the $1$D case, let $l\in\mathcal{S}$ denote the leftmost agent, and $r\in\mathcal{S}$ the rightmost one. Let $\mathcal{N}_{i}$ denote the spatial neighborhood of agent $i$, which comprises those agents located inside a small ball centered at $i$. A set-valued mapping, denoted by $f:\rightrightarrows^{2}$, maps the set of real numbers onto subsets of 2. For a bounded open set $\Omega\subset^{n}$, $\partial\Omega$ denotes its boundary, $\bar{\Omega}=\Omega\cup\partial\Omega$ its closure and $\mathring{\Omega}=\Omega\setminus\partial\Omega$ its interior with respect to the standard Euclidean topology. The set of smooth real-valued functions on $\Omega$ is denoted by $C^{\infty}(\Omega)$. We let $\mu$ (or $dx$ in 1D) denote the standard Lebesgue measure; with a slight abuse of notation, we sometimes omit $d\mu$ (resp. $dx$ in 1D) from long integrals. The Dirac measure $\delta$ on $\Omega$ defined for any $x\in\Omega$ and any measurable set $D\subseteq\Omega$ is given by $\delta_{x}(D)=1$ for $x\in D$, and $\delta_{x}(D)=0$ for $x\notin D$.

For two non-empty subsets $M_{1}$ and $M_{2}$ of a metric space $(M,d)$, the Hausdorff distance $d_{H}(M_{1},M_{2})$ between them is defined as:

[TABLE]

On a measurable space $U$, let $L^{p}(U)=\{f:U\rightarrow\,|\,\|f\|_{L^{p}(U)}=\left(\int_{U}|f|^{p}d\mu\right)^{1/p}<\infty\}$ constitute the $L^{p}$ space, where $\|\cdot\|_{L^{p}(U)}$ is the $L^{p}$ norm. Of particular interest is the $L^{2}$ space, or the space of square-integrable functions. In this paper, we denote by $\|f\|_{L^{2}(U)}$ the $L^{2}$ norm of $f$ with respect to the Lebesgue measure, and by $\|f\|_{L^{2}(U,\rho)}$ the weighted $L^{2}$ norm (with the strictly positive weight $\rho$ on $U$). The Sobolev space $W^{1,p}(U)$ over a measurable space $U$ is defined as $W^{1,p}(U)=\{f:U\rightarrow\,|\,\|f\|_{W^{1,p}}=\left(\int_{U}|f|^{p}+\int_{U}|\nabla f|^{p}\right)^{1/p}<\infty\}$. Of particular interest is the space $W^{1,2}$, also called the $H^{1}$ space. For two functions $f(t,\cdot)$ and $g(\cdot)$, we denote by $f\rightarrow_{L^{2}}g$ the convergence in $L^{2}$ norm (over the domain $U$ of the functions) of $f(t,\cdot)$ to $g(\cdot)$ as $t\rightarrow\infty$, that is, $\lim_{t\rightarrow\infty}\|f(t,\cdot)-g(\cdot)\|_{L^{2}}=0$. Convergence in $H^{1}$ norm is denoted similarly by $f\rightarrow_{H^{1}}g$.

We now state some well-known results that we will be used in the subsequent sections of this paper.

Lemma 2.1**.**

(Divergence Theorem [10]). For a smooth vector field $\mathbf{F}$ over a bounded open set $\Omega\subseteq^{n}$ with boundary $\partial\Omega$, the volume integral of the divergence $\nabla\cdot\mathbf{F}$ of $\mathbf{F}$ over $\Omega$ is equal to the surface integral of $\mathbf{F}$ over $\partial\Omega$:

[TABLE]

where $\mathbf{n}$ is the outward normal to the boundary and $dS$ the measure on the boundary. For a scalar field $U$ and a vector field $\mathbf{F}$ defined over $\Omega\subseteq^{n}$:

[TABLE]

Lemma 2.2**.**

*(Leibniz Integral Rule [10]).

Let $f\in\mathcal{C}^{\infty}(\times^{n})$ and $\Omega:\rightrightarrows^{n}$ be a smooth one-parameter family of bounded open sets in n generated by the flow corresponding to the smooth vector field $\mathbf{v}$ on n. Then:*

[TABLE]

Corollary 2.3**.**

*(Derivative of Energy Functional).

Let $U$ be an energy functional defined as follows:*

[TABLE]

for some function $f:\Omega\rightarrow$. Then,

[TABLE]

where $\frac{d}{dt}=\partial_{t}+\mathbf{v}\cdot\nabla$ is the total derivative.

Proof 2.4**.**

We have included the proof for this corollary for the sake of completeness. Using the Leibniz integral rule and the Divergence theorem, we have (it is understood that the integrations are with respect to the measure $\mu$):

[TABLE]

Lemma 2.5**.**

*(Poincaré-Wirtinger Inequality [26]).

For $p\in[1,\infty]$ and $\Omega$, a bounded connected open subset of n with a Lipschitz boundary, there exists a constant $C$ depending only on $\Omega$ and $p$ such that for every function $u$ in the Sobolev space $W^{1,p}(\Omega)$:*

[TABLE]

*where $u_{\Omega}=\frac{1}{|\Omega|}\int_{\Omega}ud\mu$, and $|\Omega|$ is the Lebesgue measure of $\Omega$. *

Lemma 2.6**.**

*(Rellich-Kondrachov Compactness Theorem [15]).

Let $U\subset^{n}$ be open, bounded and such that $\partial U$ is $C^{1}$. Suppose $1\leq p<n$, then $W^{1,p}(U)$ is compactly embedded in $L^{q}(U)$ for each $1\leq q<\frac{pn}{n-p}$. In particular, we have $W^{1,p}(U)$ is compactly contained in $L^{p}(U)$. *

Lemma 2.7**.**

*(LaSalle Invariance Principle [20, 34, 35]).

Let $\{\mathcal{P}(t)\,|\,t\in{\mathbb{R}}_{\geq 0}\}$ be a continuous semigroup of operators on a Banach space $U$ (closed subset of a Banach space with norm $\|\cdot\|$), and for any $u\in U$, define the positive orbit starting from $u$ at $t=0$ as $\Gamma_{+}(u)=\{\mathcal{P}(t)u\,|\,t\in{\mathbb{R}}_{\geq 0}\}\subseteq U$. Let $V:U\rightarrow$ be a continuous Lyapunov functional on $\mathcal{G}\subset U$ for $\mathcal{P}$ (such that $\dot{V}(u)=\frac{d}{dt}V(\mathcal{P}(t)u)\leq 0$ in $\mathcal{G}$). Define $E=\{u\in\bar{\mathcal{G}}\,|\,\dot{V}(u)=0\}$, and let $\tilde{E}$ be the largest invariant subset of $E$. If for $u_{0}\in\mathcal{G}$, the orbit $\Gamma_{+}(u_{0})$ is pre-compact (lies in a compact subset of $U$), then $\lim_{t\rightarrow+\infty}d_{U}(\mathcal{P}(t)u_{0},\tilde{E})=0$, where $d_{U}(y,\tilde{E})=\inf_{x\in\tilde{E}}\|y-x\|_{U}$ (where $d_{U}$ is the distance in $U$). *

2.1 Continuum model of the swarm

Given that $N$, the number of agents in the swarm, is very large, we will analyze the swarm dynamics through a continuum approximation. Let $t\in{\mathbb{R}}_{\geq 0}$, and let $M:\rightrightarrows^{n}$ be a smooth one-parameter family of bounded open sets, such that the agents are deployed over $\bar{M}(t)$ at time $t$. We denote by $\dot{\mathbf{r}}_{i}(t)=\mathbf{v}_{i}$, $\forall i\in\mathcal{S}$, where $\mathbf{r}_{i}(t)\in\bar{M}(t)$ is the position of the $i$th agent in the swarm at time $t$. Let $\rho:_{\geq 0}\times^{n}\rightarrow_{\geq 0}$ be the spatial density function supported on $\bar{M}(t)$ for all $t\geq 0$ (with $\rho(t,\mathbf{r})>0$ for $\mathbf{r}\in\bar{M}(t)$), such that $\int_{M(t)}\rho(t,\mathbf{r})d\mathbf{\mu}=1$. We assume that $M(t)$ is simply connected and that the boundary $\partial M(t)$ does not self-intersect for all $t\geq 0$.

Assuming that $\rho$ is smooth, the macroscopic dynamics can now be described by the continuity equation [10], assuming that the total number of agents is conserved:

[TABLE]

where $\mathbf{v}:_{\geq 0}\times^{n}\rightarrow^{n}$ is the velocity field with $\mathbf{v}_{i}(t)=\mathbf{v}(t,\mathbf{r}_{i})$, such that the one-parameter family $M$ is generated by the flow associated with $\mathbf{v}$.

2.2 Harmonic maps and diffeomorphisms

Let $(M,g)$ and $(N,h)$ be two Riemannian manifolds of dimensions $m$ and $n$, and Riemannian metrics $g$ and $h$, respectively. A map $\phi:M\rightarrow N$ is called harmonic if it minimizes the functional:

[TABLE]

where $dv_{g}$ is the Riemannian volume form on $M$. The Euler-Lagrange equation for the functional $E$, which also yields the minimum energy, is given by $\Delta\phi=0$, the Laplace equation [22]. It is useful to note that the solutions to the heat equation, in the limit $t\rightarrow\infty$, approach the harmonic map. This is proved later in Lemma 5.1, and forms the basis for the design of the distributed pseudo-localization algorithm. We now state a lemma on harmonic diffeomorphisms of Riemann surfaces (i.e., $m=n=2$ above).

Lemma 2.8**.**

*(Harmonic diffeomorphism [13]).

Let $(M,g)$ be a compact surface with boundary and $(N,h)$ a compact surface with non-positive curvature. Suppose that $\psi:M\rightarrow N$ is a diffeomorphism onto $\psi(M)$. Assume that $\psi(M)$ is convex. Then there is a unique harmonic map $\phi:M\rightarrow N$ with $\phi=\psi$ on $\partial M$, such that $\phi:{M}\rightarrow\phi(M)$ is a diffeomorphism. *

We note that the non-positive curvature constraint in the lemma is essentially a constraint on the metric $h$ on $N$, and the curvature is zero for the Euclidean metric.

3 Problem description and conceptual approach

In this section, we provide a high-level description of the proposed problem and explain the conceptual idea behind our approach. The technical details can be found in the following sections.

The problem at hand is to ultimately design a distributed control law for a swarm to converge to a desired configuration. Here, a swarm configuration is a density function $\rho$ of the multi-agent system and the objective is that agents reconfigure themselves into a desired known density $\rho^{*}$. To do this, an agent at position $x$ is able to measure the current local density value, $\rho(t,x)$; however, its position $x$ within the swarm is unknown. Thus, given $\rho^{*}$, an agent at $x$ cannot directly compute $\rho^{*}(x)$ nor a feedback law based on $\rho-\rho^{*}$. To solve this problem, we devise a mechanism that allows agents to determine their coordinates in a distributed way in an equivalent coordinate system.

Note that, given a diffeomorphism $\Theta^{*}$ from the spatial domain of the swarm onto the unit interval or disk (i.e. a coordinate transformation), we can equivalently provide the agents with a transformed density function $p^{*}$, such that $p^{*}=\rho^{*}\circ(\Theta^{*})^{-1}$. In this way, instead of $\rho^{*}$ the agents are given $p^{*}$, but still do not have access to $\Theta^{*}$. The pseudo-localization algorithm is a mechanism that agents employ to progressively compute an appropriate (configuration-dependent) diffeomorphism by local interactions.

In 1D, the pseudo-localization algorithm is a continuous-time PDE system in a new variable or pseudo-coordinate $X$ which plays the role of an “approximate $x$ coordinate” that agents can use to know where they are. The input to this system is the current density value $\rho$, see Figure 1 for an illustration, and the objective is that $X$ converges to a $\rho$-dependent diffeomorphism. On the other hand, the variable $X$ and the function $p^{*}$ are used to define the control input of another PDE system in the density $\rho$. In this way, we have a feedback interconnection of two systems, one in $X$ and one in $\rho$, with the goal to achieve $X\rightarrow\Theta^{*}$ (the pseudo-coordinate $X$ converges to a true coordinate given by $\Theta^{*}$) and $\rho\rightarrow\rho^{*}$.

As for the control design methodology, we follow a constructive, Lyapunov-based approach to designing distributed control laws for the swarm dynamics modeled by PDEs. For this, we define appropriate non-negative energy functionals that encode the objective and choose control laws that keep the time derivative of the energy functional non-positive. This, along with well-known results on the precompactness of solutions as in Lemma 2.6, the Rellich Kondrachov compactness theorem, allows us to apply the LaSalle Invariance Principle in Lemma 2.7 and other technical arguments to establish the convergence results that we seek.

In the 1D case, we can identify a set of diffeomorphisms $\Theta$ associated with any $\rho$ that eventually converge to $\Theta^{*}$, and simultaneously control boundary agents into a desired final domain (the support of $\rho^{*}$). These are given by the cumulative distribution function associated with the density function; see Section 4.1. The 2D case is more complex, and analogous results could not be derived in their full generality. Unlike the 1D case, estimating the cumulative distribution is not straightforward in the 2D case. Instead, we set out to find diffeomorphisms as the result of a distributed algorithm. Given that the discretization of heat flow naturally leads to distributed algorithms, we investigate under what conditions this is the case via harmonic map theory. On the control side, there also are additional difficulties, and because of this, we simplify the control strategy into three stages. In the first stage, the boundary agents are re-positioned onto the boundary of the desired domain while containing the others in the interior. Once this is achieved, the second and third stages can be seen again as the interconnection of two systems in pseudo-coordinates $R=(X,Y)$ (instead of $X$) and $\rho$, analogously to Figure 1. However, we apply a two time-scale separation for analysis by which coordinates are computed in a fast-time scale and reconfiguration is done in a slow-time scale, which allows for a sequential analysis of the two stages. We then study the robustness of this approach. We refer the reader to the extended version of this paper [24] for further description of the discrete implementation.

4 Self-organization in one dimension

In this section, we present our proposed pseudo-localization algorithm and the distributed control law for the $1$D self-organization problem.

For each $t\in{\mathbb{R}}_{\geq 0}$, let $M(t)=(0,L(t))\subset$ be the interval (with boundary $\{0,L(t)\}$) in which the agents are distributed in 1D, and let $\rho:\times\rightarrow_{\geq 0}$ be the normalized density function supported on $\bar{M}(t)$, for all $t\geq 0$ (with $\rho(t,x)>0$, $\forall x\in\bar{M}(t)$), describing the swarm on that interval. Without loss of generality, we place the origin at the leftmost agent of the swarm. We also assume that the leftmost and the rightmost agents, $l$ and $r$, are aware that they are at the boundary. Let $\rho^{*}:\bar{M}^{*}=[0,L^{*}]\rightarrow{\mathbb{R}}_{>0}$ be the desired normalized density distribution.

Since a direct feedback control law can not be implemented by agents because they do not have access to their positions, we introduce an equivalent representation of the density $\rho^{*}$, $p^{*}$, depending on a particular diffeomorphism $\Theta^{*}$. First, define $\Theta^{*}:\bar{M}^{*}\rightarrow[0,1]$ such that $\Theta^{*}(x)=\int_{0}^{x}\rho^{*}(\bar{x})d\bar{x}$ and $\Theta^{*}(L^{*})=1$.

Now, let $p^{*}:[0,1]\rightarrow{\mathbb{R}}_{>0}$, and $\theta^{*}\in\Theta^{*}(\bar{M}^{*})=[0,1]$, be such that $p^{*}(\theta^{*})=\rho^{*}((\Theta^{*})^{-1}(\theta^{*}))=\rho^{*}(x)$.

The function $p^{*}$, which represents the desired density distribution mapped onto the unit interval $[0,1]$, is computed offline and is broadcasted to the agents prior to the beginning of the self-organization process. We use $p^{*}$ to derive the distributed control law which the agents implement. We assume that $p^{*}$ is a Lipschitz function in the sequel. {assumption}(Uniform boundedness of density function). We assume that the density function and its derivative are uniformly bounded in its support, that is, for $\rho(t,\cdot)$ and $\partial_{x}\rho(t,\cdot)$ there exist uniform lower bounds $d_{l},D_{l}$ and uniform upper bounds $d_{u},D_{u}$ (where $0<d_{l}\leq d_{u}<\infty$ and $0<D_{l}\leq D_{u}<\infty$) (that is, $d_{l}\leq\rho(t,x)\leq d_{u}$ for all $t\in{\mathbb{R}}_{\geq 0}$ and $x\in[0,L(t)]$ and $D_{l}\leq\partial_{x}\rho(t,x)\leq D_{u}$ for all $t\in{\mathbb{R}}_{\geq 0}$ and $x\in(0,L(t))$).

4.1 Pseudo-localization algorithm in one dimension

We first consider the static case, that is, the design of the pseudo-localization dynamics on $X$ of the upper block in Figure 1, when the agents and $\rho$ are stationary. We define $\Theta:\bar{M}=[0,L]\rightarrow[0,1]$ as:

[TABLE]

such that $\Theta(L)=1$. In other words, $\Theta$ is the cumulative distribution function (CDF) associated with $\rho$. (Note that the domains are static and hence the argument $t$ has been dropped, which will be reintroduced later.)

Lemma 4.1**.**

*(The CDF diffeomorphism).

Given $\rho:\bar{M}\rightarrow{\mathbb{R}}_{>0}$, a $C^{1}$ function, the mapping $\Theta:\bar{M}\rightarrow[0,1]$ as defined above, is a diffeomorphism and $\Theta(\bar{M})=[0,1]$. *

Proof 4.2**.**

*Since $\rho(x)>0$, $\forall x\in\bar{M}$, it follows that $\Theta$ is a strictly increasing function of $x$, and is therefore a one-to-one correspondence on $\bar{M}$. Moreover, $\Theta$ is atleast $C^{1}$ and has a differentiable inverse, which implies it is a diffeomorphism. Finally, since $\Theta(L)=1$, we have $\Theta(\bar{M})=[0,1]$. *

Our goal here is to set up a partial differential equation with appropriate boundary conditions that yield the diffeomorphism $\Theta$ as its asymptotically stable steady-state solution. We begin by setting up the pseudo-localization dynamics for a stationary swarm (for which the spatial domain $M$ and the density distribution $\rho$ are fixed). Let $X:\times\bar{M}\rightarrow$ be such that $(t,x)\mapsto X(t,x)\in$, with:

[TABLE]

where $\alpha:\rightarrow$ is a control input at the boundary $x=0$ and $\beta:\rightarrow$ is a control input at the boundary $x=L$. From (5), we observe that $\partial_{x}\left(\frac{\partial_{x}\Theta}{\rho}\right)=0$. Letting $w=X-\Theta$ denote the error, we obtain:

[TABLE]

{assumption}

(Well-posedness of the pseudo-localization dynamics).

We assume that the pseudo-localization dynamics (6) (and (7)) is well-posed, that the solution is sufficiently smooth (at least $\mathcal{C}^{2}$ in the spatial variable, even as $t\rightarrow\infty$) and belong to the Sobolev space $H^{1}(M)$.

Lemma 4.3**.**

*(Pointwise convergence to diffeomorphism). Under Assumption 4.1, on the well-posedness of the pseudo-localization dynamics, and Assumption 4 on the boundedness of $\rho$, the solutions to PDE (6) converge pointwise to the CDF diffeomorphism $\Theta$ defined in (5), as $t\rightarrow\infty$, for all $C^{2}$ initial conditions $X_{0}$. *

In this case, the swarm is stationary, which implies that the distribution $\rho$ is fixed (and so is its support $\bar{M}$), and the uniform boundedness assumption 4 simply becomes a boundedness assumption.

Proof 4.4**.**

We prove that the solutions to the PDE (6) converge pointwise to the diffeomorphism $\Theta$ by showing that $w\rightarrow 0$, as $t\rightarrow\infty$, pointwise for (7). For this, we consider a functional $V$, given by (integrations are with respect to the Lebesgue measure):

[TABLE]

The time derivative $\dot{V}$ is given by:

[TABLE]

Here, replace $\partial_{t}w$ in the first integral with the dynamics in (7), and then use $\partial_{t}\partial_{x}=\partial_{x}\partial_{t}$ in the second integral together with the Divergence Theorem in Lemma 2.1. We obtain:

[TABLE]

(After the second equal sign, apply again the Divergence Theorem on the first integral of the previous line, and replace $\partial_{t}w$ from (7).) Substituting from (7), we have:

[TABLE]

*Clearly, $\dot{V}\leq 0$, and $w(t,\cdot)\in H^{1}(M)$, for all $t$. Moreover, since $V(t)\leq V(0)$ and since $\rho$ is uniformly bounded according to Assumption 4, we have that $w(t,\cdot)$ is bounded in $H^{1}(M)$. Moreover, by the Rellich-Kondrachov Theorem of Lemma 2.6, $H^{1}(M)$ is compactly contained in $L^{2}(M)$. Then it follows that the solutions $w(t,\cdot)$ are precompact. Thus, by the LaSalle Invariance Principle of Lemma 2.7, the solution to (7) converges in $L^{2}$-norm to the largest invariant subset of $\dot{V}^{-1}(0)$. Note that $\dot{V}=0$ implies $\int_{M}\frac{1}{\rho}|\partial_{x}w|^{2}=0$. Thus, $\lim_{t\rightarrow\infty}\int_{M}\frac{1}{\rho}|\partial_{x}w|^{2}=0$. Since $\rho$ is bounded ($\sup\rho<\infty$), we have $\lim_{t\rightarrow\infty}\frac{1}{\sup\rho}\int_{M}|\partial_{x}w|^{2}\leq\lim_{t\rightarrow\infty}\int_{M}\frac{1}{\rho}|\partial_{x}w|^{2}=0$, which implies $\lim_{t\rightarrow\infty}\int_{M}|\partial_{x}w|^{2}=\lim_{t\rightarrow\infty}\|\partial_{x}w\|_{L^{2}(M)}^{2}=0$. Now, $\lim_{t\rightarrow\infty}|w(t,x)|=\lim_{t\rightarrow\infty}|w(t,0)+\int_{0}^{x}\partial_{x}w(t,\cdot)|\leq\lim_{t\rightarrow\infty}|w(t,0)|+\int_{0}^{x}|\partial_{x}w(t,\cdot)|\leq\lim_{t\rightarrow\infty}|w(t,0)|+\sqrt{L(t)}\|\partial_{x}w(t,\cdot)\|_{L^{2}(M)}=0$ (since $\lim_{t\rightarrow\infty}w(t,0)=0$ and $\lim_{t\rightarrow\infty}\|\partial_{x}w(t,\cdot)\|_{L^{2}(M)}=0$). Thus, $\lim_{t\rightarrow\infty}w(t,x)=0$, for all $x\in M$. Therefore, the solutions to (7) converge to $w\equiv 0$ pointwise, as $t\rightarrow\infty$, from any smooth initial $w_{0}=X_{0}-\Theta$. *

We now have that the solution to the pseudo-localization dynamics converges to the diffeomorphism $\Theta$ in the stationary case. For the dynamic case, we modify (6) to account for agent motion. Let $X:\times\rightarrow$ be supported on $\bar{M}(t)=[0,L(t)]$ for all $t\geq 0$. Using the relation $\frac{dX}{dt}=\partial_{t}X+v\partial_{x}X$, where $v$ is the velocity field on the spatial domain, we consider:

[TABLE]

In the dynamic case, and w.l.o.g. we have set $\alpha(t)=0$ for all $t\geq 0$, for simplicity. We will use the above PDE system in the design of the distributed motion control law, redesigning the boundary control $\beta$ to achieve convergence of the entire system. We now discretize (8) to obtain a distributed pseudo-localization algorithm. Let $X_{i}(t)=X(t,x_{i})$, where $x_{i}\in\bar{M}(t)$ is the position of the $i^{\textup{th}}$ agent. We identify the agent $i$ with its desired coordinate in the unit interval at time $t$, i.e., $\Theta(t,x)=\theta\in[0,1]$, where $\Theta(t,x)=\int_{0}^{x}\rho(t,\bar{x})d\bar{x}$ from (5), which now shows the time dependency of $\rho$. In this way, $\rho(t,x)=\partial_{x}\Theta(t,x)$. It follows that $\partial_{x}(\cdot)=\partial_{\theta}(\cdot)\partial_{x}\theta=\partial_{\theta}(\cdot)\rho$. Therefore, $\frac{1}{\rho}\partial_{x}(\cdot)=\partial_{\theta}(\cdot)$. From (8), we have:

[TABLE]

Now, we discretize (9) with the consistent finite differences $\frac{dX}{dt}\approx\frac{X_{i}(t+1)-X_{i}(t)}{\Delta t}$ and $\frac{\partial^{2}X}{\partial\theta^{2}}\approx\frac{X_{i+1}-2X_{i}+X_{i-1}}{(\Delta\theta)^{2}}$ (that is, we have that $\lim_{\Delta t\rightarrow 0}\frac{X_{i}(t+1)-X_{i}(t)}{\Delta t}=\frac{dX}{dt}$ and that $\lim_{\Delta\theta\rightarrow 0}\frac{X_{i+1}-2X_{i}+X_{i-1}}{(\Delta\theta)^{2}}=$ $\frac{\partial^{2}X}{\partial\theta^{2}}$). Now, with the choice $3\Delta t=(\Delta\theta)^{2}$, and from (8), we obtain for $i\in\mathcal{S}\setminus\left\{l,r\right\}$:

[TABLE]

Equation (10) is the discrete pseudo-localization algorithm to be implemented synchronously by the agents in the swarm, starting from any initial condition $X_{0}$. The leftmost agent holds its value at zero while the rightmost agent implements the boundary control $\beta$. In the following section we analyze its behavior together with that of the dynamics on $\rho$.

4.2 Distributed density control law and analysis

In this subsection, we propose a distributed feedback control law to achieve $\rho\rightarrow\rho^{*}$ and $w\rightarrow 0$, as $t\rightarrow\infty$, through a distributed control input $v$ and a boundary control $\beta$. We refer the reader to [25] for an overview of Lyapunov-based methods for stability analysis of PDE systems.

From (3) and (8), we have the dynamics:

[TABLE]

This realizes the feedback interconnection of Figure 1. {assumption}(Well-posedness of the full PDE system).

We assume that (11) is well posed, and that the solutions $(\rho(t,\cdot),X(t,\cdot))$ are sufficiently smooth (both in $t$ and $x\in[0,L(t)]$), satisfy Assumption 4 on the uniform boundedness of $\rho$ and $\partial_{x}\rho$, and are bounded in the Sobolev space $H^{1}((0,1/d_{l}))$.

We also assume that the agent at position $x$ at time $t$ is able to measure $\rho(t,x)$. However, the agents in the swarm do not have access to their positions, and therefore cannot access $\rho^{*}(x)$, which could be used to construct a feedback law. To circumvent this problem, we propose a scheme in which the agents use the position identifier or pseudo-localization variable $X$ to compute $p^{*}\circ X(t,x)$, using this as their dynamic set-point. The idea is to then design a distributed control law and a boundary control law such that $\rho\rightarrow p^{*}\circ X$ and $X\rightarrow\Theta^{*}$, as $t\rightarrow\infty$, to obtain $\rho\rightarrow p^{*}\circ\Theta^{*}=\rho^{*}$. Recall that the function $p^{*}$ is computed offline and is broadcasted to the agents prior to the beginning of the self-organization process, and that $p^{*}$ is assumed to be a Lipschitz function. Consider the distributed control law, defined as follows for all time $t$:

[TABLE]

together with the boundary control law:

[TABLE]

We remark again that the agents implementing the control laws (12) and (13) do not require position information, because for the agent at position $x$ at time $t$, $\rho(t,x)$ is a measurement, $X(t,x)$ is the pseudo-localization variable, through which $p^{*}\circ X(t,x)$ can be computed.

Theorem 4.5**.**

*(Convergence of solutions). Under the well-posedness Assumption 4.2, the solutions $(\rho(t,\cdot),X(t,\cdot))$ to (11), under the control laws (12) and (13), converge to $(\rho^{*},\Theta^{*})$, $\rho\rightarrow\rho^{*}$ in $L^{2}-norm$ and $X\rightarrow\Theta^{*}$ pointwise as $t\rightarrow\infty$, from any smooth initial condition $(\rho_{0}$, $X_{0})$. *

Proof 4.6**.**

Consider the candidate control Lyapunov functional $V$:

[TABLE]

Taking the time derivative of $V$ along the dynamics (11), using Lemma 2.2 on the Leibniz integral rule, and applying Corollary 2.3 on the derivative of energy functionals, we obtain:

[TABLE]

*Now, $\frac{d\rho}{dt}=\partial_{t}\rho+v\partial_{x}\rho=-\rho\partial_{x}v$ (since $\partial_{t}\rho=-\partial_{x}(\rho v)$, from (11)), and $\partial_{t}w=\frac{1}{\rho}\partial_{x}\left(\frac{\partial_{x}w}{\rho}\right)-v\partial_{x}w$. Thus, we obtain: *

[TABLE]

Now, using the above equation, applying the Divergence theorem (2) (integration by parts) and rearranging the terms, we obtain:

[TABLE]

Since $w(0)=0$, the above equation reduces to:

[TABLE]

From (12) and (13), we have $\partial_{x}v=(\rho-p^{*}\circ X)-\frac{\partial_{X}p^{*}}{\rho(\rho+p^{*}\circ X)}\partial_{x}\left(\frac{\partial_{x}X}{\rho}\right)$, and \frac{dw}{dt}\bigg{|}_{L(t)}=-\left(\frac{\partial_{x}w}{\rho}+kw\right)\bigg{|}_{L(t)}, and we obtain:

[TABLE]

Clearly, $\dot{V}\leq 0$, and $\rho(t,\cdot),w(t,.)\in H^{1}((0,1/d_{l}))$, for all $t$. By Lemma 2.6, the Rellich-Kondrachov Compactness Theorem, the space $H^{1}((0,1/d_{l}))$ is compactly contained in $L^{2}((0,1/d_{l}))$, and the bounded solutions (by Assumption 4.2) in $H^{1}((0,1/d_{l}))$ are then precompact in $L^{2}((0,1/d_{l}))$. Moreover, the set of $(\rho,X)$ satisfying Assumption 4.2 is dense in $L^{2}((0,1/d_{l}))$. Then, by the LaSalle Invariance Principle, Lemma 2.7, we have that the solutions to (11) converge in the $L^{2}$-norm to the largest invariant subset of $\dot{V}^{-1}(0)$. This implies that:

[TABLE]

*Thus, we have: *

[TABLE]

Using the Poincaré-Wirtinger inequality, Lemma 2.5, again, we note that this implies $\lim_{t\rightarrow\infty}\|w-\int_{0}^{L(t)}w\|_{L^{2}((0,L(t)))}=0$. We have $\lim_{t\rightarrow\infty}|\int_{0}^{L(t)}w|=|\int_{0}^{L(t)}\int_{0}^{x}\partial_{x}w|\leq L(t)^{3/2}\|\partial_{x}w\|_{L^{2}((0,L(t)))}=0$, which implies that $\lim_{t\rightarrow\infty}\int_{0}^{L(t)}w=0$ and therefore $\lim_{t\rightarrow\infty}\|w\|_{L^{2}((0,L(t)))}=0$. Thus, we get $\lim_{t\rightarrow\infty}\|w(t,\cdot)\|_{H^{1}((0,L(t)))}=0$, or in other words, $w\rightarrow_{H^{1}}0$. Now, $\lim_{t\rightarrow\infty}|w(t,x)|=\lim_{t\rightarrow\infty}|w(t,0)+\int_{0}^{x}\partial_{x}w(t,\cdot)|\leq\lim_{t\rightarrow\infty}|w(t,0)|+\int_{0}^{x}|\partial_{x}w(t,\cdot)|\leq\lim_{t\rightarrow\infty}|w(t,0)|+\sqrt{L(t)}\|w(t,\cdot)\|_{H^{1}((0,L(t)))}=0$, which implies that $w\rightarrow 0$ pointwise. Given that $w=X-\Theta$, we have $\lim_{t\rightarrow\infty}X(t,\cdot)-\Theta(t,\cdot)=0$. Let $\lim_{t\rightarrow\infty}L(t)=L$ and $\lim_{t\rightarrow\infty}\Theta(t,\cdot)=\bar{\Theta}(\cdot)$, which implies that $X\rightarrow\bar{\Theta}$ pointwise.

From the above, we have $\lim_{t\rightarrow\infty}\|\rho(t,\cdot)-p^{*}\circ\bar{\Theta}\|_{L^{2}((0,L(t)))}=\lim_{t\rightarrow\infty}\|\rho(t,\cdot)-p^{*}\circ X(t,\cdot)+p^{*}\circ X(t,\cdot)-p^{*}\circ\bar{\Theta}\|_{L^{2}((0,L(t)))}\leq\lim_{t\rightarrow\infty}\|\rho(t,\cdot)-p^{*}\circ X(t,\cdot)\|_{L^{2}((0,L(t)))}+\|p^{*}\circ X(t,\cdot)-p^{*}\circ\bar{\Theta}\|_{L^{2}((0,L(t)))}=0$ (this follows from the assumption that $p^{*}$ is Lipschitz, since $\|p^{*}\circ X-p^{*}\circ\bar{\Theta}\|_{L^{2}}\leq c\|X-\bar{\Theta}\|_{L^{2}}$ for some Lipschitz constant $c$). Thus, we have $\rho\rightarrow_{L^{2}}p^{*}\circ\bar{\Theta}$.

Now, we are interested in the limit density distribution $\bar{\rho}=p^{*}\circ\bar{\Theta}$, and by the definition of $\bar{\Theta}$ we have $\bar{\Theta}(x)=\int_{0}^{x}\bar{\rho}$. We now prove that this limit $(\bar{\rho},\bar{\Theta})$ is unique, and that $(\bar{\rho},\bar{\Theta})=(\rho^{*},\Theta^{*})$. From the definition of $\bar{\Theta}$, we get $\frac{d\bar{\Theta}}{dx}(x)=\bar{\rho}(x)=p^{*}(\bar{\Theta}(x))>0$, $\forall\bar{\Theta}(x)\in[0,1]$. We therefore have:

[TABLE]

Recall from the definition of $p^{*}$ and (5) that $p^{*}\circ\Theta^{*}(x)=\rho^{*}(x)$, and $\frac{d}{dx}\Theta^{*}(x)=\rho^{*}(x)=p^{*}\circ\Theta^{*}(x)$, which implies that $\frac{d\Theta^{*}}{dx}=p^{*}(\theta^{*})>0$, where $\theta^{*}=\Theta^{*}(x)$. Therefore:

[TABLE]

From the above two equations, we get:

[TABLE]

*for all $x$, and since $p^{*}$ is strictly positive, it implies that $\bar{\Theta}=\Theta^{*}$, and we obtain $\bar{\rho}=p^{*}\circ\bar{\Theta}=p^{*}\circ\Theta^{*}=\rho^{*}$. And we know that $\rho\rightarrow_{L^{2}}p^{*}\circ\bar{\Theta}=p^{*}\circ\Theta^{*}=\rho^{*}$. In other words, $\rho$ converges to $\rho^{*}$ in the $L^{2}$ norm. *

4.2.1 Physical interpretation of the density control law

For a physical interpretation of the control law, we first rewrite some of the terms in a suitable form. From (11), we know that:

[TABLE]

The second term in the expression for $\partial_{x}v$ in the law (12) can thus be rewritten as:

[TABLE]

Now, from above and (12), we obtain:

[TABLE]

Equation (15) gives the velocity of the agent at $x$ at time $t$. Now, to interpret it, we first consider the case where the pseudo-localization error is zero, that is, when $X=\Theta^{*}$. This would imply that $p^{*}\circ X=p^{*}\circ\Theta^{*}=\rho^{*}$, $\frac{dX}{dt}=\frac{d\Theta^{*}}{dt}=0$, and we obtain:

[TABLE]

The term $\int_{0}^{x}(\rho-\rho^{*})=\int_{0}^{x}\rho-\int_{0}^{x}\rho^{*}$ is the difference between the number of agents in the interval $[0,x]$ and the desired number of agents in $[0,x]$. If the term is positive, it implies that there are more than the desired number of agents in $[0,x]$ and the control law essentially exerts a pressure on the agent to move right thereby trying to reduce the concentration of agents in the interval $[0,x]$, and, vice versa, when the term is negative. This eventually accomplishes the desired distribution of agents over a given interval. This would be the physical interpretation of the control law for the case where the pseudo-localization error is zero (that is, the agents have full information of their positions).

However, in the transient case when the agents do not possess full information of their positions and are implementing the pseudo-localization algorithm for that purpose, the control law requires a correction term that accounts for the fact that the transient pseudo coordinates $X(t,x)$ cannot be completely relied upon. This is what the second term $\int_{0}^{x}\frac{1}{(\rho+p^{*}\circ X)}\frac{dp^{*}}{dt}$ in (15) corrects for. When this term is positive, that is, $\int_{0}^{x}\frac{1}{(\rho+p^{*}\circ X)}\frac{dp^{*}}{dt}>0$, it roughly implies that the “estimate” of the desired number of agents in the interval $[0,x]$ is increasing (indicating that an increase in the concentration of agents in $[0,x]$ is desirable), and the term essentially reduces the “rightward pressure” on the agent (note that this term will have a negative contribution to the velocity (15)).

4.3 Discrete implementation

In this section, we present a scheme to compute $p^{*}$ (the transformed desired density profile) and a consistent discretization scheme for the distributed control law. We follow that up with a discussion on the convergence of the discretized system and a pseudo-code for the implementation.

4.3.1 On the computation of $p^{*}$

We now provide a method for computing $p^{*}$ from a given $\rho^{*}$ via interpolation. Let the desired domain $M^{*}=[0,L^{*}]$ be discretized uniformly to obtain $M^{*}_{d}=\{0=x_{1},\ldots,x_{m}=L^{*}\}$ such that $x_{j}-x_{j-1}=h$ (constant step-size). Note that $m$ is the number of interpolation points, not equal to the number of agents. The desired density $\rho^{*}:[0,L^{*}]\rightarrow{\mathbb{R}}_{>0}$ is known, and we compute the value of $\rho^{*}$ on $M^{*}_{d}$ to get $\rho^{*}(x_{1},\ldots,x_{m})=(\rho^{*}_{1},\ldots,\rho^{*}_{m})$. We also have $\Theta^{*}(x)=\int_{0}^{x}\rho^{*}d\mu$, for all $x\in[0,L^{*}]$. Now, computing the integral with respect to the Dirac measure for the set $M^{*}_{d}$, we obtain $\Theta_{d}^{*}(x_{1},\ldots,x_{m})=(\theta^{*}_{1},\ldots,\theta^{*}_{m})$, where $\theta^{*}_{1}=0$ and $\theta^{*}_{k}=\frac{1}{2}\sum_{j=1}^{k}(\rho^{*}_{j-1}+\rho^{*}_{j})h$, for $k=2,\ldots,m$ (note that $0=\theta^{*}_{1}\leq\theta^{*}_{2}\leq\ldots\leq\theta^{*}_{m}\leq 1$ and $\lim_{h\rightarrow 0}\theta^{*}_{m}=\Theta^{*}(L^{*})=1$). Now, the value of the function $p^{*}$ at any $X\in[0,1]$ can be now obtained from the relation $p^{*}(\theta^{*}_{k})=\rho^{*}_{k}$, for $k=1,\ldots,m$, by an appropriate interpolation.

4.3.2 Discrete control law

A discretized pseudo-localization algorithm is given by (10). We now discretize (12) to obtain an implementable control law for a finite number of agents $i\in\mathcal{S}$, and a numerical simulation of this law is later presented in Section 6.

Let $i\in\mathcal{S}\setminus\{l,r\}$. First note that \partial_{x}v=(\partial_{\theta}v)\bigg{|}_{\theta=\Theta(x)}(\partial_{x}\Theta)=(\partial_{\theta}v)\bigg{|}_{\theta=\Theta(x)}\rho (where $v\equiv v(\Theta(x))$). Using a consistent backward differencing approximation, and recalling that $\Delta\theta=\epsilon$, we can write:

[TABLE]

where $\rho_{i}$ is agent $i$’s density measurement.

From Section 4.1, recall the consistent finite-difference approximation:

[TABLE]

With $\kappa=\frac{1}{2\epsilon}$, from (12) and the above equation, we obtain the law for agent $i$ as:

[TABLE]

with $v_{l}=0$. The computation in $v$ can be implemented by propagating from the leftmost agent to the rightmost agent along a line graph $\mathcal{G}_{line}$ (with message receipt acknowledgment). Note that this propagation can alternatively be formulated by each agent averaging appropriate variables with left and right neighbors, which will result in a process similar to a finite-time consensus algorithm. Now, the boundary control (13) is discretized (with $\partial_{t}\beta\approx\frac{\beta(t+1)-\beta(t)}{\Delta t}$), with the choice $k=\frac{1}{\epsilon}$ to:

[TABLE]

4.3.3 On the convergence of the discrete system

The discretized pseudo-localization algorithm (10) with the boundary control law (13), can be rewritten as:

[TABLE]

where $X(t)=(X_{l}(t),\ldots,X_{r}(t))$, $L$ is the Laplacian of the line graph $\mathcal{G}_{line}$ and the input $u(t)=\left(0,\ldots,0,\frac{\epsilon}{3}(2-\beta(t))\right)$. This discretized system is stable and we thereby have that the discretized pseudo-localization algorithm is consistent and stable. Thus, by the Lax Equivalence Theorem [33], the solution of (19) converges to the solution of (8) with the boundary control (13) as $N\rightarrow\infty$. Due to the nonlinear nature of the discrete implementation of the equation in $\rho$, we are only certain that we have a consistent discrete implementation in this case (no similar convergence theorem exists for discrete approximations of nonlinear PDEs.)

5 Self-organization in two dimensions

In this section, we present the two-dimensional self-organization problem. Although our approach to the  $2$D problem is fundamentally similar to the  $1$D case, we encounter a problem in the two-dimensional case that did not require consideration in one dimension, and it is the need to control the shape of the spatial domain in which the agents are distributed. We overcome this problem by controlling the shape of the domain with the agents on the boundary, while controlling the density distribution of the agents in the interior.

Let $M:\rightrightarrows^{2}$ be a smooth one-parameter family of bounded open subsets of 2, such that $\bar{M}(t)$ is the spatial domain in which the agents are distributed at time $t\geq 0$. Let $\rho:\times^{2}\rightarrow_{\geq 0}$ be the spatial density function with support $\bar{M}(t)$ for all $t\geq 0$; that is, $\rho(t,x)>0$, $\forall\,x\in\bar{M}(t)$, and $t\geq 0$. Without loss of generality, we shift the origin to a point on the boundary of the family of domains, such that $(0,0)\in\partial M(t)$, for all $t$. Let $\rho^{*}:M^{*}\rightarrow{\mathbb{R}}_{>0}$ be the desired density distribution, where $M^{*}$ is the target spatial domain. From here on, we view $\bar{M}$ as a one-parameter family of compact $2$-submanifolds with boundary of 2. Just as in the $1$D case, the agents do no have access to their positions but know the true $x$- and $y$-directions.

In what follows we present our strategy to solve this problem, which we divide into three stages for simplicity of presentation and analysis. In the first stage, the agents converge to the target spatial domain $M^{*}$ with the boundary agents controlling the shape of the domain. In stage two, the agents implement the pseudo-localization algorithm to compute the coordinate transformation. In the third stage, the boundary agents remain stationary and the agents in the interior converge to the desired density distribution. This simplification is performed under the assumption that, once the agents have localized themselves at a given time, they can accurately update this information by integrating their (noiseless) velocity inputs. Noisy measurements would require that these phases are rerun with some frequency; e.g. using fast and slow time scales as described in Section 3.

5.1 Pseudo-localization algorithm for boundary agents

To begin with, we propose a pseudo-localization algorithm for the boundary agents which allows for their control in the first stage. To do this, we assume that the agents have a boundary detection capability (can approximate the normal to the boundary), the ability to communicate with neighbors immediately on either side along the boundary curve, and can measure the density of boundary agents.

Let $M_{0}\subset^{2}$ be a compact $2$-manifold with boundary $\partial M_{0}$ and let $(0,0)\in\partial M_{0}$. To localize themselves, the agents on $\partial M_{0}$ implement the distributed $1$D pseudo-localization algorithm presented in Section 4.1. This yields a parametrization of the boundary $\Gamma:\partial M_{0}\rightarrow[0,1)$, with $\Gamma(0,0)=0$, such that the closed curve which is the boundary $\partial M_{0}$ is identified with the interval $[0,1)$. We have that, for $\gamma\in[0,1)$, $\Gamma^{-1}(\gamma)\in\partial M_{0}$. For $\gamma\in[0,1)$, let $s(\gamma)$ be the arc length of the curve $\partial M_{0}$ from the origin, such that $s(0)=0$ and $\lim_{\gamma\rightarrow 1}s(\gamma)=l$. We assume that the boundary agents have access to the unit outward normal $\mathbf{n}(\gamma)$ to the boundary, and thus the unit tangent $\mathbf{s}(\gamma)$.

Let $q:[0,l)\rightarrow{\mathbb{R}}_{>0}$ denote the normalized density of agents on the boundary, such that we have $\int_{0}^{l}q(s)ds=1$. Now the 1D pseudo-localization algorithm of Section 4.1 serves to provide a 2D boundary pseudo-localization as follows. Note that $\frac{ds}{d\gamma}=\frac{1}{q(\gamma)}$, and $(dx,dy)=\mathbf{s}ds$, which implies $(dx,dy)=\frac{1}{q(\gamma)}\mathbf{s}(\gamma)d\gamma$. Therefore, we get the position of the boundary agent at $\gamma$, $(x(\gamma),y(\gamma))$, as $(x(\gamma),y(\gamma))=\int_{0}^{\gamma}\frac{1}{q(\bar{\gamma})}\mathbf{s}(\bar{\gamma})d\bar{\gamma}$, and the arc-length $s(\gamma)=\int_{0}^{\gamma}\frac{1}{q(\bar{\gamma})}d\bar{\gamma}$, which is discretized by a consistent scheme to obtain:

[TABLE]

and we recall that the agents have access to $q$ and $\mathbf{s}$. The computation of $(x_{i},y_{i})$ can be implemented by propagating from the agent with $\gamma_{i}=0$ along the boundary agents in the direction as $\gamma_{i}\rightarrow 1$, along a line graph $\mathcal{G}_{\text{line}}$ (with message receipt acknowledgment). Note that this propagation can alternatively be formulated by each agent averaging appropriate variables with left and right neighbors, which will result in a process similar to a finite-time consensus algorithm.

This way, the boundary agents are localized at time $t=0$, and they update their position estimates using their velocities, for $t\geq 0$.

5.2 Pseudo-localization algorithm in two dimensions

In this subsection, we present the pseudo-localization algorithm for the agents in the interior of the spatial domain. We first describe the idea of the coordinate transformation (diffeomorphism) we employ and construct a PDE that converges asymptotically to this diffeomorphism. We then discretize the PDE to obtain the distributed pseudo-localization algorithm.

The main idea is to employ harmonic maps to construct a coordinate transformation or diffeomorphism from the spatial domain of the swarm onto the unit disk. We begin the construction with the static case, where the agents are stationary. Let $M\subseteq^{2}$ be a compact, static $2$-manifold with boundary and $N=\{(x,y)\in^{2}\,|\,(x-1)^{2}+y^{2}\leq 1\}$ be the unit disk. The manifolds $M$ and $N$ are both equipped with a Euclidean metric $g=h=\delta$.

First, we define a mapping for the boundary of $M$. Let $\Gamma:\partial M\rightarrow[0,1)$ be a parametrization of the boundary of $M$, as outlined in Section 5.1. Let $\xi:\bar{M}\rightarrow N$ be any diffeomorphism that takes the following form on the boundary of $M$:

[TABLE]

and we know that $\Gamma^{-1}[0,1)=\partial M$.

Now, from Lemma 2.8, on harmonic diffeomorphisms, there is a unique harmonic diffeomorphism, $\Psi:M\rightarrow N$, such that $\Psi=\xi$ on $\partial M$. We know that, by definition, the mapping $\Psi=(\psi_{1},\psi_{2})$ satisfies:

[TABLE]

where $\Delta$ is the Laplace operator. Let $\Psi^{*}$ be the corresponding map from the target domain $M^{*}$ to the unit disk $N$. Now, we define a function $p^{*}:N\rightarrow_{>0}$ by $p^{*}=\rho^{*}\circ({\Psi^{*}})^{-1}$, the image of the desired spatial density distribution on the unit disk, which is computed offline and is broadcasted to the agents prior to the beginning of the self-organization process. We later use $p^{*}$ to derive the distributed control law which the agents implement.

We now construct a PDE that asymptotically converges to the harmonic diffeomorphism, which we then discretize to obtain a distributed pseudo-localization algorithm. We use the heat flow equation as the basis to define the pseudo-localization algorithm, which yields a harmonic map as its asymptotically stable steady-state solution. We begin by setting up the system for a stationary swarm, for which the spatial domain is fixed.

Let $M\subset^{2}$ be a compact $2$-manifold with boundary, $N$ be the unit disk of 2, and $\mathbf{R}=(X,Y):M\rightarrow N$. The heat flow equation is given by:

[TABLE]

The heat flow equation has been studied extensively in the literature. For well-known existence and uniqueness results, we refer the reader to [13].

Lemma 5.1**.**

*(Pointwise convergence of the heat flow equation to a harmonic diffeomorphism). The solutions of the heat flow equation (23) converge pointwise to the harmonic map satisfying (22), exponentially as $t\rightarrow\infty$, from any smooth initial $\mathbf{R}_{0}\in H^{1}(M)\times H^{1}(M)$. *

Proof 5.2**.**

Let $\Psi$ be the solution to (22), which is a harmonic map by definition. Let $\tilde{\mathbf{R}}=\mathbf{R}-\Psi$ be the error where $\mathbf{R}=(X,Y)$ is the solution to (23). Subtracting (22) from (23), we obtain:

[TABLE]

*The Laplace operator $\Delta$ with the Dirichlet boundary condition in (24) is self-adjoint and has an infinite sequence of eigenvalues $0<\lambda_{1}<\lambda_{2}<\ldots$, with the corresponding eigenfunctions $\{\phi_{i}\}_{i=1}^{\infty}$ forming an orthonormal basis of $L^{2}(M)$ (where $\phi_{i}\in L^{2}(M)$ and $\Delta\phi_{i}=\lambda_{i}\phi_{i}$ for all $i$, with $\phi_{i}=0$ on the boundary) [15]. Let the initial condition be $\tilde{X}_{0}=\sum_{i=1}^{\infty}a_{i}\phi_{i}$ and $\tilde{Y}_{0}=\sum_{i=1}^{\infty}b_{i}\phi_{i}$ (where $a_{i}$ and $b_{i}$ are constants for all $i$). The solution to (24) is then given by $\tilde{X}(t,\mathbf{r})=\sum_{i=1}^{\infty}a_{i}e^{-\lambda_{i}t}\phi_{i}(\mathbf{r})$ and $\tilde{Y}(t,\mathbf{r})=\sum_{i=1}^{\infty}b_{i}e^{-\lambda_{i}t}\phi_{i}(\mathbf{r})$. Since $\lambda_{i}>0$, for all $i$, we obtain $\lim_{t\rightarrow\infty}\tilde{X}(t,\mathbf{r})=0$ and $\lim_{t\rightarrow\infty}\tilde{Y}(t,\mathbf{r})=0$, for all $\mathbf{r}\in\bar{M}$. Therefore, $\lim_{t\rightarrow\infty}\mathbf{R}(t,\mathbf{r})=\Psi(\mathbf{r})$, for all $\mathbf{r}\in\bar{M}$, and the convergence is exponential. *

We now have a PDE that converges to the diffeomorphism given by (22) for the stationary case (agents in the swarm are at rest). For the dynamic case, and to describe the algorithm while the agents are in motion, we modify (23) as follows. Let $\mathbf{R}=(X,Y):\times^{2}\rightarrow$. We are only interested in the restriction to $M(t)$, $\mathbf{R}|_{M(t)}$, at any time $t$, so we drop the restriction and just identify $\mathbf{R}\equiv\mathbf{R}_{|_{M(t)}}$. Using the relation $\frac{dX}{dt}=\partial_{t}X+\nabla X\cdot\mathbf{v}$, where $\mathbf{v}$ is a velocity field, we obtain:

[TABLE]

We now discretize (25) to derive the distributed pseudo-localization algorithm. Now, we have $\rho:\times^{2}\rightarrow_{\geq 0}$ with support $M(t)$, the density distribution of the swarm on the domain $M(t)$. We view the swarm as a discrete approximation of the domain $M(t)$ with density $\rho$, and the PDE (25) as approximated by a distributed algorithm implemented by the swarm.

Here, we propose a candidate distributed algorithm, which would yield the heat flow equation via a functional approximation. Our candidate algorithm is a time-varying weighted Laplacian-based distributed algorithm, owing to the connection between the graph Laplacian and the manifold Laplacian [4]:

[TABLE]

and a similar equation for $Y$. We show how to derive next the values for the weights $w_{ij}(t)\in$, for all $t$. First, the set of neighbors, $j\in{\mathcal{N}}_{i}(t)$, of $i$ at time $t$, are the spatial neighbors of $i$ in $M(t)$, that is, ${\mathcal{N}}_{i}(t)=\{j\in\mathcal{S}\,|\,\|\mathbf{r}_{j}(t)-\mathbf{r}_{i}(t)\|\leq\epsilon\}\equiv B_{\epsilon}(\mathbf{r}_{i}(t))$. Using $X_{i}(t+1)-X_{i}(t)=\frac{dX}{dt}\delta t$, for a small $\delta t$, we make use of a functional approximation of (26):

[TABLE]

where $d\nu=\rho~{}d\mu$ is a density-dependent measure on the manifold, and the weighting function $w$ satisfies $w(t,\mathbf{r}_{i}(t),\mathbf{r}_{j}(t))=w_{ij}(t)$, for all $i,j\in\mathcal{S}$. We note that the summation term in (26) is a special form of the integral in (27) with a Dirac measure $d\nu$ supported on the set $\{\mathbf{r}_{1}(t),\ldots,\mathbf{r}_{N}(t)\}$ at time $t$. Now, with the choice $w(t,\mathbf{r}_{i},\mathbf{s})=\frac{1}{\int_{B_{\epsilon}(\mathbf{s}(t))}\rho(t,\mathbf{\bar{s}})d\mu}$ and for very small $\epsilon$ (making $\mathcal{O}(\epsilon^{3})$ terms negligible), (27) reduces to:

[TABLE]

where $a=\frac{1}{4\epsilon}\int_{B_{\epsilon}(\mathbf{r}_{i}(t))}(\mathbf{s}-\mathbf{r}_{i}(t))\cdot(\mathbf{s}-\mathbf{r}_{i}(t))d\mu$ is a constant. Now, with the choice $\delta t=a$, we obtain:

[TABLE]

which is the PDE (25). Let $d(t,\mathbf{r}_{i}(t))=\int_{B_{\epsilon}(\mathbf{r}_{i}(t))}\rho(t,\mathbf{s})d\mu$ and $d_{i}(t)=|\mathcal{N}_{i}(t)|$, for $i\in\mathcal{S}$. Substituting $w_{ij}(t)=w(t,\mathbf{r}_{i}(t),\mathbf{r}_{j}(t))=\frac{1}{\int_{B_{\epsilon}(\mathbf{r}_{j}(t))}\rho(t,\mathbf{\bar{s}})d\mu}=\frac{1}{d(t,\mathbf{r}_{j}(t))}\approx\frac{1}{d_{j}(t)}$, in (26), we get the distributed pseudo-localization algorithm for the agents in the interior of the swarm to be:

[TABLE]

For the agents on the boundary $\partial M(t)$, we have:

[TABLE]

where $\xi_{i}=\xi(\mathbf{r}_{i}(t))$, for $\mathbf{r}_{i}(t)\in\partial M(t)$. Note that the discretization scheme is consistent, in that as the number of agents $N\rightarrow\infty$, the discrete equation (28) converges to the PDE (25). In this way, from (28), the pseudo-localization algorithm is a Laplacian-based distributed algorithm, with a time-varying weighted graph Laplacian.

5.3 Distributed density control law and analysis

In this section, we derive the distributed feedback control law to converge to the desired density distribution over the target domain in the two-dimensional case. The swarm dynamics are given by:

[TABLE]

{assumption}

(Well-posedness of the PDE system). We assume that (29) is well-posed, and that its solution $\rho(t,\cdot)$ is sufficiently smooth and is bounded in the Sobolev space $H^{1}(\cup_{t}M(t))$, the components of the velocity field $\mathbf{v}$ are bounded in the Sobolev space $H^{1}(\cup_{t}M(t))$ and of the parametrized velocity on the boundary are bounded in the Sobolev space $H^{1}((0,1))$.

In what follows, we describe the control strategy based on three different stages.

5.3.1 Stage $1$

In this stage, the objective is for the swarm to converge to the target spatial domain $M^{*}$.

Let $\mathbf{r}^{*}:[0,1]\rightarrow\partial M^{*}$ be the closed curve describing the desired boundary. Let $\mathbf{e}(\gamma)=\mathbf{r}(\gamma)-\mathbf{r}^{*}(\gamma)$ be the position error of agent $\gamma$ on the boundary, where $\mathbf{r}(\gamma)$ is the actual position of agent $\gamma$ computed as presented in Section 5.1. We define a distributed control law for swarm motion as follows:

[TABLE]

Theorem 5.3**.**

*(Convergence to the desired spatial domain).

Under the well-posedness Assumption 5.3, the domain $M(t)$ of the system (29), with the distributed control law (30) converges to the target spatial domain $M^{*}$ as $t\rightarrow\infty$, from any initial domain $M_{0}$ with smooth boundary. *

Proof 5.4**.**

We consider an energy functional $E$ given by:

[TABLE]

Its time derivative, $\dot{E}$, using (30), is given by:

[TABLE]

Clearly, $\dot{E}\leq 0$, and considering a parametrization of $\partial M(t)$ by the interval $[0,1)$, we have  $\mathbf{v}(t,\cdot)\in H^{1}((0,1))$ and bounded. By Lemma 2.6, the Rellich-Kondrachov Compactness theorem, $H^{1}((0,1))$ is compactly contained in $L^{2}((0,1))$ (and we also have that $H^{1}((0,1))$ is dense in $L^{2}((0,1))$). Thus, by the LaSalle Invariance Principle, Lemma 2.7, we have that the solutions to (29) with the control law (30) converge in the $L^{2}$-norm to the largest invariant subset of $\dot{E}^{-1}(0)$, which satisfies:

[TABLE]

The set $\dot{E}^{-1}(0)$ is characterized by the first equality above and the second equality is further satisfied by the invariant subset of $\dot{E}^{-1}(0)$. We know from (30) that $\partial_{t}\mathbf{v}=-\mathbf{e}-\mathbf{v}$ on $\partial M(t)$, which upon multiplying on both sides by $\mathbf{v}$, integrating over $\partial M(t)$ and applying the previous equality on the integral of $\mathbf{v}\cdot\partial_{t}\mathbf{v}$, yields $\lim_{t\rightarrow\infty}\int_{\partial M(t)}\mathbf{e}\cdot\mathbf{v}=0$. Now, we have $|\partial_{t}\mathbf{v}|^{2}=|\mathbf{e}|^{2}+|\mathbf{v}|^{2}+2\mathbf{e}\cdot\mathbf{v}$, which on integrating over $\partial M(t)$ yields $\lim_{t\rightarrow\infty}\||\partial_{t}\mathbf{\mathbf{v}}|\|_{L^{2}(\partial M(t))}=\lim_{t\rightarrow\infty}\||\mathbf{e}|\|_{L^{2}(\partial M(t))}$. By multiplying $\partial_{t}\mathbf{v}=-\mathbf{e}-\mathbf{v}$ on both sides by $\partial_{t}\mathbf{v}$, integrating over $\partial M(t)$, and using the Cauchy-Schwarz inequality, we obtain:

[TABLE]

*In this way, the Cauchy-Schwarz inequality becomes an equality, which implies that $\lim_{t\rightarrow\infty}\int_{\partial M(t)}\left[|\mathbf{e}||\partial_{t}\mathbf{v}|-(-\mathbf{e})\cdot\partial_{t}\mathbf{v}\right]=0$ (since the integrand is non-negative and its integral is zero, it is zero almost everywhere), thus $\lim_{t\rightarrow\infty}\partial_{t}\mathbf{v}=-\lim_{t\rightarrow\infty}\mathbf{e}$ almost everywhere (a.e.) on the boundary, and, in turn, implies that $\lim_{t\rightarrow\infty}\mathbf{v}=0$ a.e. on the boundary (since $\partial_{t}\mathbf{v}=-\mathbf{e}-\mathbf{v}$ and $\lim_{t\rightarrow\infty}\partial_{t}\mathbf{v}=-\lim_{t\rightarrow\infty}\mathbf{e}$). From here, and owing to the Invariance Principle, we have $\lim_{t\rightarrow\infty}\partial_{t}\mathbf{v}=0=\lim_{t\rightarrow\infty}\mathbf{e}$ a.e. on the boundary. Thus, we have that $\lim_{t\rightarrow\infty}M(t)=M^{*}$. *

5.3.2 Stage $2$

Here, the agents in the swarm implement the pseudo-localization algorithm presented in Section 5.2. Since the agents are distributed across the target spatial domain $M^{*}$, implementing the pseudo-localization algorithm yields the coordinate transformation $\Psi^{*}$ characteristic of the domain $M^{*}$. We therefore have $\partial_{t}\Psi^{*}=0$, which implies that $\frac{d\Psi^{*}}{dt}=\partial_{t}\Psi^{*}+\nabla(\Psi^{*})\mathbf{v}=\nabla(\Psi^{*})\mathbf{v}$, which will be used in Stage $3$.

5.3.3 Stage $3$

In this stage, the boundary agents of the swarm remain stationary and interior agents converge to the desired density distribution.

Consider the distributed control law, defined as follows for all time $t$:

[TABLE]

where $\frac{d\mathbf{v}}{dt}$ at $\mathbf{r}\in M$ is the acceleration of the agent at $\mathbf{r}$, the control input. Using the relation $\frac{d}{dt}=\partial_{t}+\mathbf{v}\cdot\nabla$, it follows from (31) that $\partial_{t}\mathbf{v}=-\rho\nabla(\rho-p^{*}\circ\Psi^{*})+\Delta\mathbf{v}-\mathbf{v}$.

Theorem 5.5**.**

*(Convergence to the desired density).

The solutions $\rho(t,\cdot)$ to (29) for the fixed domain $M^{*}$, under the distributed control law (31) and the well-posedness Assumption 5.3, converge to the desired density distribution $\rho^{*}$ in the $L^{2}$-norm as $t\rightarrow\infty$.*

Proof 5.6**.**

We consider an energy functional $E$ given by:

[TABLE]

Using Corollary 2.3, to compute the derivative of energy functionals, we obtain $\dot{E}$ (letting $\bar{\nabla}=\left(\partial_{X},\partial_{Y}\right)$) as follows:

[TABLE]

where, to obtain the third equality, we expand the square $|\rho-p^{*}\circ\Psi^{*}|^{2}$ in the second integral of the second equality. Since $\mathbf{v}=0$ on $\partial M^{*}$ and from Section 5.3.2, we have $\frac{d\Psi^{*}}{dt}=\nabla(\Psi^{*})\mathbf{v}$, we obtain:

[TABLE]

We have $\bar{\nabla}p^{*}\nabla\Psi^{*}=\nabla(p^{*}\circ\Psi^{*})$, and $\nabla(\rho^{2}-(p^{*}\circ\Psi^{*})^{2})=(\rho-p^{*}\circ\Psi^{*})\nabla(\rho+p^{*}\circ\Psi^{*})+(\rho+p^{*}\circ\Psi^{*})\nabla(\rho-p^{*}\circ\Psi^{*})$. Thus, we get:

[TABLE]

*We therefore get: *

[TABLE]

From (31), we have $\partial_{t}\mathbf{v}=-\rho\nabla(\rho-p^{*}\circ\Psi^{*})+\Delta\mathbf{v}-\mathbf{v}$, and we obtain:

[TABLE]

Clearly, $\dot{E}\leq 0$, with $\rho(t,.),\mathbf{v}\in H^{1}(M^{*})$ and bounded (by Assumption 5.3). By Lemma 2.6, the Rellich-Kondrachov Compactness theorem, $H^{1}(M^{*})$ is compactly contained in $L^{2}(M^{*})$ (and we also know that the set of all $(\rho,\mathbf{v})$ satisfying Assumption 5.3 is dense in $L^{2}(M^{*})$). Thus, by the Invariance Principle, Lemma 2.7, we have that the solution to (29) converges in the $L^{2}$-norm to the largest invariant subset of $\dot{E}^{-1}(0)$, which satisfies:

[TABLE]

The set $\dot{E}^{-1}(0)$ is characterized by the first equality above and the second equality is further satisfied by the invariant subset of $\dot{E}^{-1}(0)$. We know from (31) that

[TABLE]

*which substituted in (32) yields $\int_{M^{*}}\rho\mathbf{v}\cdot\nabla(\rho-p^{*}\circ\Psi^{*})=0$. Now, from (33), we obtain $\||\partial_{t}\mathbf{v}|\|_{L^{2}(M^{*})}^{2}=\int_{M^{*}}|\rho\nabla(\rho-p^{*}\circ\Psi^{*})|^{2}+\int_{M^{*}}|\mathbf{v}|^{2}+2\int_{M^{*}}\rho\mathbf{v}\cdot\nabla(\rho-p^{*}\circ\Psi^{*})=\int_{M^{*}}|\rho\nabla(\rho-p^{*}\circ\Psi^{*})|^{2}$; that is, $\||\partial_{t}\mathbf{v}|\|_{L^{2}(M^{*})}=\||\rho\nabla(\rho-p^{*}\circ\Psi^{*})|\|_{L^{2}(M^{*})}$. By multiplying (33) by $\partial_{t}\mathbf{v}$ on both sides and applying the Cauchy-Schwarz inequality, we can also get that $\||\partial_{t}\mathbf{v}|\|_{L^{2}(M^{*})}^{2}=-\int_{M^{*}}\rho\partial_{t}\mathbf{v}\cdot\nabla(\rho-p^{*}\circ\Psi^{*})\leq\int_{M^{*}}|\partial_{t}\mathbf{v}||\rho\nabla(\rho-p^{*}\circ\Psi^{*})|\leq\||\partial_{t}\mathbf{v}|\|_{L^{2}(M^{*})}||\rho\nabla(\rho-p^{*}\circ\Psi^{*})|\|_{L^{2}(M^{*})}=\||\partial_{t}\mathbf{v}|\|_{L^{2}(M^{*})}^{2}$. Thus, the Cauchy-Schwarz inequality is in fact an equality, which implies that $\partial_{t}\mathbf{v}=-\rho\nabla(\rho-p^{*}\circ\Psi^{*})$ almost everywhere in $M^{*}$, which, from (33) implies in turn that $\mathbf{v}=0$ a.e. in $M^{*}$. It thus follows that $\partial_{t}\mathbf{v}=0$ and $\nabla(\rho-p^{*}\circ\Psi^{*})=0$ a.e in $M^{*}$, and therefore $\rho-p^{*}\circ\Psi^{*}$ is constant a.e. in $M^{*}$. Using the Poincare-Wirtinger inequality, Lemma 2.5, we obtain that $\|(\rho-p^{*}\circ\Psi^{*})-(\rho-p^{*}\circ\Psi^{*})_{M^{*}}\|\leq C\|\nabla(\rho-p^{*}\circ\Psi^{*})\|=0$, where $(\rho-p^{*}\circ\Psi^{*})_{M^{*}}=\frac{1}{|M^{*}|}\int_{M^{*}}(\rho-p^{*}\circ\Psi^{*})$. Since $\int_{M^{*}}\rho=\int_{N}p^{*}=\int_{M^{*}}p^{*}\circ\Psi^{*}=1$, we have that $(\rho-p^{*}\circ\Psi^{*})_{M^{*}}=0$, and therefore $\|\rho-p^{*}\circ\Psi^{*}\|_{L^{2}(M^{*})}=0$. *

5.3.4 Robustness of the distributed control law

The self-organization algorithm in $2$D has been divided into three stages, where asymptotic convergence is achieved in each stage (with exponential convergence in the second stage). We now present a robustness result for convergence in Stage $3$ under incomplete convergence in the preceding stages.

Lemma 5.7**.**

*(Robustness of the control law). For every $\delta>0$, there exist $T_{1},T_{2}<\infty$ such that when Stages $1$ and $2$ are terminated at $t_{1}>T_{1}$ and $t_{2}>T_{2}$ respectively, we have that $\lim_{t\rightarrow\infty}\|\rho(t,\cdot)-\rho^{*}\|_{L^{2}(M(t_{1}))}<\delta$. *

Proof 5.8**.**

*In Stage $1$, it follows from Theorem 5.3 on the convergence to the desired spatial domain that $\lim_{t\rightarrow\infty}M(t)=M^{*}$. Then for every $\epsilon_{1}>0$, we have $T_{1}<\infty$, such that $d_{H}(M(t),M^{*})<\epsilon_{1}$ for all $t>T_{1}$, where $d_{H}$ is the Hausdorff distance between two sets; see (1). (Note that any appropriate notion of distance can alternatively be used here.) Let Stage $1$ be terminated at $t_{1}>T_{1}$, which implies that the swarm is distributed across the domain $M(t_{1})$. In Stage $2$, it follows from Lemma 5.1 on the convergence of the heat flow equation to the harmonic map, that for a domain $M(t_{1})$, we have that $\lim_{t\rightarrow\infty}\mathbf{R}(t,\cdot)=\Psi_{M(t_{1})}$ pointwise, where $\Psi_{M(t_{1})}$ is the harmonic map from $M(t_{1})$ to $N$ (the unit disk). Then, for every $\epsilon_{2}>0$, we have a $T_{2}<\infty$, such that $\|\mathbf{R}(t,\cdot)-\Psi_{M(t_{1})}\|_{\infty}<\epsilon_{2}$ for all $t>T_{2}$. Let Stage $2$ be terminated at $t_{2}>T_{2}$, which implies that the map from the spatial domain to the disk is $\mathbf{R}(t_{2},\cdot)$. In Stage $3$, it follows from the arguments in the proof of Theorem 5.5 (on the convergence to the desired density distribution) that $\lim_{t\rightarrow\infty}\rho(t,\cdot)=p^{*}\circ\mathbf{R}(t_{2},\cdot)$ a.e. in $M(t_{1})$ if the map at the end of Stage $2$ is $\mathbf{R}(t_{2},\cdot)$. We characterize the error as $\lim_{t\rightarrow\infty}\|\rho-\rho^{*}\|_{L^{2}(M(t_{1}))}=\|p^{*}\circ\mathbf{R}(t_{2},\cdot)-p^{*}\circ\Psi^{*}\|_{L^{2}(M(t_{1}))}=\|p^{*}\circ\mathbf{R}(t_{2},\cdot)-p^{*}\circ\Psi_{M(t_{1})}+p^{*}\circ\Psi_{M(t_{1})}-p^{*}\circ\Psi^{*}\|_{L^{2}(M(t_{1}))}\leq\|p^{*}\circ\mathbf{R}(t_{2},\cdot)-p^{*}\circ\Psi_{M(t_{1})}\|_{L^{2}(M(t_{1}))}+\|p^{*}\circ\Psi_{M(t_{1})}-p^{*}\circ\Psi^{*}\|_{L^{2}(M(t_{1}))}$. Recall that $\|\mathbf{R}(t_{2},\cdot)-\Psi_{M(t_{1})}\|_{\infty}<\epsilon_{2}$, and since $p^{*}$ is Lipschitz, we can get the bound $\|p^{*}\circ\mathbf{R}(t_{2})-p^{*}\circ\Psi_{M(t_{1})}\|_{L^{2}(M(t_{1}))}<\delta_{1}=c\epsilon_{2}$ (where $c$ is the Lipschitz constant times the area of $M(t_{1})$). The harmonic map also depends continuously on its domain [19], which yields the bound $\|\Psi_{M(t_{1})}-\Psi^{*}\|_{\infty}<\epsilon_{3}$, since $d_{H}(M(t_{1}),M^{*})<\epsilon_{1}$. Thus, we get another bound $\|p^{*}\circ\Psi_{M(t_{1})}-p^{*}\circ\Psi^{*}\|_{L^{2}(M(t_{1}))}<\delta_{2}=c\epsilon_{3}$, and that $\|\rho-\rho^{*}\|_{L^{2}(M(t_{1}))}<\delta_{1}+\delta_{2}=\delta$. Therefore, going backwards, for all $\delta>0$, we can find $T_{1}$ and $T_{2}$ such that the density error is bounded by $\delta$, when the Stages $1$ and $2$ are terminated at $t_{1}>T_{1}$ and $t_{2}>T_{2}$ respectively. *

5.4 Discrete implementation

In this section, we present consistent schemes for discrete implementation of the distributed control laws (30) and (33), where the key aspect is the computation of spatial gradients (of $\rho$ in Stage $1$, and of $\rho$, $\Psi^{*}$ and the components of velocity $\mathbf{v}$ in Stage $3$). The network graph underlying the swarm is a random geometric graph, where the nodes are distributed according to the density distribution over the spatial domain. According to this, every agent communicates with other agents within a disk of given radius (say $r$) determined by the hardware capabilities, which reduces to the graph having an edge between two nodes if and only if the nodes are separated by a distance less than $r$. We recall the earlier stated assumption that the agents know the true $x$- and $y$-directions.

5.4.1 On the computation of $p^{*}$

We first begin with an approach to compute offline the map $p^{*}$ via interpolation. Let the desired domain $M^{*}\in^{2}$ be discretized into a uniform grid to obtain $M^{*}_{d}=\{\mathbf{r}_{1},\ldots,\mathbf{r}_{m}\}$ (the centers of finite elements, where $\mathbf{r}_{k}=(x_{k},y_{k})$). The desired density $\rho^{*}:M^{*}\rightarrow{\mathbb{R}}_{>0}$ is known, and we compute the value of $\rho^{*}$ on $M^{*}_{d}$ to get $\rho^{*}(\mathbf{r}_{1},\ldots,\mathbf{r}_{m})=(\rho^{*}_{1},\ldots,\rho^{*}_{m})$. We also have $\Psi^{*}(x,y)=(X^{*},Y^{*})\in N$, for all $(x,y)\in M^{*}$. Now, computing the integral with respect to the Dirac measure for the set $M^{*}_{d}$, we obtain $\Psi^{*}(\mathbf{r}_{1},\ldots,\mathbf{r}_{m})=(\Psi^{*}_{1},\ldots,\Psi^{*}_{m})$. The value of the function $p^{*}$ at any $(X,Y)\in N$ can be obtained from the relation $p^{*}(\Psi^{*}_{1},\ldots,\Psi^{*}_{m})=\rho^{*}(\mathbf{r}_{1},\ldots,\mathbf{r}_{m})$ for $k=1,\ldots,m$ by an appropriate interpolation.

5.4.2 Discrete control law

As stated earlier, for the discrete implementation of the distributed control laws (30) and (33), the key aspect is the computation of spatial gradients (of $\rho$ in Stage $1$, and of $\rho$, $\Psi^{*}$ and the components of velocity $\mathbf{v}$ in Stage $3$). In the subsequent sections we present two alternative, consistent schemes for computing the spatial gradient (of any smooth function, with the above being the ones of interest), one using the Jacobian of the harmonic map and the other without it.

Computing the Jacobian of the harmonic map

Let $J(\mathbf{r})=\nabla\Psi(\mathbf{r})$ be the (non-singular) Jacobian of the harmonic diffeomorphism $\Psi:M\rightarrow N$. When the steady-state is reached in the pseudo-localization algorithm (28) (i.e., $X_{i}(t+1)=X_{i}(t)=\psi_{1}^{i}$ and $Y_{i}(t+1)=Y_{i}(t)=\psi_{2}^{i}$), we have, $\forall\,i\in\mathcal{S}$:

[TABLE]

where $i$ is the index of the agent located at $\mathbf{r}\in M$ and $\mathcal{N}_{i}$ is the set of agents in a disk-shaped neighborhood $B_{\epsilon}(\mathbf{r})$ of area $\epsilon$ centered at $\mathbf{r}$. Rewriting the above, we get, $\forall\,i\in\mathcal{S}$:

[TABLE]

We assume that the agents have the capability in their hardware to perturb the disk of communication $B_{\epsilon}(\mathbf{r})$ (by moving an antenna, for instance). The Jacobian $J=\nabla\Psi$ , where $\Psi=(\psi_{1},\psi_{2})$ is computed through perturbations to $\mathcal{N}_{i}$ (i.e., the neighborhood $B_{\epsilon}(\mathbf{r})$) and using consistent discrete approximations:

[TABLE]

and similarly for $\psi_{2}$. Now, $\psi_{1}(\mathbf{r}+\delta x\mathbf{e}_{1})$ is computed as in (34) for $\mathcal{N}^{\delta x}_{i}$, the set of agents in $B_{\epsilon}(\mathbf{r}+\delta x\mathbf{e}_{1})$ and $\psi_{1}(\mathbf{r}+\delta y\mathbf{e}_{2})$ from $B_{\epsilon}(\mathbf{r}+\delta y\mathbf{e}_{2})$.

Computing the spatial gradient of a smooth

function using the Jacobian of $\Psi$

Let $\nabla=\left(\partial_{x},\partial_{y}\right)$ and $\bar{\nabla}=\left(\partial_{\psi_{1}},\partial_{\psi_{2}}\right)$, where $\Psi=(\psi_{1},\psi_{2})$. We have $\partial_{x}=(\partial_{x}\psi_{1})\partial_{\psi_{1}}+(\partial_{x}\psi_{2})\partial_{\psi_{2}}$ and $\partial_{y}=(\partial_{y}\psi_{1})\partial_{\psi_{1}}+(\partial_{y}\psi_{2})\partial_{\psi_{2}}$. Therefore, $\nabla=J^{\top}\bar{\nabla}$. For a smooth function $f:M\rightarrow$, we have, $\nabla f=J^{\top}\bar{\nabla}f$, and the agents can numerically compute $\bar{\nabla}$ by:

[TABLE]

where $i$ is the index of the agent located at $\mathbf{r}\in M$ and $\mathcal{N}_{i}$ is the set of agents in a ball $B_{\epsilon}(\mathbf{r})$.

Computing the spatial gradient of a smooth function

without the Jacobian of $\Psi$

In the absence of a Jacobian estimate, we use the following alternative method for computing an approximate spatial gradient estimate of a smooth function. This is used in Stage $1$ of the self-organization process.

Let $\bar{f}(\mathbf{r})$ be the mean value of $f$ over a ball $B_{\epsilon}(\mathbf{r})$:

[TABLE]

We have:

[TABLE]

Similarly,

[TABLE]

In all, for any scalar function $f$, each agent can use the approximation:

[TABLE]

to estimate of the gradient $\nabla f$.

5.4.3 On the convergence of the discrete system

We have noted earlier that the pseudo-localization algorithm (28) satisfies the consistency condition in that as $N\rightarrow\infty$, Equation (28) converges to the PDE (25). The pseudo-localization algorithm is also essentially a weighted Laplacian-based distributed algorithm that is stable. Thus, by the Lax Equivalence theorem [33], the solution of (28) converges to the solution of (25) as $N\rightarrow\infty$. However, for the distributed control laws in Stages $1$-$3$, we are only able to provide consistent discretization schemes. The dynamics of the swarm (29) with the control laws (30) and (31) are nonlinear for which is no equivalent convergence theorem. Further analysis to determine convergence is required, which falls out the scope of this present work.

6 Numerical simulations

In this section, we present numerical simulations of swarm self-organization, that is, of the control laws presented in Sections 4.2 and of Section 5.3.

6.1 Self-organization in one dimension

In the simulation of the $1$D case, we consider a swarm of $N=10000$ agents, the desired density distribution is given by $\rho^{*}(x)=a\sin(x)+b$, where $a=1-\frac{\pi}{2N}$ and $b=\frac{1}{N}$, $x\in\left[0,\frac{\pi}{2}\right]$. We use a kernel-based method to approximate the continuous density function, which is given by:

[TABLE]

is a flat kernel and $c_{d}\in{\mathbb{R}}_{>0}$ is a constant [9]. We discretize the spatial domain with $\Delta x=0.001$ units, and use an adaptive time step. The self-organization begins from an arbitrary initial density distribution. Figure 2 shows the initial density distribution, an intermediate distribution and the final distribution. We observe that there is convergence to the desired density distribution, even with noisy density measurements.

6.2 Self-organization in two dimensions

In the simulation of the $2$D case, we first present in Figure 3 the evolution of the boundary of the swarm in Stage $1$, where the swarm converges to the target spatial domain $M^{*}$ from an initial spatial domain. The target spatial domain, a circle of radius $0.5$ units, given by $M^{*}=\{(x,y)\in^{2}\,|\,(x-0.6)^{2}+y^{2}\leq 0.25\}$, with the desired density distribution $\rho^{*}$ given by $\rho^{*}(x,y)=\frac{1}{\left((x-0.4)^{2}+y^{2}\right)^{0.3}}$.

We present in Figures 4 and 5 the result of implementation of the pseudo-localization algorithm with the steady state distributions of $\Psi^{*}=(\psi^{*}_{1},\psi^{*}_{2})$ respectively. We note that the steady state distribution $\Psi^{*}$ as a function of the spatial coordinates $(x,y)$ in this case is linear.

Next, we focus on Stage $3$ of the self-organization process, where the agents already distributed over the target spatial domain, converge to the desired density distribution. The initial density distribution of the swarm is uniform, and the distributed control law of Stage $3$ in Section 5.3 is implemented. Figure 7 shows the density distribution at a few intermediate time instants of implementation and figure 7 shows the spatial density error plot, where $e(\rho)=\int_{M^{*}}|\rho-\rho^{*}|^{2}$ is the spatial density error. The results show convergence as desired.

7 Conclusions

In this paper, we considered the problem of self-organization in multi-agent swarms, in one and two dimensions, respectively. The primary contribution of this paper is the analysis and design of position and index-free distributed control laws for swarm self-organization for a large class of configurations. This was accomplished through the introduction of a distributed pseudo-localization algorithm that the agents implement to find their position identifiers, which then use in their control laws. The validation of the results for more general non-simply connected domains will be considered in the future. An extension to this work will involve the characterization of constraints on the local density function to capture finite robot sizes and collision avoidance constraints, as well as accounting for possible non-holonomic constraints on the motion of the robots.

Acknowledgments

The authors would like to thank Prof. Lei Ni at the UC San Diego Mathematics Department and the reviewers of this manuscript for their valuable inputs.