{'fraction_non_alphanumeric': 0.04764950431723697, 'fraction_numerical': 0.010375581849838326, 'mean_word_length': 4.444766879473786, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Transporting Robotic Swarms via Mean-Field Feedback Control

With the rapid development of AI and robotics, transporting a large swarm of networked robots has foreseeable applications in the near future. Existing research in swarm robotics has mainly followed a bottom-up philosophy with predefined local coordination and control rules. However, it is arduous to verify the global requirements and analyze their performance. This motivates us to pursue a top-down approach, and develop a provable control strategy for deploying a robotic swarm to achieve a desired global configuration. Specifically, we use mean-field partial differential equations (PDEs) to model the swarm and control its mean-field density (i.e., probability density) over a bounded spatial domain using mean-field feedback. The presented control law uses density estimates as feedback signals and generates corresponding velocity fields that, by acting locally on individual robots, guide their global distribution to a target profile. The design of the velocity field is therefore centralized, but the implementation of the controller can be fully distributed -- individual robots sense the velocity field and derive their own velocity control signals accordingly. The key contribution lies in applying the concept of input-to-state stability (ISS) to show that the perturbed closed-loop system (a nonlinear and time-varying PDE) is locally ISS with respect to density estimation errors. The effectiveness of the proposed control laws is verified using agent-based simulations.


I. INTRODUCTION
Transporting a large robotic swarm to form certain desired global configuration is a fundamental question for a wide range of applications, such as scheduling transportation systems [1] and employing nanorobots for drug delivery. Swarm robotic system provides superior robustness and flexibility, but also poses significant challenges in its design [2]. We pursue the design of robotic swarms as a control problem, and propose a control theory based design framework.
The major difficulty of controlling such large-scale systems results from that the control mechanism is expected to be scalable and satisfy the robots' own kinematics while their collective behaviors should be predictable and controllable. Existing work has revealed two different philosophies, termed as bottom-up and top-down, respectively [3]. Probabilistic finite state machines [4] and artificial physics [5] are two classic representatives dating back to 1980s which follow the bottomup philosophy and are known to be decentralized and scalable. However, evaluation of their stability and performance quickly becomes intractable when we increase the swarm size. In 2000s, graph theory was introduced into the multi-agent system community and has seen successful applications in designing coordination protocols [6]. Nevertheless, it needs to address the dimensionality issue due to large-size matrices arising in swarm robotic systems. In recent years, top-down design has received increased research interests, which usually employs compact models to describe the macroscopic behaviors. The challenge lies in the appropriate decomposition of the global control strategy into local commands. Markov chains approach is one representative that uses abstraction-based models for macroscopic descriptions, which partitions the workspace into a grid, over which it defines a probability distribution and designs the movement between cells to govern the evolution of the distribution [7], [8]. The major drawback is that the robots' dynamics are not considered. Potential games approach is another top-down approach that adopts a game theoretic formulation, which suggests to decompose the global objective function into local ones that align with the global objective [9]. However, finding such a decomposition is problem dependent and difficult in general.
Our work is inspired by the recent top-down design that uses PDEs for macroscopic descriptions. There exist mainly two types of PDE models in the literature. The first type is motivated from the fact that certain discretized PDEs match the dynamics of graph-based coordination algorithms [10]- [14]. Boundary control and backstepping design are popular techniques for such models. Although boundary actuation is intriguing, it has difficulty for higher-dimensional extension. Our work adopts the other type that is known as meanfield PDEs [15]- [22]. These mean-field models fill the gap between individual dynamics and their global behavior with a family of ordinary/stochastic differential equations that describe the motion of individual robots, and a PDE that models the time-evolution of the mean field (i.e. their probability distribution). A typical design process starts with specifying the task using the macrostate of the PDE, and then computes local motion commands for individuals. In [19], the authors formulate an optimization problem for a set of advection-diffusion-reaction PDEs to compute the velocity field and switching rates. In [18], the authors present a PDE-constrained optimal control problem, and microscopic control laws are derived from the optimal macroscopic description using a potential function approach. These optimization-based approaches are however computationally expensive, open-loop, and may be unstable in the presence of unknown disturbance. Mean field games incorporate the mean-field idea into large population differential games and obtain a compact model with two coupled PDEs [16]. The control strategy in our work is inspired by the recent idea that uses meanfield feedback to design appropriate velocity fields [20]- [22]. Such control laws can be computed efficiently and be formally proved to be convergent. Nevertheless, the works [20], [21] are restricted to deterministic individual dynamics, while stochasticity is ubiquitous in practice, caused possibly by sensor and actuator errors, or the inherent avoidance mechanism of the robots. The control law proposed in [22] applies to the stochastic case, but their focus is on the controllability property. In practice, mean-field feedback control relies on estimating the unknown density, which causes robustness issues in terms of estimation errors. Our work distinguishes from [20]- [22] in that we will consider the stochastic case, present general results for its solution property, and study the robustness issue of the proposed control law.
In particular, we study the problem of mean-field feedback control of swarm robotic systems modelled by PDEs. We design velocity fields (from which individual control commands can be derived) by using the real-time density as feedback signals, such that the density of the closed-loop system evolves towards a desired one. Our contribution includes three aspects. First, we present general results for the solution property (well-posedness, regularity and positivity) of the PDE system. Second, we propose mean-field feedback laws for robotic swarms that involve stochastic motions and apply the notion of input-to-state stability (ISS) to prove that the closed-loop system is locally ISS with respect to density estimation errors. Third, in terms of theoretical contribution to PDE control systems, our results apply the concept of ISS to nonlinear and time-varying PDEs with unbounded operators. Most existing work that studies ISS for PDE systems is restricted to the linear case or one-dimensional case (e.g. boundary control). However, the control input in our problem is a vector field that couples with the system state (which thus makes the perturbed closed-loop system nonlinear), and acts on the system through unbounded operators.
The rest of the paper is organized as follow. Section II introduces some preliminaries and useful lemmas. Problem formulation is given in Section III. Section IV is our main results, in which we present general results for the solution property of the PDE system, present a mean-field feedback control law and then study its robustness issue with respect to density estimation errors. Section V performs an agent-based simulation to verify the effectiveness of the control law. Section VI summarizes the contribution and points out future research.

A. Notations and useful lemmas
Let E ⊂ R n be a measurable set and k ∈ N.
with the norm f L ∞ (E) . We use D α f to represent the weak derivatives of f for all multi-indices α of length |α|.
|α|≤k D α f L p (E) < ∞}, endowed with the norm · W k,p (E) . Analogously, W k,∞ (E) is defined, equipped with the norm · W k,∞ (E) . We also denote H k = W k,2 . The spaces W k,p (E) are referred to as Sobolev spaces.
Let ω be a bounded and connected C 1 -domain in R n and T > 0 be a constant. Denote by ∂ω the boundary of ω. Set Ω = ω × (0, T ). For a function u(x, t) : Ω → R, we call x the spatial variable and t the time variable. We denote ∂tu = ∂u/∂t and ∂iu = ∂u/∂xi, where xi is the i-th coordinate of x. The gradient and Laplacian of a scalar function f are denoted by ∇f and ∆f , respectively, and the divergence of a vector field F is denoted by ∇·F . The differentiation operation of these operators are only taken with respect to the spatial variable x if f and F are also functions of t.
The norm on M is defined by Throughout this paper, when we say f (x) is a density function on ω, we mean f (x) is a probability density function, i.e. f (x) ≥ 0, ∀x ∈ ω and ω f (x)dx = 1.
Lemma 1: (Poincaré inequality [23]). For p ∈ [1, ∞) and ω, a bounded connected open set of R n with a Lipschitz boundary, there exists a constant C depending only on ω and p such that for every function f ∈ W 1,p (ω), |ω| ω f dx, and |ω| is the Lebesgue measure of ω.

B. Input-to-state stability
Input-to-state stability is a stability notion widely used to study stability of nonlinear control systems with external inputs [24]. We introduce its extension to infinite-dimensional systems presented in [25]. Let (X, · X ) and (U, · U ) be the state space and the space of input values, endowed with norms · X and · U , respectively. Denote by P C(I; Y ) the space of piecewise rightcontinuous functions from I ⊂ R to Y , equipped with the standard sup-norm. Define the following classes of comparison functions: We use the following axiomatic definition of a control system [25]. Definition 1: The triple Σ = (X, Uc, φ), consisting of the state space X, the space of admissible input functions Uc ⊂ {f : R+ → U }, both of which are linear normed spaces, equipped with norms · X and · Uc , respectively, and of a transition map φ : is called a control system, if the following properties hold: Here, φ(t, s, x, u) denotes the system state at time t ∈ R+, if its state at time s ∈ R+ was x ∈ X and the input u ∈ Uc was applied.
The control system is called input-to-state stable (ISS), if in the above definition ρx and ρu can be chosen equal to ∞. If Uc = P C (R+; U ), then u Uc = sup 0≤s≤∞ u(s) U , and due to the causality property of Σ, one can obtain an equivalent definition of (L)ISS by replacing (1) with the following inequality [25]: To verify the ISS property, Lyapunov functions can be exploited.
if there exist constants ρx, ρu > 0, functions ψ1, ψ2 ∈ K∞, χ ∈ K, and a continuous positive definite function W , such that: where the derivative of V corresponding to the input u is given bẏ If in the previous definition D = X, ρx = ∞ and ρu = ∞, then the function V is called ISS-Lyapunov function. If Uc = P C (R+; U ), then condition (ii) in Definition 3 can be replaced by the following condition (ii) due to the causality property [25]:

III. PROBLEM FORMULATION
This work studies the transport problem of robotic swarms. Specifically, we want to design velocity commands for individual robots such that the swarm evolves to certain global distribution. The robots are assumed homogeneous whose motions satisfy: where N is robots' population, Xi ∈ ω is the position of the i-th robot, v(Xi, t) = (v1, . . . , vn) ∈ R n is the velocity field that acts on the robots, Bt ∈ R n is an n-dimensional Wiener process which represents stochastic motions, and 2σ(Xi, t) ∈ R is the standard deviation of the stochastic motion at position Xi. Their macroscopic state can be described by the following meanfield PDE, also known as the Fokker-Planck equation, which models the evolution of the swarm's mean-field density p(x, t) on ω: where n is the unit inner normal to the boundary ∂ω, and p0(x) is the initial density. The last equation is the reflecting boundary condition to confine the swarm within the domain ω.
Remark 1: We point out that (4) holds regardless of the number of robots. However, if N is small, using the swarm's (probability) density to represent its global state doesn't make much sense. Hence, we usually assume N is large. Note that (3) and (4) share the same set of coefficients, which means that the macroscopic velocity field we design for the PDE system can be easily transmitted to individual robots. Note that individual robots need to derive their own lowlevel controller to track the reference velocity command, which can however be done in a distributed way. The velocity tracking problem has been widely studied in literature especially for mobile robots, and hence is not studied in this paper.
Problem 1: Given a desired density p * (x), we want to design the velocity field v(x, t) such that the solution of (4) converges to p * (x).

A. Well-posedness and regularity
Before presenting the velocity laws, we shall study the solution property of (4), including its well-posedness, regularity, mass conservation and positivity preservation. We point out that (4) is a special case of the so-called conormal derivative problem for parabolic equations of divergent form [23]. Relevant results from Chapter VI in [23] are summarized in Appendices. Considering the initial/boundary value problem (4), we have the following result for its weak solutions (see Definition 4 in Appendices).

Remark 2:
The regularity conditions for σ and p0 in (5) and (6) can be easily satisfied, while the regularity condition for v depends on the velocity field we design. We shall further study this problem in subsequent sections. We point out that such a weak solution has L 2 spatial derivatives. It can be derived from Definition 4 and (23) that such a weak solution satisfies ∂tp ∈ L 2 (Ω), which implies that p(t, ·) is absolutely continuous from (0, T ) to L 2 (ω). This time regularity will enable us to use Lyapunov functions to study its stability.

B. Exponentially stable mean-field feedback control
First, we present a mean-field feedback law with exponential convergence assuming p(x, t) is available. Given a desired density Our main idea is to design v(x, t) such that Φ(x, t) satisfies the following diffusion equation: where α(x, t) > 0 is the diffusion coefficient. It is known that under mild conditions on α(x, t), the solution of (7) evolves towards a constant function in ω, which will be 0 because Φ is the difference of two density functions and, for any t, The idea of using diffusion/heat equations for designing velocity fields is originally from [20]. Our work extends the original work in three aspects. First, we generalize the design to PDEs that contains stochastic motions and rigorously study its solution property to justify the Lyapunov-based stability analysis. Second, the control law given in [20] can be problematic if the density becomes zero. We will show how to avoid this issue by appropriately constructing the density estimate in the mean-field feedback law. Third, we continue to study the robustness of this modified feedback law with respect to density estimation errors (which includes not only the inherent error of any estimation algorithm, but also the "artificial error" introduced to ensure that the feedback law remains bounded).
Our first result is to enhance the stability result in [20] assuming that the density is strictly positive and can be perfectly measured.
Proof: Substituting (9) into (4), we obtain the closed-loop PDE which is a diffusion equation with Neumann boundary condition. Its exponential convergence is well-known [26]. We include the proof for completeness. Consider a Lyapunov function where we used divergence theorem for the third equality, the boundary condition for the forth equality, and Poincaré inequality (for which we also use the fact that ω Φdx ≡ 0) for the inequality, and C > 0 is a constant depending on ω. Since α has a uniform positive lower bound, we obtain exponential stability.
The control law (9) generates a dynamic velocity field on ω based on the real-time density (in a centralized way), where α is a design parameter for adjusting the local velocity magnitude. By following this velocity command, the swarm density is guaranteed to evolve towards p * (x). In implementation, each robot computes its velocity command v(Xi, t) using only function values around its position Xi and then derives its own velocity tracking controller in a distributed way. Thus, this control strategy is computationally efficient and scalable to swarm sizes. Note that for (9) to be well-defined, we require its denominator p > 0. This problem will be addressed when we replace p with an estimated density later.

C. Mean-field feedback control using density estimates
The exponential stability result requires that p(x, t) is known and positive for all x ∈ ω. In this section, we use kernel density estimation (KDE) to obtain an estimate of p that is always positive, and study its robustness with respect to estimation errors.
KDE is a non-parametric way to estimate an unknown density [27]. The robots' positions {Xi(t)} N i=1 can be seen as a set of N samples of the common density p(x, t). The density estimator is given bŷ where K(x) is a kernel function chosen to be the Gaussian kernel and h is the bandwidth, usually chosen as a function of N such that limN→∞ hN = 0 and limN→∞ N hN = ∞. (Many boundary correction methods exist for refining the density estimate to have compact support [27], so we shall not worry about this issue.) These estimates (and their derivatives) are in general uniformly consistent in the sense that limN→∞ pN − p L ∞ = 0 with probability 1.
With the density estimate, the control law is changed to which is well-defined sincep > 0 with our choice of kernels. Remark 3: Since the Gaussian kernels are positive C ∞ functions and the Wiener processes have continuous paths, we have thatp ∈ C ∞ (ω)×C([0, T ]). Moreover, since ω is bounded and N is finite, we have that inf (x,t)∈Ωp (x, t) > 0 and sup (x,t)∈Ω ∂ k ip (x, t) < ∞ for any k ∈ N as long as we fix h.
with ∂ivi ∈ L ∞ (Ω), i.e. v satisfies the regularity conditions (5) and (6) in Theorem 1. Now we study the robustness issue in terms of density estimation errors. Such errors can arise from not only the inherent error of any estimation algorithm, but also some "artificial error" we impose onp to ensure that the feedback law (11) remains bounded. Since p > 0, we can define (x, t) :=p(x, t)/p(x, t) − 1, or equivalentlyp = p 1 + ). Then = 0 if and only ifp = p, for which we view as estimation errors. We also have > −1 sincep > 0. Our idea is to treat a functional of (x, t), denoted by d(t) (defined later), as external input and establish ISS property with respect to d(t). In this way, the perturbed closed-loop system will be bounded by a function of d(t) and be asymptotically stable when d(t) = 0.

By defining
we can rewrite (14) in a form of an abstract bilinear control system: where A1, B1, B2 are linear operators and A2 is bilinear. Hence, this system is essentially nonlinear. To study its ISS property, we first present the following theorem which exploits Lyapunov functions to verify the ISS property for nonlinear and time-varying infinitedimensional control systems. Theorem 3: Let Σ = (X, Uc, φ) be a control system, and x ≡ 0 be its equilibrium point. Assume for all u ∈ Uc and for all s ≥ 0 a functionũ, defined byũ(τ ) = u(τ + s) for all τ ≥ 0, belongs to Uc and ũ Uc ≤ u Uc . If Σ possesses an (L)ISS-Lyapunov function, then it is (L)ISS.
Proof: The proof is included in the Appendices. It is based on the proof in [25], but extends it to time-varying control systems.
Remark 4: We point out that in the development of the ISS notion, there is no reference to specific notion of solution. Instead, it is based on the concept of an abstract control system defined in Definition 1. Hence, as long as the notion of solution (e.g. weak, mild, strong, classical) of the infinite-dimensional problem is selected such that the properties (especially the continuity and semigroup property) in Definition 1 are satisfied, and as long as the derivativė Vu defined in Definition 3 exists for almost all t ≥ 0, then we can exploit (L)ISS-Lyapunov functions to study the ISS property. In the existing literature, the notion of mild solution, defined using C0 semigroups, is more commonly used [25]. It however may lose many useful structures and properties of the specific equation (especially PDEs) under study, and requires more complicated techniques to characterize time-varying and nonlinear systems. The notion of weak solution adopted in this work is standard for parabolic PDEs in the PDE literature, which satisfies the properties in Definition 1 when it exists and is unique. (In fact, for linear PDEs with time-independent coefficients, these two notions are equivalent; see page 105 in [28].) By using weak solutions, we are able to obtain the necessary solution properties for studying the ISS property of system (14) even with time-varying coefficients and unbounded control operators, which could have been difficult to study if using mild solutions.
Theorem 4 indicates that the closed-loop system using control law (11) remains bounded as long as the density estimation error satisfies the constraint (17). Note that we can drop the constraint (17) and obtain (global) ISS property if we let α = σ, in the price of possibly slow convergence since σ is usually small. Therefore, the design parameter α yields a trade-off between convergence speed and robustness in the sense that a larger α produces faster convergence but also reduces the convergent domain in (17). Nevertheless, one can also increase the probability of satisfying (17) by increasing the number of robots N , which is the consistency property of KDE.
Remark 5: The term ∇ in (16) is caused by the gradient operator ∇ in (11), which is unavoidable because the gradient operator is unbounded, that is, we cannot bound ∇ using . In fact, also acts on the system (12) through the divergence operator ∇· on the righthand side. The reason why it does not show up in (16) is that in the formulation of weak solutions (23), by using integration by parts, the divergence actually acts on the test functions η (corresponding to Φ in (18)) and produces a boundary term on S(Ω) which eventually disappears due to the reflecting boundary condition in (4). It would have been difficult to deal with the unbounded divergence operator if we adopt mild solutions for ISS analysis.
Remark 6: We shall clarify that our swarm control strategy essentially consists of two parts: centralized velocity field design and distributed velocity tracking. This work mainly focuses on the design of velocity field, which is centralized because it requires knowing the positions of all the robots to estimate their density. This can be implemented by a monitoring system which collects the robots' positions to estimate the global density and then broadcasts the velocity field to the robots. The individual velocity tracking control is however distributed because each robot receives its reference velocity command and then derive its own control signal accordingly (which is a well-studied control problem for mobile robots).

V. SIMULATION STUDIES
An agent-based simulation using 1024 robots is performed on Matlab to verify the proposed control law. We set ω = (0, 1) 2 , σ = 0.0005 and α = 0.03. Each robot is simulated by a Langevin equation (3) under the velocity command (11). The robots' initial positions are drawn from a uniform distribution. The desired density p * (x) is illustrated in Fig. 1a (which is C ∞ and lower bounded by a very small positive constant due to smoothing preprocessing). KDE is used to obtain the density estimatep(x, t), in which we set h = 0.045. Numerical computation of the velocity field (11) is based on finite difference. Specifically, ω is discretized into a 64 × 64 grid, and the time difference is 0.02s.
(b) The convergence error. , the estimated densityp(x, t) of the swarm, and the velocity field v(x, t) generated by (11), which suggests that the swarm is able to evolve towards the desired configuration. The convergence error p−p * is given in Fig. 1b, which shows that the error converges exponentially to a small neighbourhood around 0 and remains bounded, which verifies the ISS property of the proposed algorithm.

VI. CONCLUSIONS
This paper studied controlling the density of a swarm of robots using velocity fields that are computed in a feedback manner. The resulting closed-loop system was proven to be LISS with respect to density estimation errors. The presented framework filled the gap between local kinematics of individual robots and their emergent behaviors in swarm robotic systems. It was top-down and computationally efficient. With the feedback technique, the global performance was guaranteed to be convergent and robust to estimation errors when performing in real-time. Our future work includes studying the distributed density estimation problem and considering more general robotic dynamics.

A. Conormal derivative problems
Equation (4) is a special case of the so-called conormal derivative problem for parabolic equations of divergent form [23]. We summarize (and modify appropriately) main results from Chapter VI [23].
We use the same notations as in Section II-A. We follow the summation convention that any term with a repeated index i is summed over i = 1 to n. For example, bi∂iu = n i=1 bi∂iu. For bounded functions aij, bi, ci, and c0 in ω, define the operator Lu := −∂tu + ∂i(aij∂ju + biu) + ci∂iu + c0u. We always assume that {aij} is uniformly elliptic, i.e., for some positive constant λ, aij(x, t)ξiξj ≥ λ|ξ| 2 for any (x, t) ∈ Ω and any ξ ∈ R n . where ν = (ν1, · · · , νn) is the unit inner normal to the boundary.
For given functions fi and g on Ω, ϕ on ω, and ψ on S(Ω), the conormal derivative problem has the following form: In this paper, we only need the case fi = g = 0 in Ω and b0 = ψ = 0 on S(Ω). We present the general form for completeness. Take any test function η ∈ C 1 (Ω). Multiplying the first equation of (22) by −η and integrating by parts, we obtain where ds is the area form of the boundary ∂ω.
We note that a weak solution is simultaneously a weak subsolution and a weak supersolution. We now discuss the well-posedness and some properties of the weak solution. The following result is based on Theorem 6.38 and Theorem 6.39 in [23].
We have the following energy identity for the weak solution u ∈ M: for almost every t ∈ (0, T ], The proof is by an approximation argument, i.e., take η = u ∈ M and show it is the limit in M of a sequence of H 1 functions [23]. From now on, we assume ω is a connected domain. The following result is based on Theorem 6.43 in [23], which is for subsolutions. Theorem 6: (Strong maximum principle [23]). Assume fi = g = 0 in Ω, ψ = 0 and b0 ≤ 0 on S(Ω), ϕ ∈ L ∞ (ω), and, for any v ∈ C 1 (Ω) with v ≥ 0, Ω (−bi∂iv + c0v)dxdt ≤ 0.
Let u ∈ M be a weak subsolution of the problem (22). Then, Moreover, u is constant if the equality holds at some (x, t) ∈ ω × (0, T ]. If ∂ibi ∈ L ∞ (Ω), the condition (25) can be substituted by its pointwise form ∂ibi + c0 ≤ 0 in Ω, and is not needed if we compare u with 0. Specifically, we have the following positivity result.

1
• ψ2 • χ ∈ K. Our estimates hold for arbitrary control u : u Uc ≤ ρu; thus, combining (26) and (27), we obtain the claim of the theorem. To prove the ISS of Σ from existence of ISS-Lyapunov function, one can argue as above but with ρx = ρu = ∞.