Collective motion planning for a group of robots using intermittent diffusion

TL;DR

This paper introduces a novel robot group motion planning method based on optimal transport theory, enabling complex formations with proven collision avoidance and convergence guarantees.

Contribution

It presents a new approach leveraging optimal transport for collective robot motion planning, with rigorous proofs of safety and convergence.

Findings

Effective shape formation and assembly achieved

Collision avoidance is guaranteed

Convergence of the algorithm is rigorously proven

Abstract

In this work we establish a simple yet effective strategy, based on optimal transport theory, for enabling a group of robots to accomplish complex tasks, such as shape formation and assembly. We demonstrate the feasibility of this approach and rigorously prove collision avoidance and convergence properties of the proposed algorithms.

Tables

Table 1. Table 1: Simulation parameters

Symbol	Description	Value
$G_0$	Repelling function amplitude	.01
$R$	Robot sensor radius	$10r$
$\Delta t$	Time step	$.1r$
$\alpha$	ID Diffusion scale	$r$
$\beta$	ID Time scale	10
$M$	Computational domain size	6

Table 2. Table 2: Final objective function value for both target shapes and varied robot radii ( r 𝑟 r ), starting from a random initialization.

$r$	$N$	$\Psi \left(X_{ID}\right)$	$\Psi \left(X_{GD}\right)$
‘Q’
.1	$50$	$0.02537$	$0.02808$
.05	$150$	$0.02905$	$0.02931$
.01	$1,000$	$0.00066$	$0.00223$
‘Jie’
.1	$200$	$0.12891$	$0.13273$
.05	$400$	$0.05964$	$0.06545$
.01	$3,000$	$0.00702$	$0.01648$

Equations

$$dY(t)=-\nabla \Psi (Y(t))dt+\sigma dW(t),$$

$$\frac{\partial \rho }{\partial t}=\nabla \cdot (\nabla \Psi (y)\rho )+\frac{1}{2}{\sigma }^2\Delta \rho .$$

$${\rho }^{\ast }(y)$$

$$\text{where}P$$

$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(W_{2}(\rho^{1},\rho^{2}))^{2}}=\inf_{v}{\int_{0}^{1}\int_{\Omega}v^{2}\rho dydt},$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (v\rho )=0,$$

$$\rho (0,y)={\rho }^1(y),\rho (1,y)={\rho }^2(y),$$

$$\mathcal{G}(\rho )={\int }_{\Omega }\Psi (y)\rho (y)dy+\frac{1}{2}{\sigma }^2{\int }_{\Omega }\rho (y)\log \rho (y)dy,$$

$$X(t)=(X^1(t),\dots ,X^N(t)),$$

$$F(X)=\frac{1}{N}\sum _{i=1}^N\mu (X^i);\mu (X^i)=\underset{X^{\prime }\in \Gamma }{\min }\Vert X^i-X^{\prime }{\Vert }^2,$$

$$G(X)=\{\begin{array}{ll}{G}_0\sum _{i=1}^N\sum _{j\ne i}\phi (\Vert X^i-X^j\Vert /2),& \text{if }\Vert X^i-X^j\Vert <R,\\ 0,& \text{ otherwise}\end{array},$$

$$\underset{x\to R^{-}}{\lim }\phi (x)=\underset{x\to R^{-}}{\lim }{\phi }^{\prime }(x)=0.$$

$$\Psi (X)=F(X)+G(X).$$

$$\frac{d{X}^i(t)}{dt}=-(\nabla \Psi (X(t)){)}_i.$$

$$\{\begin{array}{ll}d{Y}^i(t)=-(\nabla \Psi (Y(t)){)}_{i}dt+\sigma (t)dW(t),& t>0,\\ Y(0)=Y_0,& \end{array}$$

$$\sigma (t)=\{\begin{array}{ll}0& \text{if t∈[Sk,Tk]}\\ {\sigma }_k& \text{if t∈[Tk−1,Sk]}.\end{array}$$

$$\hat{F}(X)=\frac{1}{N}\sum _{i=1}^N\Vert X^i-{\hat{Y}}^i{\Vert }^2.$$

$$\frac{d{X}^i(t)}{dt}=-{\left(\nabla \left(\hat{F}(X(t))+G(X(t))\right)\right)}_i.$$

$$X_{n+1}^i=X_n^i-\Delta t(\nabla \Psi (X_n)){)}_i,$$

$$Y_{n+1}^i=Y_n^i-\Delta t(\nabla \Psi (Y_n){)}_i+{\xi }_n\sqrt{\Delta t},$$

$$\displaystyle{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}Y_{m+1}^{i}=Y_{m}^{i}-(\frac{1}{N}\nabla\mu(Y_{m}^{i})+(\nabla G(Y_{m}))_{i})\Delta t+\sigma\xi^{i}_{m}\sqrt{\Delta t},}$$

$${\hat{F}}_i(X^i)$$

$$X_{n+1}^i$$

$$G(X)=G_0\sum _{i=1}^N\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^N\cot (\pi /2(\Vert X^i-X^j{\Vert }^2)/R^2).$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathcal{X}=\{(X^{1},\ldots,X^{N})\mid X^{i}\in\Omega,\inf_{i,j\neq i}\|X^{{i}}-X^{j}\|>r>0\}.}$$

$$m_0:=\underset{i\ne j}{\min }\{\Vert X_0^i-X_0^j\Vert \}>r.$$

$$G_0\phi (r)>E_0:=\Psi (X_0).$$

$$\phi (r)>N^2\phi (m_0)+\frac{2M}{G_0}.$$

$$\underset{i,j\ne i}{\inf }\Vert X^i(t)-X^j(t){\Vert }^2>r^2$$

$$\frac{d}{dt}\left(\Psi (X)\right)$$

Full text

{AMS}

68Q25, 68R10, 68U05

1 Introduction

Motion planning for multi-robot systems has drawn significant attention in recent years due to the emergence of a number of new application scenarios, e.g., [40, 69]. Compared to single robot systems, multi-robot systems have many benefits, including spatial distribution, efficiency and robustness at completing a task due to division of labor, localization, information-sharing, redundancy, and potentially lower cost. On the other hand, motion planning for multi-robot systems must face significant challenges, such as collisions, deadlock due to the presence of local minima in the multi-objective functions from which the controllers are derived, and uncertainty introduced from the environment and stochastic effects in the system [38]. Computationally, the motion planning problem can be NP-hard and not solvable in polynomial time even for some two -dimensional cases [60]. Furthermore, all of these difficulties are exacerbated when the robots are limited in capabilities, for example, short -range communications. Addressing those existing challenges and satisfying the ever -growing desire for new missions demand novel strategies and developments in both control engineering and their underlying mathematical theory.

There is a vast literature for path planning that spans widely known methods, including graph -based approaches such as A*, D*, or D*-lite, [13, 18, 24, 37, 35, 36, 65, 45], randomized algorithms such as Probabilistic Road Maps (PRM) [55, 33, 2, 56], and tree-search algorithms, including Rapidly-exploring Random Tree (RRT) [52, 41, 20, 57, 32]. These methods find trajectories, often optimal ones, by generating feasible paths defined by nodes on a lattice or random tree that characterizes the space of possible configurations.

Much progress has been made in adapting existing methods to cooperative path planning problems for relatively small groups of robots [4, 7, 19, 27, 42, 54, 60, 23, 61, 48, 1, 69, 25, 49] or the design of cooperative motion strategies without explicit preplanning of optimal paths [64, 47, 14, 42]. Readers are referred to a few survey papers on aerial swarm robots [10] and collective behavior of multi-agent algorithms [63] that provide extensive lists of papers and summaries of many methods appeared in recent years. It is worth noting that one of the conclusions in [63] highlights the artificial potential functions (APF) method for its versatility, simplicity, scalability and high expressivity in swarm robots, and calls for new developments in both theory and algorithms that share the key properties of APF.

APF is proposed in [34]. It formulates the shape-formation problem as a problem of minimizing a potential composed of an attractive field, based on the desired shape, and a repelling field based on obstacles. Designed originally for single-robot trajectories [8, 67], these theories and methods have been extended and improved upon over the past several decades, including the addition of simulated annealing and an extension to dynamic environments [62, 22, 58, 68, 59, 51]. Due to its simplicity and scalability, APF methods can handle large groups of robots, in which each robot regards others as obstacles, and higher dimensional problems efficiently. Recently, in [26], potential based methods were succesfully used to develop decentralized controllers for shape formation of a swarm of robots. However, a well-known limitation of APF is the presence of local minima caused by the repelling forces of obstacles, leading to potential deadlocks.

In this paper, we advocate designing motion planning methods for multi-robot systems by equipping APF with new ideas, such as intermittent diffusion, in recent developments in stochastic differential equations (SDEs) and global optimization, and Wasserstein gradient flows in probability space. We cast the motion planning for a group of robots as transporting one point-mass distribution (initial shape) to another point-mass distribution (target shape). Unlike many existing motion planning problems in which each robot knows its target configuration, we do not assume that the robots know their precise destination, rather they must form the desired shape or distribution in the end. We propose a strategy that produces algorithms to control the group dynamics using carefully designed potentials and stochasticity. Our contributions include

Design two dynamical systems, based on the idea of intermittent diffusion, that alternately produce the motion trajectories for a group of robots. 2. 2.

Prove that our strategy produce s collision-free motions in both continuous and discrete settings. 3. 3.

Prove the convergence to the desired shape by using optimal transport theory. Demonstrating the approach can overcome the problem of local minima and deadlocks.

It is worth mentioning that our approach is closely related to the theory of optimal transport [31, 66], a mathematical branch that finds many successful applications in optics, econometrics, and computer graphics [3, 12, 17, 21, 50, 70], just to name a few. The connection between our approach and optimal transport theory has two different aspects. In theory, our proof for the convergence is through intermittent diffusion, whose proof relies on optimal transport theory. In our algorithm, the paths produced can be viewed as randomly perturbed particle motions, whose distribution density satisfies the well-known Fokker-Planck equation, which is regarded as a gradient flow of the relative entropy in optimal transport theory. This gradient flow viewpoint, which ensures its convergence to the Gibbs distribution, inspired our design. We use the target shape to create a Gibbs distribution, which guides the particle motions to produce the path.

We also want to mention that the proposed method differs from similar applications of optimal transport to robot path-planning [5, 11, 39], in which either linear programming, quadratic programming, or primal-dual method is used to identify the transport map. Instead of resorting to optimization methods in computations, we directly prescribe the gradient like dynamics for each robot to generate its trajectory using local information. The resulting equations can be simulated by robust numerical algorithms and executed efficiently. Furthermore, although our method shares a lot of similarities with APF, there are key differences. We add intermittent random perturbations in our dynamics to avoid deadlock, which overcomes the main limitation of APF at a moderately increased computation cost. As in the method of evolving junctions (MEJ) [46, 44], it can be shown that the intermittent dynamics converges to the desired shape much quicker than the continuous white noise perturbations. In addition, the repelling fields from obstacles in APF methods affect the potential everywhere in the domain, while in the proposed method, each robot is viewed as a dynamically moving obstacle to the other robots, and its repelling effect is restricted to a small, local region.

When viewing motion planning for multi-robot systems as transport of distributions, we note that there is recent work inspired by statistical physics [6, 71], in which rigorous error estimates have been obtained between partial differential equations (PDEs) that model the swarm dynamics and the target distribution, enabling desired coverage performance. PDEs have also been used in [16] to generate velocity fields that govern the motion-planning and incorporate collision avoidance. In [15], the controllability properties of the advection-diffusion equation are used to derive conditions on the target probability distribution that guarantee convergence in finite time for certain control inputs. In addition, there are also other stochastic methods for path planning and control [29, 30, 43].

The paper is organized as follows. In §2, we present the basic optimal transport theory and Fokker-Planck equation that inspire us to design the dynamics. In §3, we formulate the continuous problem in terms of a system of SDEs. The discretized problem is described in §4. In §5 we provide numerical simulations of the shape formation problem for different shapes and different size groups. We provide theoretical guarantees for global convergence of the system and collision avoidance, both in the continuous and discrete settings in §6.

2 Relations between SDEs and optimal transport

In this section, we briefly review the connections among stochastic differential equations (SDEs), Fokker-Planck equation, Gibbs distribution, free energy, and optimal transport distance. These relations provide the theoretical foundation on which we design the dynamics for the motion planning of a group of robots.

Let us consider a potential function $\Psi(y)$, in which $y\in{\mathbb{R}}^{Nd}$ represents the locations of $N$ robots in a bounded domain $\Omega\subset{\mathbb{R}}^{d}$. The white noise perturbed gradient flow refers to a SDE

[TABLE]

where $W(t)$ is the standard {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}Nd}-dimensional Brownian motion and $\sigma$ a given constant. Denoting $\rho(t,y)$ the density function for the random variable $Y(t)$, the evolution of $\rho$ is governed by the Fokker-Planck equation according to the classical diffusion theory, i.e.

[TABLE]

By directly plugging in the well known Gibbs distribution

[TABLE]

we see that $\rho^{*}$ is a steady state solution of (2), because it satisfies the equation and it is time independent. In other words, Gibbs distribution is an invariant measure of the system (1). From the exponential form of $\rho^{*}(y)$, we also observe that the density $\rho^{*}(y)$ takes the largest value when $\Psi(y)$ reaches its global minimum.

We would like to remark that in order for the Gibbs distribution to be well-defined, $P$ in (3) must be a finite number. This can be guaranteed by requiring that the potential function $\Psi(y)$ grows quadratically when $|y|$ tends to infinity. In this study, our interest is within a bounded region. Therefore, we can assume that $\Psi(y)$ is defined with such a property at infinity.

The understanding about the connections between Fokker-Planck equation and Gibbs distribution has been greatly enriched in the past few decades, thanks to the new developments in optimal transport theory. In short, defining the 2-Wasserstein distance between any two density functions $\rho^{1}(y)$ and $\rho^{2}(y)$ by

[TABLE]

where the velocity field $v(t,y)$ and density $\rho(t,y)$ satisfy the transport equation

[TABLE]

with boundary values given by

[TABLE]

induces a metric in the probability density space and turns the space into a Riemannian manifold to which one can apply various geometric operations. One of the most impactful results reveals that the Fokker-Planck equation is the gradient flow, with respect to the 2-Wasserstein metric, of a free energy given by

[TABLE]

in which the first term is the potential energy while the second is called entropy [28]. Following the properties of gradient flow, one can prove that Gibbs distribution is the unique attractor of the Fokker-Planck equation (2), and its convergence rate to the Gibbs distribution is exponential, see Theorem 24.7 in [66] for details.

Our idea for motion planning is finding a potential such that the target shape is where the potential attains its global minima if shape formation is the task, or the target distribution is the Gibbs distribution if the goal is to move a group of robots to a given distribution, The exponential convergence of (2) to the Gibbs distribution from any initial distribution forms the basis that guarantees the success of planned motions. The potential is used in conjunction with (1) to create two deterministic dynamics that are used alternately to produce the trajectories for all robots. In the design, we must ensure that (a) the motions are collision-free; (b) there is no deadlock , and; (c) the dynamics converge to the desir ed shape. In the rest of this paper, we use shape formation as the task to illustrate our strategy. Its extension to the distribution case requires simple modifications, which will be omitted in the paper.

Remark: Besides the strategy that we propose here, there are different ways to apply optimal transport theory for motion planning. For example, one may view the robots as a collection of point mass and move them according to the transport equation (5) for which the initial and target distributions are used as $\rho^{1}$ and $\rho^{2}$ respectively. This amounts to finding a velocity $v$ while maintaining point masses throughout the optimization procedure. We do not adopt this view in this paper. Instead, we directly design dynamics based on formulation (1), because the resulting algorithm is simple and efficient in implementation, yet has desirable properties that can be rigorously proved.

It is worth noting that the convergence to the Gibbs distribution for the solution of the Fokker-Planck equation does not have a direct guarantee for the convergence of the SDE (1) to a desirable shape. The subtlety lies in the fact that the convergence for the SDE is only in the distribution sense. The solution of (1) with a positive constant $\sigma$ never settles down asymptotically. To make the solution converge, one has to reduce the value of $\sigma$ gradually to zero, which is precisely the idea used in simulated annealing. However, it is well known that the reduction rate of $\sigma$ must be slower than a logarithmic function in time to avoid local traps. To speed up the convergence, we borrow ideas from intermittent diffusion [9], a stochastic strategy developed for global optimization that can improve the convergence with the probability of success increased to 1 as a geometric sequence, a rate that is much faster than logarithmic functions. Besides, directly applying random perturbations to the motions can cause wasteful jittering effects which we want to avoid in our design.

3 Model setup

Suppose $\Gamma\subset{\mathbb{R}}^{2}$ is a set of spatial locations that form a desired shape, and consider the trajectories of $N$ robots given by

[TABLE]

where $X^{{i}}(t)$ is a curve in ${\mathbb{R}}^{2}$ describing the position of the $i^{\text{th}}$ robot at time $t\geq 0$. The objective is to produce paths $\{X(t)\}_{0\leq t\leq T}$ from an initial state $X(0)=X_{0}$ to a final state $X(T)$ such that $X(T)\in\Gamma$. Our strategy is to design modified gradient flows whose solutions prescribe the path $X(t)$ for all robots.

In order to do so, we first introduce a shape function $F(X)$ that is smooth and has a global minimum only for $X\in\Gamma$. A convenient choice, among many candidates, is the distance function,

[TABLE]

where $\|\cdot\|$ is the Euclidean norm: $\|x\|=\sqrt{x_{1}^{2}+x_{2}^{2}}$. Then, $F(X)$ is a non-negative function achieving its minimum only when $\cup_{i=1}^{N}X^{i}\subset\Gamma$. Figure 1 illustrates the level-sets of $F(X)$ corresponding to two different target shapes.

We also introduce a penalty function $G(X)$ that takes a large value when $X$ exhibits undesirable behavior. In multi-robot systems, one of the main objectives is to ensure that the trajectories are collision-free, meaning the pairwise distances $\|X^{{i}}(t)-X^{j}(t)\|$, $j\neq i$ must be larger than a given positive value $r$, for all $t>0$. For example, we can select the penalty $G(X)$ as the following smooth, “repelling” function that peaks when the pairwise distances are small,

[TABLE]

where the function $\varphi\in C^{1}(0,\infty)$ can be chosen as a decreasing function having the following properties

[TABLE]

For instance, $\varphi(x)=\frac{1}{x}\exp{(\frac{-1}{R^{2}-x^{2}}})$ can be picked to satisfy the requirements. Here, the constant $R>r$ is related to the sensing radius of each robot, and the constant $G_{0}$ is calibrated according to the initial positions of robots and to achieve desirable dynamics. Further constraints on the system, such as obstacle avoidance, can also be easily included. To simplify the presentation, we do not consider obstacle avoidance in this paper.

Combining the shape function (6) with the penalty function (7), we obtain the potential function

[TABLE]

Then the trajectories of the robots are primarily generated by the gradient flow that minimizes $\Psi(X)$, i.e.

[TABLE]

Following it, the robots get to the desired shape when $F(X)=0$, while minimizing $G(X)$ helps to spread out their locations in addition to avoid collisions.

However, the path generated by such a simple gradient flow may suffer a well known shortcoming that the trajectories can get trapped in locations corresponding to local minimizers. To overcome this limitation, we use ideas from intermittent diffusion. More precisely, we intermittently add random perturbations to (10), leading to the following SDEs,

[TABLE]

where $W(t)$ is the standard Brownian motion in ${\mathbb{R}}^{2}$ and $\sigma(t)$ is a piecewise constant function alternating between zero and a positive value, i.e.

[TABLE]

Here we partition $[0,T]$ as $\cup_{k=1}^{K}([T_{k-1},T_{k}])$ with $T_{0}=0$, $T_{N}=T$ and $S_{k}\in[T_{k-1},T_{k}]$.

We want to highlight that the random perturbations are added to the gradient flow to avoid trajectories being trapped at local minimizers. Therefore, the constant $\sigma_{k}$ doesn’t have to be small. This is different from the choice used in simulated annealing, in which the corresponding coefficient, also called temperature, must go to zero asymptotically. The effectiveness of random perturbations can be verified by numerical experiments and comes with guarantees based on optimal transport theory. More precisely, the solution of (11) converges to the global minimizer in the distribution sense according to the Gibbs distribution. The Gibbs distribution is an invariant measure of the system (11), and $\rho^{*}(X)$ takes the largest value when $\Psi(X)$ reaches its global minimum. Further details of the theory are provided in §6.

Unlike many other applications of SDEs, it is important to emphasize that the random portion of the solution $Y(t)$ when $t\in[T_{k-1},S_{k}]$ is not used as the trajectories for the robots due to inefficient jittering motions. Instead, $Y(t)$ is only computed virtually to create the vector $Y(S_{k})$, denoted as $\hat{Y}$ in the rest of the paper, of intermediate positions to move the robots to. Once this position $\hat{Y}$ is computed, we define another objective function

[TABLE]

Using it together with $G(X)$, we create another gradient flow

[TABLE]

In the end, the path $X(t)$ of the robots is generated by alternating between two gradient flows (10) and (14). The implementation of the method is given in the next section.

4 Implementation

The gradient flows and the SDEs presented in the previous section must be solved numerically when calculating the path. We employ the simple Euler scheme to do so in this paper. More precisely, we compute

[TABLE]

where $\Delta t$ is the step size, $\Psi(X)$ takes $F(X)+G(X)$ for (10) and $\hat{F}(X)+G(X)$ for (14) respectively. The SDEs (11) is discretized as

[TABLE]

where $\xi_{n}\in{\mathbb{R}}^{2}$ is a normally distributed random vector generated at each iteration.

As mentioned in the previous section, the path is generated by alternating between (10) and (14). This is implemented by repeating a 2-step strategy. In the first step, the robots are moved, using (14), toward temporary destinations computed by a simulation of (16). After the temporary locations are reached, the second step has the robots follow (10) toward the desired shape. The robots then repeat the two steps until the task is accomplished. Details are presented in Algorithm 1, and the computed descent directions in two different iterations are plotted in Figure 2. Again, we want to re-iterate that $Y_{n}$ is not part of the trajectories. They are computed only virtually to generated the intermediate positions $\hat{Y}$.

This algorithm is a practical modification of the theory developed in the later sections. The main difference is in the diffusion stage of the algorithm, the aim is to produce trajectories that are influenced by both the desired shape and random noise. This procedure is performed offline to save resources; a random path simulated by a robot may be costly even if the ending location is close to the starting position of the robot. Instead of a random path, it is more efficient for the robots to move directly toward these temporary locations. Therefore, in the implementation, each robot moves to its computed destination following a gradient flow, without regard for the shape density function. By doing this, the energy of the system will possibly be increased. This is reflected in the variance of the energy functional in Figure 4.

In §6, we shall prove that the Algorithm 1 generates a guaranteed collision -free path for each robot that converges to the desired shape. Before doing so, we present a few numerical experiments to illustrate the performance in the next section.

5 Numerical Results

In our numerical experiments, we confine the robots in a square domain given by $\Omega=[-M,M]\times[-M,M]$. We assume that each robot has knowledge of its location $X^{i}$, the gradient of the shape function $(\nabla F(X))_{i}=\nabla F_{i}(X)$ , and a sensing radius $R$, meaning that a robot can only detect other robots if they are within a circular region centered at $X^{i}$ with radius $R$. This $R$ is also the parameter we use in $G(X)$:

[TABLE]

We note that this choice of $G$ is different from the function we mentioned in Section 3, demonstrating the flexibility of choosing $G$.

We evaluate the success of the algorithm by determining if the robots are in the desired region, distributed uniformly, and if the nearest -neighbors difference is minimized.

The numerical tests are performed on two shapes. The first shape consisting of points in the set $\Gamma_{1}$, corresponds to a handwritten letter ‘Q’. In this case, the closed loop feature poses difficulties. The second shape consisting of points in the set $\Gamma_{2}$, is a Chinese character, pronounced as ‘JIE’, with multiple complicated strokes and two disconnected components. The initial positions for the robots are either clustered at a corner (demonstrated for shape $\Gamma_{1}$) or randomly distributed in the domain (demonstrated for shape $\Gamma_{2}$). The time evolution, shown in two cases in Figure 3, indicates that the robot trajectories driven by our proposed algorithm drive the robots to the desired shapes without suffering from congestion or getting stuck at local minimizers.

To test the scalability of our algorithm, we varied the size of the robot radius (resulting in different values of $N$). The choice of $N$ is based on a-priori knowledge that there is a global minimum with $N$ robots positioned entirely in the desired shape, determined by trial and error.

From the experiments, we observe that the faster convergence occurs with a random initial configuration that minimizes congestion from the start and provides the robot group immediate access to all sides of the target shape. When robots are initialized in a cluster near one end of the domain, they risk stagnating near the corner of the shape and missing entire sections of the shape unless intermittent diffusion becomes active.

We compared our method to a standard gradient descent with the potential $\Psi$, which is the result of APF. From the energy plots shown in Figure 4, it is clear that gradient descent (APF) alone leaves some robots trapped in local minima. After about 2000 iterations, the congestion caused by the gradient descent iterations is not resolved. Furthermore, the energy decays at a much slower rate than in the iterations produced by Algorithm 1.

6 Mathematical Underpinnings

In this section, we justify theoretically that the generated path using the proposed method can achieve the desired shape while maintaining collision-free motions. We start with the collision-free property first.

Our model determines the trajectories of the robots based on two different gradient flows, (10) and (14) respectively. In both cases, the energy functional $\Psi(X)$ consists of a potential $F(X)$ (or $\hat{F}(X)$) that attracts the robots to the destinations, and the repelling function $G(X)$ that keeps them away from each other. In our theoretical study, it suffices to consider a general potential $F$ that is differentiable, is bounded, and has minimizers only at the desired regions. In this general setting, the governing equation for the path is still given by the gradient flow presented in (10).

6.1 Continuous time collision avoidance

Recall that the location of the $i^{\text{th}}$ robot is given by $X^{{i}}$, and the set of admissible robot coordinates is $\mathcal{X}$, where

[TABLE]

We note that the repelling function $G(X)$ satisfies (7), for a function $\varphi$ satisfying (8), which implies $G(X)$ is a $C^{1}$ function on $\mathcal{X}$.

Let $X_{0}\in\mathcal{X}$ be the initial robot locations with smallest pairwise distance given by

[TABLE]

Suppose that the decreasing function $\varphi(x)$ satisfies

[TABLE]

This can be achieved for $\varphi(x)$ satisfying the property

[TABLE]

Then we have the following theorem.

Theorem 6.1**.**

For any trajectory $X(t)=(X^{1}(t),\ldots,X^{N}(t))$, $N>1$ generated by (10) with initial position $X(0)=X_{0}\in\mathcal{X}$, the inequality

[TABLE]

*holds for all $t>0$. *

Proof 6.2**.**

The function $\Psi(X(t))$ is non-increasing along the solution of (10) since it satisfies

[TABLE]

Assume there is a time $t^{*}>0$ such that $\|X^{{i}}(t^{*})-X^{j}(t^{*})\|^{2}\leq{r^{2}}$ for some $i,j\neq i$, then

[TABLE]

Here we used that $F$ is non-negative by construction. This is a contradiction, because $\Psi(X(t))$ is non-increasing, so we must have $\Psi(X({t}^{*}))\leq E_{0}$.

6.2 Discrete time collision avoidance

Equation (10) and (14) are solved in discrete time using the iterations

[TABLE]

where $X_{n}\simeq X({t_{n}})$ and $t_{n}=n\Delta t$ for some fixed time step $\Delta t>0$. It is known that the Euler scheme converges to the continuous solution if $\nabla\Psi$ is $L-$Lipschitz continuous in space. This ensures no collision in the discrete case when the step size is small enough. In the next theorem, we present such a result, and prove it by using a standard argument from [53].

Theorem 6.3**.**

*Suppose $\Psi\in C^{1}(\mathcal{X})$ is a positive function that is bounded below and $\nabla\Psi$ is $L-$ Lipschitz continuous in space. Then, if $\Delta t\leq\frac{1}{L}$, one step of the gradient method (20) will not increase the objective function $\Psi$, that is $\Psi(X_{n+1})\leq\Psi(X_{n})$. *

Proof 6.4**.**

Denote the Euclidean inner product by $\langle X,Z\rangle=\left(\sum_{i=1}^{N}X^{i}\cdot Z^{i}\right)^{1/2}$ For $X,Z\in\mathcal{X}$, we can express $\Psi(Z)-\Psi(X)$ by

[TABLE]

This results in

[TABLE]

Taking $Z^{i}=X^{i}_{n+1}=X_{n}^{i}-(\nabla\Psi(X_{n}))_{i}\Delta t$, we have

[TABLE]

*Therefore $\Psi(X_{n+1})\leq\Psi(X_{n})$ if $\Delta t\leq\frac{2}{L}$. *

We remark that $\nabla\Psi(X)=F(X)+G(X)$, with $G(X)$ being defined through $\varphi(x)$, satisfies the $L$-Lipschitz condition in the domain of interest, because $\varphi(x)$ is a $C^{1}$ function on the closed interval $[r,2M]$. The Lipschitz constant $L$ depends on the choice of $\varphi$, the size of computational domain $\Omega$, and the number of robots in the group.

Corollary 6.5**.**

The discrete trajectory $X_{n}$ computed by (15) satisfies

[TABLE]

*for all $n\geq 0$, provided $E_{0}=\Psi(X_{0})<G_{0}\varphi(m_{0})$. *

The proof of this corollary follows directly from the proof of Theorem 6.1 and the result of Theorem 6.3.

6.3 Convergence to the global minima in probability

As described in the model, the goal of introducing (14) is to move the robots to the intermediate locations generated by the SDEs (11). Therefore, the convergence of the trajectories to the desired shape means that the solutions of (11) march to the global minima of $\Psi(x)$, which is guaranteed by the theory of optimal transport. More precisely, the idea of combining (10) and (14) comes from the intermittent diffusion. Together, the dynamics can be equivalently described by a uniform formula given in (11), in which (10) is performed when $\sigma=0$, and (14) reaches the same spatial locations as (11) when $\sigma$ is not zero. Hence the question of whether or not $X(t)$ converges to the desired shape can be investigated by examining the distributions of trajectories in (11).

We recall from Section 2 that the probability density function $\rho(y,t)$ of the stochastic process $Y(t)$ from (11) evolves according to the Fokker-Planck equation, which is a transport equation when $\sigma=0$, and a diffusion equation when $\sigma>0$. In the diffusion case, the asymptotic solution, also called the steady state, is the Gibbs distribution defined in (3), suggesting the probability that $X(t)$ is within the attractive neighborhood, denoted by $\hat{U}$, of the global minimum of $\Psi$ is positive if $t$ is large enough. Here $\hat{U}$ is defined as the neighborhood of $\Gamma$, in which the trajectory of the gradient flow (10) with any initial configuration $X_{0}\in\hat{U}$ satisfies $\lim_{t\rightarrow\infty}\mu(X^{i}(t))=0$, where $\mu$ is the distance function to $\Gamma$ defined in (6). By the subsequent gradient flow (10), $X(t)$ remains inside of the target $\Gamma$ or moves arbitrarily close to it. This suggests that there is a positive probability that $X(t)$ is within a small neighborhood of $\Gamma$ after one cycle of intermittent diffusion ($\sigma$ taking positive and then zero values once) is also positive. Repeating the cycle of intermittent diffusion, we obtain the following convergence theorem.

Theorem 6.6**.**

Suppose $\Psi(x)$ attains its global minima on a set $\Gamma$ of positive Lebesgue measure, and let $U\subseteq\hat{U}$ be a small neighborhood of $\Gamma$. Then for any $0<\eta<1$ there exist constants $T^{*}>0$, $\sigma_{0}>0$ and $K_{0}>0$ such that if $(T_{i}-S_{i})>T^{*}$, $\sigma_{i}<\sigma_{0}$ for $1\leq i\leq K$ and $K>K_{0}$, the solution $X_{opt}$ calculated by Algorithm 1 satisfies

[TABLE]

*where $\mathbb{P}$ is the probability function. *

The proof of this theorem essentially follows the same steps as the proof in [9], in which the convergence of the density is considered in the $L^{1}$ sense by using the Csiszar-Kullback inequality (see Remark 22.12 in [66]). For the completeness of this paper, we present a sketch of the proof modified to guarantee convergence in the $W_{2}$ sense.

Proof 6.7**.**

By the construction of $\Phi(X)$, its value is non-negative and reaches the minimum [math] only if $X\in\Gamma$. This suggests that the Gibbs distribution $\rho^{*}(y)$ attains its maximum when $y\in\Gamma$. Since $\Gamma$ has a positive Lebesgue measure, so does $\hat{U}\supset\Gamma$. Hence there exits a positive constant $\nu$ such that

[TABLE]

In fact, $2\nu\in(0,1]$ approaches $1$ when $\sigma$ tends to [math] according to the property of Gibbs distribution.

From Theorem 24.7 and discussions following in Example 24.8 and Remark 24.12 in [66], we have

[TABLE]

where $C$ and $\lambda$ are constants, $\lambda$ is related to the well-known Log-Sobolev inequality (see Definition 21.1 in [66]), and $\rho(y,t)$ is the solution of Fokker-Planck equation (2) with $\sigma>0$. This implies that there exists a constant $T^{*}>0$ such that

[TABLE]

for arbitrary $t>T^{*}$. It suggests that there is a positive probability greater than $\nu$, that $Y(T_{i})$ is in $\hat{U}$ when $(T_{i}-S_{i})>T^{*}$ in the virtual diffusion process (11). Because $\hat{Y}=Y(T_{i})$ is used in (14), we conclude that the initial position $X(T_{i})=\hat{Y}$ for the gradient flow (10) belongs to $\hat{U}$ with a positive probability. The trajectories $X(t)$ that start from $X(T_{i})$ are pushed into the neighborhood $U$ exponentially fast due to the definition of $\hat{U}$ and gradient flow properties.

In other words, the probability that the trajectory $X(t)$ does not end in $U$ is at most $(1-\nu)$ every time when one virtual diffusion and gradient flow cycle is completed. If such a cycle is performed $K$ times, the probability that $X(t)$ does not reach $U$ is $(1-\nu)^{K}$. Since $0\leq(1-\nu)<1$, there exist a $K_{0}>0$ such that $(1-\nu)^{K}<\eta$ for any $K>K_{0}$. Therefore,

[TABLE]

*which completes the proof. *

The proof also indicates that $(1-\eta)$ approaches $1$ in the manner of

[TABLE]

which forms a geometric sequence in term s of $K$. This is a much quick convergence rate than the logarithm function owned by the simulated annealing. We would like to point out that the convergence result presented in Theorem 6.6 is in the sense of probability, which is different from the usual deterministic convergence results given in the $L^{p}$-norm or maximum norm, but our numerical experiments show that $X_{opt}$ always reaches the desired shape $\Gamma$ without failure if the parameters are selected properly.

7 Conclusions and Future Work

We present a motion planning strategy for a large group of robots to accomplish shape formation, one of the fundamental tasks in many applications that employ multi-robot systems. Typical challenges include how to avoid collisions and deadlocks in motion planning and how to achieve the desired shape with assurance. Those challenges become more significant for large groups of of robots and robots with low functionality. In our method, we calculate the individual robot trajectories by alternating two gradient flows that involve an attractive potential, a repelling function , and a process of intermittent diffusion. The potential attracts robots to form the targeted shape, while the repelling function is designed to ensure collision-free motions. The intermittent diffusion, originally a stochastic approach but here realized by deterministic means, overcomes situations with deadlocks. Our strategy is inspired by recent developments in the theory of optimal transport which in turn provides the basis for theoretical guarantees of collision avoidance and global convergence. Numerical experiments confirm that the proposed algorithm is simple, yet effective in achieving the desired objectives.

The presentation here in the two-dimensional setting can be extended to higher dimensions with straight forward adaptations. The proposed strategy can also be adapted to accommodate inhomogeneous multi-robot systems, in which robots may have different functionalities. In this scenario, the differences among robots must be reflected throughout the selections of the potential functions, including both $F(x)$ and $G(x)$. On the technical side, this may not be easy to accomplish and it is in our plan for further investigation.