# Density Stabilization Strategies for Nonholonomic Agents on Compact
Manifolds

Karthik Elamvazhuthi and Spring Berman KE and SB were partially supported by
ONR Young Investigator Award N000141612605. KE was also partially supported by
AFOSR MURI FA9550-18-1-0502 and ONR Grant No. N000142012093.Karthik
Elamvazhuthi is with the Department of Mechanical Engineering, University of
California, Riverside, CA 92521, USA<EMAIL_ADDRESS>Berman
is with the School for Engineering of Matter, Transport and Energy, Arizona
State University, Tempe, AZ, 85281 USA<EMAIL_ADDRESS>
###### Abstract

In this article, we consider the problem of stabilizing a class of degenerate
stochastic processes, which are constrained to a bounded Euclidean domain or a
compact smooth manifold, to a given target probability density. This
stabilization problem arises in the field of swarm robotics, for example in
applications where a swarm of robots is required to cover an area according to
a target probability density. Most existing works on modeling and control of
robotic swarms that use partial differential equation (PDE) models assume that
the robots’ dynamics are holonomic, and hence, the associated stochastic
processes have generators that are elliptic. We relax this assumption on the
ellipticity of the generator of the stochastic processes, and consider the
more practical case of the stabilization problem for a swarm of agents whose
dynamics are given by a controllable driftless control-affine system. We
construct state-feedback control laws that exponentially stabilize a swarm of
nonholonomic agents to a target probability density that is sufficiently
regular. State-feedback laws can stabilize a swarm only to target probability
densities that are positive everywhere. To stabilize the swarm to probability
densities that possibly have disconnected supports, we introduce a semilinear
PDE model of a collection of interacting agents governed by a hybrid switching
diffusion process. The interaction between the agents is modeled using a
(mean-field) feedback law that is a function of the local density of the
swarm, with the switching parameters as the control inputs. We show that under
the action of this feedback law, the semilinear PDE system is globally
asymptotically stable about the given target probability density.  The
stabilization strategies with and without agent interactions are verified
numerically for agents that evolve according to the Brockett integrator; the
strategy with interactions is additionally verified for agents that evolve
according to an underactuated system on the sphere $S^{2}$.

###### Index Terms:

Mean-field control, Hypoelliptic operators, Nonholonomic systems, Multi-agent
systems, Swarms, Semilinear PDE

##  I Introduction

In recent years, there has been much work on the construction of decentralized
control laws for multi-robot systems using mean-field models [21], in which a
large collective of agents is treated as a continuum. Mean-field based
approaches for modeling and control have been used in the multi-agent control
and swarm robotics literature for problems such as consensus [59], flocking
[30], task allocation [8, 45], and cooperative transport [61]. Similar
problems have also been considered over the last two decades in the control
and mathematics literature in the context of mean-field games [38, 31], mean-
field control [28, 51, 12], and optimal transport [7].

One advantage of mean-field based approaches for decentralized control is that
the constructed control laws are identity-free, that is, the control laws do
not depend on the agents’ identities. The identity-free nature of the control
laws simplifies some aspects of their implementation on large swarms of
homogeneous agents when compared to control laws that are identity-dependent.
For instance, suppose that a central supervisor observes the states of the
agents via an overhead camera and uses these measurements to update state-
feedback control laws, which it periodically broadcasts to the agents. If the
control laws are identity-free, then the supervisor does not need to expend
computational power to distinguish between individual agents. Another
advantage of mean-field control approaches, from a theoretical point of view,
is that as the number of agents tends to infinity, the mean-field behavior of
the swarm is governed by a deterministic differential equation or difference
equation model, even though each agent might exhibit stochasticity in its
dynamics. Such models are more analytically tractable than models describing
the dynamics of a large number of individual agents.

In this article, we consider a mean-field stabilization problem motivated by
coverage problems in multi-agent control. A classical approach to multi-agent
coverage is described in [17], which presents a distributed method for
implementing Lloyd’s algorithm for positioning multiple agents in a domain
according to a given probability density function. The stochastic task
allocation problem considered in [8], where the goal is to stabilize a swarm
of agents evolving according to a continuous-time Markov chain to a target
distribution among a set of states (e.g., physical locations), can be viewed
as a mean-field version of this coverage problem. Stochastic task allocation
approaches have also been developed for swarms of agents that evolve according
to discrete-time Markov chains [1] and that follow control laws which depend
on the local density of agents [44, 22]. A drawback of Markov chain-based
approaches to this problem is that the state space of the agents needs to be
discretized beforehand.

In [48, 47], the authors consider the problem of stabilizing a swarm of agents
with general controllable dynamics to a target distribution. This approach has
been extended to the case of holonomic agents on bounded domains [19, 24] and
compact manifolds without boundary [20]. In these extensions, the control is
either a diffusion coefficient or a velocity field in the Fokker-Planck
partial differential equation (PDE) [53] that determines the spatio-temporal
evolution of the probability density of the agents, each of which is governed
by a reflected stochastic differential equation [57]. In [10], similar
stochastic control laws are constructed for agents evolving according to a
discrete-time deterministic nonlinear control system on a bounded Euclidean
domain.

An advantage of the stochastic coverage approaches in [8, 1, 48, 47, 19, 24,
20, 10] over the classical coverage strategy in [17] is that they only require
each agent’s control action to depend on its own state. The stochasticity
designed into the system ensures that the swarm reaches a target probability
density. A disadvantage of these approaches is that agents do not stop
switching between states even after the swarm has reached the target
probability density, resulting in an unnecessary expenditure of energy.
Moreover, a large number of agents is required for the swarm density to
stabilize close to the target density. One way to resolve these issues is to
design control laws that are functions of the local swarm density. Such
density-dependent or mean-field feedback laws have been proposed for agents
that evolve according to Markov chains on discrete state spaces [44, 22], for
agents that evolve according to ordinary differential equations on Euclidean
domains [26, 37], and for agents that evolve according to stochastic or
ordinary differential equations on compact manifolds without boundary [20].
The works [26, 37, 20] assume that the agents are holonomic, and  in order to
achieve global asymptotic stability, the target distribution is required to be
strictly positive everywhere on the domain.

In the context of these previous works, the main contributions of this article
are the following:

  1. 1.

Extension of the stochastic multi-agent coverage approach developed by the
authors in [19, 24] to non-interacting agents with nonholonomic dynamics that
evolve on a bounded subset of $\mathbb{R}^{d}$ or a compact manifold without
boundary, given a target probability density for the swarm that is bounded
from below by a positive constant.

  2. 2.

Development of a stochastic multi-agent coverage approach for interacting
agents governed by the hybrid switching diffusion model introduced by the
authors in [23]. In this approach, a mean-field feedback law (i.e., a control
law that depends on the local swarm density) is constructed to  globally
asymptotically stabilize the swarm to any given target probability density,
which is not necessarily positive everywhere.

Contribution 1 is partially motivated by the fact that extension of multi-
agent control strategies designed for Euclidean state spaces to general
manifolds [56, 9] is important, given that many mechanical systems are
naturally modeled on manifolds [15]. While the work [47], as in Contribution
1, does consider agents with general nonholonomic dynamics, due to the
assumption that the domain is unbounded, this work requires assumptions on the
behavior of the target probability density at infinity. The extension of the
coverage strategy presented in [19, 24] to the case of nonholonomic agents
evolving on manifolds is complicated by the fact that the associated PDEs are
not elliptic. Contribution 2 improves over existing work on mean-field
feedback laws [26, 37], since our approach does not make strong assumptions on
the regularity of solutions of the associated PDEs. Instead, we prove all the
regularity required to enable the stability analysis, which makes the analysis
much more technically involved. Moreover, we are able to stabilize a larger
class of probability densities than those considered in [19, 24, 20, 26, 37],
which require that the target probability density is strictly bounded from
below by a positive constant everywhere on the domain. A control law similar
to the one described in Contribution 2 was constructed by the authors in [11]
for a discrete-time control system. However, it is assumed in [11] that the
system is locally controllable within one time-step, and hence is fully
actuated. Moreover, in contrast to the mean-field feedback laws constructed in
[26, 37, 20] in which the control input is the agents’ velocity field, the
control inputs that we design in Contribution 2 are the transition rates of
the hybrid switching diffusion process that describes the agents’ dynamics.

This article is organized as follows. In Section II, we establish notation and
provide some definitions that are used throughout the article. In Section III,
we present and analyze the properties of the degenerate PDEs that describe the
mean-field model in the case where the agents do not interact with one
another. In Section IV, we present a semilinear PDE mean-field model for
stabilizing the density of a swarm of interacting agents and establish global
asymptotic stability properties of the model. In Section V, we validate the
control strategies presented in Sections III and IV with numerical
simulations.

##  II NOTATION

We denote the $n$-dimensional Euclidean space by $\mathbb{R}^{n}$.
$\mathbb{R}^{n\times m}$ refers to the space of $n\times m$ matrices, and
$\mathbb{R}_{+}$ refers to the set of non-negative real numbers. Given a
vector $\mathbf{x}\in\mathbb{R}^{n}$, $x_{i}$ denotes the $i^{th}$ coordinate
value of $\mathbf{x}$. For a matrix $\mathbf{A}\in\mathbb{R}^{n\times m}$,
$A_{ij}$ refers to the element in the $i^{th}$ row and $j^{th}$ column of
$\mathbf{A}$. For a subset $B\subset\mathbb{R}^{M}$, ${\rm int}(B)$ refers to
the interior of the set $B$. $\mathbb{C}$, $\mathbb{C}_{-}$, and
$\bar{\mathbb{C}}_{-}$ denote the set of complex numbers, the set of complex
numbers with negative real parts, and the set of complex numbers with non-
positive real parts, respectively. $\mathbb{Z}_{+}$ refers to the set of
positive integers.

We denote by $\Omega$ an open, bounded, and connected subset of an
$N$-dimensional smooth Riemannian manifold $M$ [39, 40] with a Riemannian
volume measure $d\mathbf{x}$. The boundary of $\Omega$ is denoted by
$\partial\Omega$. We denote by $\int_{\Omega}f(\mathbf{x})d\mathbf{x}$ the
integral of a function $f:\Omega\rightarrow\mathbb{R}$ with respect to the
Riemannian volume.

For each $1\leq p<\infty$, we define $L^{p}(\Omega)$ as the Banach space of
complex-valued measurable functions over the set $\Omega$ whose absolute value
raised to $p^{{\rm th}}$ power has finite integral.  The space $L^{1}$ is
equipped with the norm $\|f\|_{1}:=\int_{\Omega}|f(x)|dx$. We define
$L^{\infty}(\Omega)$ as the space of essentially bounded measurable functions
on $\Omega$. The space $L^{\infty}(\Omega)$ is equipped with the norm
$\|z\|_{\infty}=~{}$ess sup${}_{\mathbf{x}\in\Omega}|z(\mathbf{x})|,$ where
ess sup${}_{\mathbf{x}\in\Omega}(\cdot)$ denotes the essential supremum
attained by its argument over the domain $\Omega$.

For a given real-valued function $a\in L^{\infty}(\Omega)$,
$L^{2}_{a}(\Omega)$ refers to the set of all functions $f$ such that the norm
of $f$ is defined as
$\|f\|_{a}:=(\int_{0}^{1}|f(\mathbf{x})|^{2}a(\mathbf{x})d\mathbf{x})^{{1/2}}<\infty.$
We will always assume that the associated function $a$ is uniformly bounded
from below by a positive constant, in which case the space $L^{2}_{a}(\Omega)$
is a Hilbert space with respect to the weighted inner product
$\langle\cdot,\cdot\rangle_{a}:L^{2}_{a}(\Omega)\times
L^{2}_{a}(\Omega)\rightarrow\mathbb{R}$, given by $\langle
f,g\rangle_{a}=\int_{\Omega}f(\mathbf{x})\bar{g}(\mathbf{x})a(\mathbf{x})d\mathbf{x}$
for each $f,g\in L^{2}_{a}(\Omega)$, where $\bar{g}$ is the complex conjugate
of the function $g.$ When $a=\mathbf{1}$, where $\mathbf{1}$ is the function
that takes the value $1$ almost everywhere on $\Omega$, the space
$L^{2}_{a}(\Omega)$ coincides with the space $L^{2}(\Omega)$. For a function
$f\in L^{2}(\Omega)$ and a given constant $c$, we write $f\geq c$ to imply
that $f$ is real-valued and $f(\mathbf{x})\geq c$ for almost every (a.e.)
$\mathbf{x}\in\Omega$.

Suppose $e^{Xt}$ is the flow generated by a vector field $X$. Then $X$ defines
a differential operator on the set of smooth functions $C^{\infty}(M)$ through
the following action,

$(Xf)(\mathbf{x})=\lim_{t\rightarrow
0}\frac{f(e^{tX}(\mathbf{x}))-f(\mathbf{x})}{t}$ (1)

for all $\mathbf{x}\in\Omega$ and all $f\in C^{\infty}(\Omega)$.  This is the
differential geometric definition from [39] of a vector field $X$ as an
associated differential operator acting on the space of smooth functions.

Let $\mathcal{V}=\\{X_{1},...,X_{m}\\}$, $m\leq N$, be a collection of smooth
vector fields $X_{i}$, each defined as in (1). Let $[X,Y]$ denote the Lie
bracket of the vector fields $X$ and $Y$.  We define
$\mathcal{V}^{0}=\mathcal{V}$. For each $i\in\mathbb{Z}_{+}$, we define in an
iterative manner the set of vector fields
$\mathcal{V}^{i}=\\{[X,Y];~{}X\in\mathcal{V},~{}Y\in\mathcal{V}^{j-1},~{}j=1,...,i\\}$.
We will assume that the collection of vector fields $\mathcal{V}$ satisfies
the Chow-Rashevsky condition [4] (also known as Hörmander’s condition [13]),
i.e., the Lie algebra generated by the vector fields $\mathcal{V}$, given by
$\cup_{i=0}^{r}\mathcal{V}^{i}$, has rank $N$, for sufficiently large $r$.

A horizontal curve $\boldsymbol{\gamma}:[0,1]\rightarrow\Omega$ connecting two
points $\mathbf{x},\mathbf{y}\in M$ is a Lipschitz curve in $\Omega$ for which
there exist essentially bounded functions $a_{i}(t)$ such that

$\dot{\boldsymbol{\gamma}}(t)=\sum_{i=1}^{m}a_{i}(t)X_{i}(\boldsymbol{\gamma}(t))$
(2)

for almost every $t\in[0,1]$,  where $X_{i}\in\mathcal{V}$,
$\boldsymbol{\gamma}(0)=\mathbf{x}$, and $\boldsymbol{\gamma}(1)=\mathbf{y}$.
Then $\mathcal{V}$ defines a distance $d:\Omega\rightarrow\mathbb{R}_{\geq 0}$
on $M$ as

$\displaystyle d(\mathbf{x},\mathbf{y})={\rm
inf}~{}\\{\int_{0}^{1}|\dot{\boldsymbol{\gamma}}(t)|dt;$ $\boldsymbol{\gamma}$
is a horizontal curve $\displaystyle\text{ connecting $\mathbf{x}$ and
$\mathbf{y}$}\\}$

######  Definition II.1.

The domain $\Omega\subset M$ is said to be $\epsilon-\delta$ if there exist
$\delta>0$, $0<\epsilon\leq 1$ such that for any pair of points
$\mathbf{p},\mathbf{q}\in\Omega$, if $d(\mathbf{p},\mathbf{q})\geq\delta$,
then there exists a continuous curve
$\boldsymbol{\gamma}:[0,T]\rightarrow\Omega$ such that
$\boldsymbol{\gamma}(0)={\color[rgb]{0,0,0}\mathbf{p}}$,
$\boldsymbol{\gamma}(T)={\color[rgb]{0,0,0}\mathbf{q}}$, and

$\displaystyle\int_{0}^{1}|\dot{\boldsymbol{\gamma}}(t)|dt~{}\geq~{}\frac{1}{\epsilon}d({\color[rgb]{0,0,0}\mathbf{p},\mathbf{q}})$
$\displaystyle
d(\mathbf{{\color[rgb]{0,0,0}z}},\partial\Omega)~{}\geq~{}\epsilon~{}{\rm
min}(d({\color[rgb]{0,0,0}\mathbf{p},\mathbf{z}}),d({\color[rgb]{0,0,0}\mathbf{z},\mathbf{q}}))~{}~{}\forall{\color[rgb]{0,0,0}\mathbf{z}}\in\\{\boldsymbol{\gamma}(t):t\in[0,T]\\}$

The metric $d$ on $\Omega$ is known as the sub-Riemannian or Carnot-
Caratheodory metric [2, 13]. The topology induced by this metric on $d$
coincides with the usual Riemannian metric. We will assume that the radius
$r(\Omega)$ of $\Omega$, given by $r(\Omega)={\rm
sup}\\{d(\mathbf{x},\mathbf{y});\mathbf{x},\mathbf{y}\in M\\}$, is finite.

Given $a\in L^{\infty}(\Omega)$, with $a\geq c$ for a positive parameter
$c>0$, we define the weighted horizontal Sobolev space
$WH^{1}_{a}(\Omega)=\big{\\{}f\in L^{2}(\Omega):X_{i}(af)\in
L^{2}(\Omega)\text{ for }1\leq i\leq m\big{\\}}$. We equip this space with the
weighted horizontal Sobolev norm $\|\cdot\|_{WH^{1}_{a}}$, given by
$\|f\|_{WH^{1}_{a}}=\Big{(}\|f\|^{2}_{2}+\sum_{i=1}^{n}\|X_{i}(af)\|^{2}_{2}\Big{)}^{1/2}$
for each $f\in WH^{1}_{a}(\Omega)$. Here, the derivative action of $X_{i}$ on
a function $f$ is to be understood in the distributional sense. When
$a=\mathbf{1}$, where $\mathbf{1}$ is the constant function that is equal to
$1$ everywhere, we will denote $WH_{a}^{1}$ by $WH^{1}$.

Let $X$ be a Hilbert space with the norm $\|\cdot\|_{X}$. The space
$C([0,T];X)$ consists of all continuous functions $u:[0,T]\rightarrow X$ for
which $\|u\|_{C([0,T];X)}:=\max_{0\leq t\leq T}\|u(t)\|_{X}~{}<~{}\infty$. If
$Y$ is a Hilbert space, then $\mathcal{L}(X,Y)$ will denote the space of
linear bounded operators from $X$ to $Y$. We will also use the multiplication
operator $\mathcal{M}_{a}:L^{2}(\Omega)\rightarrow L^{2}(\Omega)$, defined as
$(\mathcal{M}_{a}u)(\mathbf{x})=a(\mathbf{x})u(\mathbf{x})$ for a.e.
$\mathbf{x}\in\Omega$ and each $u\in L^{2}(\Omega)$.

We will need an appropriate notion of a solution of the PDEs considered in
this paper. Toward this end, let $A$ be a closed linear operator that is
densely defined on a subset $\mathcal{D}(A)$, the domain of the operator, of a
Hilbert space $H$. We will define ${\rm spec}(A)$ as the set
$\\{\lambda\in\mathbb{C}:\lambda\mathbb{I}-A$ is not invertible in $X\\}$,
where $\mathbb{I}$ is the identity map on $X$. If $A$ is a bounded operator,
then $\|A\|_{op}$ will denote the operator norm induced by the norm defined on
$H$. From [25], we have the following definition.

######  Definition II.2.

For a given time $T>0$, a mild solution of the ODE

$\dot{u}(t)=Au(t);\hskip 5.69054ptu(0)=u_{0}\in H$ (3)

is a function $u\in C([0,T];X)$ such that $u(t)=u_{0}+A\int_{0}^{t}u(s)ds$ for
each $t\in[0,T]$.

Under appropriate conditions satisfied by $A$, the mild solution of a PDE is
given by a strongly continuous semigroup of linear operators,
$(\mathcal{T}(t))_{t\geq 0}$, that are generated by the operator $A$ [25].

######  Definition II.3.

A strongly continuous semigroup of linear operators $(\mathcal{T}(t))_{t\geq
0}$ on a Hilbert space $X$ is called positive if $u\in X$ such that $u\geq 0$
implies that $\mathcal{T}(t)u\geq 0$ for all $t\geq 0$.

##  III Stabilization without Agent Interactions

Given the definitions in Section II, consider the following reflected
stochastic differential equation (SDE) [52] constrained to a domain
$\Omega\subseteq M$:

$\displaystyle d\mathbf{Z}(t)$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{m}u_{i}(\mathbf{Z}(t))X_{i}dt+\sqrt{2}\sum_{i=1}^{m}v_{i}(\mathbf{Z}(t))X_{i}\circ
dW_{i}$ $\displaystyle+~{}\mathbf{n}(\mathbf{Z}(t))d\psi(t),$
$\displaystyle\mathbf{Z}(0)$ $\displaystyle=$
$\displaystyle~{}\mathbf{Z}_{0},$ (4)

where $\psi(t)\in\mathbb{R}$ is called the reflecting function or local time
[52], a stochastic process that constrains $\mathbf{Z}(t)$ to the domain
$\Omega$; $\mathbf{n}(\mathbf{x})$ is the normal to the boundary at
$\mathbf{x}\in\partial\Omega$; $W_{i}$ are $m$ copies of the one-dimensional
Wiener process; and $u_{i}$ and $v_{i}$ are $m$ feedback laws. In the above
SDE (III), the notation $\circ$ is used to mean that the SDE should be
interpreted in the sense of Stratonovich [34]. Let $y(\mathbf{x},t)$ denote
the probability density of the random variable $\mathbf{Z}(t)$, defined as
$\mathbb{P}(\mathbf{Z}(t)\in A)=\int_{A}y(\mathbf{x},t)d\mathbf{x}$. In this
section, we consider the following control problem:

######  Problem III.1.

Given a target probability density $f$ on $\Omega$, design control laws
$u_{i}(\mathbf{x})$ and $v_{i}(\mathbf{x})$ in Eq. (III) such that the
probability density $y(\mathbf{x},t)$of the stochastic process
$\mathbf{Z}(t)$, which evolves according to Eq. (III), converges
asymptotically to $f$.

The motivation for this problem comes from stochastic coverage applications in
swarm robotics that are framed as follows. Let the random variable
$\mathbf{Z}_{j}(t)$, $j\in\\{1,...,N_{p}\\}$, denote the position of the
$j^{th}$ robot in a swarm of $N_{p}$ robots at time $t$. This position evolves
according to Eq. (III), in which $u_{i}$ and $v_{i}$ are control laws that
govern each robot’s motion.  Since each robot follows the same control laws
$u_{i}$ and $v_{i}$, the random variables $\mathbf{Z}_{i}(t)$ are independent
and identically distributed. Then, denoting by $\delta_{\mathbf{x}}$ the delta
distribution at $\mathbf{x}\in M$, the empirical distribution
$\frac{1}{N_{p}}\sum_{j=1}^{N_{p}}\delta_{\mathbf{Z}_{j}(t)}$, which
represents the distribution of the robots in space, converges to the density
$y(\mathbf{x},t)$ as $N_{p}\rightarrow\infty$ due to the law of large numbers.

The stabilization problem III.1 has been considered by the authors in [19] for
the case where the system is holonomic and the vector fields
$X_{i}=\frac{\partial}{\partial x_{i}}$ are the standard coordinate vector
fields. The goal in this section is to extend the results in [19] to the
general case where the number of vector fields $X_{i}$ is possibly less than
the dimension $N$ of the state space $M$. Such density stabilization problems
were first considered for the case where the domain $\Omega$ is the whole of
the Euclidean space $\mathbb{R}^{n}$ in [48, 47]. When time is discrete, and
the system is controllable in one time step, this problem has been considered
in [10].

The main difficulty in extending the results from [19] is that when the number
of control vector fields, $m$, is less than the dimension of the state space
$M$, the generator $\sum_{i=1}^{m}(v_{i}X_{i})^{2}+u_{i}X_{i}$ of the
stochastic process $\mathbf{Z}(t)$ is not elliptic, which makes standard
results in the literature on parabolic PDEs inapplicable. In particular, let
$A=\sum_{i=1}^{m}(v_{i}X_{i})^{2}$. The associated probability density
$y(\mathbf{x},t)$ of the process $\mathbf{Z}(t)$ evolves according to the PDE

$\displaystyle
y_{t}=A^{*}y-\nabla\cdot(\sum_{i=1}^{m}{u}_{i}(\mathbf{x})X_{i}y)$
$\displaystyle~{}~{}in~{}~{}\Omega\times[0,T]$ (5) $\displaystyle
y(\cdot,0)=y^{0}$ $\displaystyle~{}~{}in~{}~{}\Omega$

with zero flux boundary conditions, where $\nabla\cdot$ denotes the divergence
operation with respect to the measure $d\mathbf{x}$, and $A^{*}$ is the
adjoint of the operator $A$. The stabilization problem III.1 is an open-loop
control problem for the PDE (5), in which the goal is to stabilize the
solution $y(\mathbf{x},t)$ of (5) to a target function $f$.

The operator $A$ is not elliptic in general, but only hypoelliptic.
Particularly, if $f\in C^{\infty}_{0}(\Omega)$ has compact support $K$, then,
due to the Chow-Rashevsky’s Lie rank condition, if $u$ is a function on
$\Omega$ such that $Au=f$, then $u$ is smooth on $K$ [13]. Using this property
of $A$, we will extend the stabilization results of [19, 24] to the case where
the agents have nonholonomic dynamics.

First, formally, we will provide a number of candidate control laws that are
solutions to Problem III.1. Given these control laws, in Section III-A we will
present a stability analysis of a class of PDEs that coincide with (5) when
the operators $X_{i}$ are formally skew-adjoint with respect to the volume
form $d\mathbf{x}$, that is, $X_{i}^{*}=-X_{i}$.

Suppose that the operators $X_{i}$ are formally skew-adjoint. Let $f\in
W^{1,\infty}(\Omega)$ be a positive function that is bounded from below by a
positive number and for which $\int_{\Omega}f(\mathbf{x})d\mathbf{x}=1$. If we
set $u_{i}(\cdot)=X_{i}g/g$ and $v_{i}(\cdot)=1$ for each $i\in\\{1,...,m\\}$
and all $t\geq 0$, then the PDE (5) becomes

$\displaystyle
y_{t}=\sum_{i=1}^{m}X_{i}^{2}y-\nabla\cdot(\sum_{i=1}^{m}\frac{X_{i}f}{f}X_{i}y)~{}$
$\displaystyle~{}in~{}~{}\Omega\times[0,T]$ (6) $\displaystyle
y(\cdot,0)=y^{0}$ $\displaystyle~{}in~{}~{}\Omega$

Let $\nabla_{H}$ be the horizontal gradient operator, which maps functions to
vector fields and is defined as $\nabla_{H}g=\sum_{i=1}^{m}(X_{i}g)X_{i}$ When
the manifold $M$ is a Lie group $G$ and it is unimodular, i.e., the left- and
right-Haar measures [39] coincide, then we have that
$\nabla\cdot\nabla_{H}(\cdot)=\sum_{i=1}^{m}X_{i}^{2}$ [3]. Hence, if we set
$y=f$, then

$\sum_{i=1}^{m}X_{i}^{2}y-\nabla\cdot\left(\sum_{i=1}^{m}\frac{X_{i}f}{f}X_{i}y\right)=\sum_{i=1}^{m}X_{i}^{2}g-\nabla\cdot(\nabla_{H}g)=0$
(7)

Thus, $f$ is an equilibrium solution of the PDE (6). Our goal in Section III-A
will be to show that $f$ is the globally exponentially stable equilibrium
solution of PDE (6) on the the set of square-integrable probability densities.

Let $a=\frac{1}{f}$. Then the operator
$\sum_{i=1}^{m}X_{i}^{2}-\nabla\cdot(\sum_{i=1}^{m}\frac{X_{i}f}{f}X_{i})$ can
be alternatively expressed as
$\nabla\cdot(\frac{1}{a(\mathbf{x})}\nabla_{H}(a(\mathbf{x})\cdot))$.
Similarly, one can also consider the feedback laws $u_{i}=0$ and
$v_{i}=\frac{1}{f}$. In this case, the corresponding operator of interest is
given by

$A^{*}=\nabla\cdot({a(\mathbf{x})}\nabla_{H}(a(\mathbf{x})\cdot))$ (8)

This control law is similar to the one presented in [19] for the case of
holonomic agents, where instead of the Stratonovich integral, we considered
the Ito integral, and the resulting generator was of the form
$\nabla\cdot(\nabla_{H}(a(\mathbf{x})\cdot))$.

###  III-A Stability Analysis

The preceding discussion motivates us to study stability properties of PDEs
associated with a class of hypoelliptic operators that have a given
probability density as their equilibrium solution. In this section, we will
provide a semigroup theoretic analysis of a class of such PDEs. There have
been a number of works on semigroups generated by hypoelliptic operators on
manifolds without boundary [33], or manifolds with boundary under the
Dirichlet boundary [60]. Due to the term $a(\mathbf{x})$, the operators that
we consider are more general than those in [60]. There has also been work on
long-term behavior of hypoelliptic diffusions to uniform distributions for the
special case of Carnot groups [54] and to more general equilibrium
distributions, as well as on more general state spaces, for certain examples
of hypoelliptic diffusions using log-Sobolev inequalities [6]. In comparison,
our results hold for a more general class of degenerate diffusions by
establishing a spectral gap for the generator, which can be guaranteed to
exist for equilibrium distributions that are bounded uniformly from above and
below by positive numbers.

Before we present our stability analysis, we give some more preliminary
definitions. Given $a,b\in L^{\infty}(\Omega)$ such that $a\geq c$ and $b\geq
c$ for some positive constant $c$, and
$\mathcal{D}(\omega^{b}_{a})=WH^{1}_{a}(\Omega)$, we define the sesquilinear
form
$\omega^{b}_{a}:\mathcal{D}(\omega^{b}_{a})\times\mathcal{D}(\omega^{b}_{a})\rightarrow\mathbb{C}$
as

$\omega^{b}_{a}(u,v)=\sum_{i=1}^{m}\int_{\Omega}b(\mathbf{x})X_{i}(a(\mathbf{x})u(\mathbf{x}))\cdot
X_{i}(a(\mathbf{x})\bar{v}(\mathbf{x}))d\mathbf{x}$ (9)

for each $u\in\mathcal{D}(\omega^{b}_{a})$. We associate with the form
$\omega^{b}_{a}$ an operator $A_{a}^{b}:\mathcal{D}(A_{a}^{b})\rightarrow
L^{2}_{a}(\Omega)$, defined as $A_{a}^{b}u=v$, if
$\omega_{a}^{b}(u,\phi)=\langle v,\phi\rangle_{a}$ for all
$\phi\in\mathcal{D}(\omega^{b}_{a})$ and for all
$u\in\mathcal{D}(A_{a}^{b})=\\{g\in\mathcal{D}(\omega_{a}^{b}):~{}\exists f\in
L^{2}_{a}(\Omega)~{}\text{s.t.}~{}\omega_{a}^{b}(g,\phi)=\langle
f,\phi\rangle_{a}~{}\forall\phi\in\mathcal{D}(\omega^{b}_{a})\\}$. When the
$X_{i}$ are formally skew-adjoint, the operator $A_{a}^{b}$ is a weak
formulation of the the second-order partial differential operator
$\sum_{i=1}^{m}X^{*}_{i}(b(\mathbf{x})X_{i}(a(\mathbf{x})~{}\cdot~{}))$. An
advantage of using this weak formulation of the operator, rather than a strong
formulation, is that one does not need to establish that the domain
$\mathcal{D}(A_{a}^{b})$ contains twice (weakly) differentiable functions,
which might not be true in general, given the very weak regularity assumed on
the boundary of the domain $\Omega$ [32]. Formally, we will be studying
properties of the following PDE,

$\displaystyle
y_{t}=\sum_{i=1}^{m}X_{i}^{*}(b(\mathbf{x})X_{i}(a(\mathbf{x})y))~{}~{}in~{}~{}\Omega\times[0,T],$
(10)

Note that for $b=1/a$, $b=a$, and $b=\mathbf{1}$, we recover the operators
introduced in the previous subsection.

The proofs of the results in this section follow closely, almost verbatim, to
those for the elliptic case considered by the authors in [24]. Therefore, due
to space limitations, we only include the proof of the stability result, which
is stated in Theorem III.6. These results will be used extensively in Section
IV, where we consider the case in which the agents have local interactions
with one another. The main technical difference between the proofs in [24] and
the proofs of the results presented in this section is that here, we use the
horizontal Sobolev spaces $WH^{1}({\Omega})$ to establish semigroup generation
properties of the generator $A$, instead of the classical Sobolev space
$H^{1}(\Omega)$. Due to the bracket generating property of the vector fields
$\mathcal{V}$, it is known that the space $WH^{1}({\Omega})$ has many
properties similar to the classical Sobolev space $H^{1}(\Omega)$ [29, 49].

We will need the following assumption for the results in this section to hold
true.

######  Assumption III.2.

Only one of the following conditions holds:

  1. 1.

The domain $\Omega$ is $\epsilon-\delta$ and $M=\mathbb{R}^{n}$.

  2. 2.

The manifold $M$ is compact and without a boundary.

  3. 3.

The boundary $\partial\Omega$ of the domain $\Omega$ is $C^{1}$ and  ${\rm
span}\\{\mathcal{V}\\}=T_{\mathbf{x}}M$ for all $\mathbf{x}\in\Omega$.

######  Lemma III.3.

  1. 1.

The operator $A_{a}^{b}:\mathcal{D}(A_{a}^{b})\rightarrow L^{2}_{a}(\Omega)$
is closed, densely defined, and self-adjoint.

  2. 2.

The operator $A_{a}^{b}:\mathcal{D}(A_{a}^{b})\rightarrow L^{2}_{a}(\Omega)$
has a purely discrete spectrum.

  3. 3.

$-A^{b}_{a}$ generates a semigroup of operators
$(\mathcal{T}^{b}_{a}(t))_{t\geq 0}$. The semigroup
$(\mathcal{T}^{b}_{a}(t))_{t\geq 0}$ is positive and a contraction.

That the operator $A_{a}^{b}$ is densely defined follows from the fact that
the domain of the form $\omega_{a}^{b}$, which is $WH_{a}^{1}(\Omega)$, is
dense in $L^{2}_{a}(\Omega)$ when equipped with the norm
$(\omega_{a}^{b}(u,u)+\|u\|^{2}_{2})^{1/2}$. The result on discreteness of the
spectrum of the operator $A^{b}_{a}$ follows from the compactness of the
embedding $WH^{1}(\Omega)$ in $L^{2}(\Omega)$. The positivity properties of
the semigroup $\mathcal{T}^{b}_{a}(t)$, we need the fact that if $f\in
WH^{1}(\Omega)$, then $|f|\in WH^{1}(\Omega)$. This result is known for the
case where $\Omega$ is a subset of $\mathbb{R}^{n}$. The proof for the general
case where $\Omega$ is a manifold follows verbatim the results of [29, 49],
since the proof only requires the density of the space $C^{\infty}(M)$ in
$WH^{1}(\Omega)$, which can be verified.  Using these facts, the proof of the
previous lemma follows the proof of [24][Lemma IV.1] very closely. The
contractivity of the semigroup follows from [50][Proposition 1.51] due to the
fact that the form $\omega_{a}^{b}$ is accretive; that is,
$\omega_{b}^{a}(u,u)\geq 0$ for all $u\in\mathcal{D}(\omega_{b}^{a})$.

######  Proposition III.4.

[29, 49] Given $f\in WH^{1}(\Omega)$, we have that $|f|\in WH^{1}(\Omega)$.

Using the above proposition and the established properties of the operator
$A_{a}^{b}$, we can prove the following results on the semigroup generated by
the operator $-A^{b}_{a}$ concerning its positivity and mass-conserving
properties. These results will be used in the analysis presented in Section
IV.

######  Corollary III.5.

The operator $-A^{b}_{a}$ generates a semigroup of operators
$(\mathcal{T}^{b}_{a}(t))_{t\geq 0}$ acting on $L^{2}_{a}({\Omega})$. The
semigroup $(\mathcal{T}^{b}_{a}(t))_{t\geq 0}$ is positive. Furthermore, if
$a=b=\mathbf{1}$ is the constant function equal to $1$ almost everywhere, then
$\|y^{0}\|_{\infty}\leq 1$ implies that
$\|\mathcal{T}^{b}_{a}(t)y^{0}\|_{\infty}\leq 1$ for all $t\geq 0$.

Additionally, this semigroup has the following mass conservation property: if
$y_{0}\in L^{2}_{a}(\Omega)$, $y^{0}\geq 0$ and
$\int_{\Omega}y^{0}(\mathbf{x})d\mathbf{x}=1$, then
$\int_{\Omega}(\mathcal{T}_{a}^{b}(t)y^{0})(\mathbf{x})d\mathbf{x}=\int_{\Omega}(\mathcal{T}^{B}_{a}(t)y^{0})(\mathbf{x})d\mathbf{x}=1$
for all $t\geq 0$.

Lastly, we establish the following important result on the long-term stability
properties of the semigroups associated with the operators in
$\sum_{i=1}^{m}X^{*}_{i}(b(\mathbf{x})X_{i}(a(\mathbf{x})~{}\cdot~{}))$.

######  Theorem III.6.

(Exponential stability of semigroup)
The semigroup $(\mathcal{T}_{a}^{b}(t))_{t\geq 0}$ generated by the operator
$-A_{a}^{b}$ is analytic. Moreover, $0$ is a simple eigenvalue of $-A_{a}^{b}$
corresponding to the eigenvector $f=1/a$. Hence, if $y^{0}\geq 0$ and
$\int_{\Omega}y^{0}(\mathbf{x})d\mathbf{x}=\int_{\Omega}f(\mathbf{x})d\mathbf{x}=1$,
then the following estimate holds for some positive constants $M_{0},\lambda$
and all $t\geq 0$:

$\displaystyle\|\mathcal{T}_{a}^{b}(t)y^{0}-f\|_{a}~{}$ $\displaystyle\leq$
$\displaystyle~{}M_{0}e^{-\lambda t}\|y^{0}-f\|_{a}$ (11)

###### Proof.

We only prove the stability result, since  the analyticity of the semigroup
and the conservation of mass properties can be established as in [24].  The
semigroup is compact due to the compactness of the embedding
$\mathcal{D}(A_{a}^{b})\subset WH^{1}(\Omega)$ in $L^{2}(\Omega)$. Since the
semigroup is analytic and compact, in order to establish  its stability
properties, it is sufficient to identify the eigenvectors associated with the
eigenvalue $0$. In the proof of the corresponding result in [24], we used the
Poincaré inequality to establish the uniqueness of the eigenvector of constant
functions, corresponding to the eigenvalue $0$ of the Laplacian $\Delta$. It
is not clear whether the Poincaré inequality holds for the operator
$-A^{b}_{a}$ for condition 2 in Assumption III.2. Hence, instead of using a
Poincaré inequality, we will prove that the kernel of the operator
$-A^{b}_{a}$ consists only of constant functions. Suppose $u\in\mathcal{D}(A)$
is such that $Au=\mathbf{0}$, where $A:=A^{\mathbf{1}}_{\mathbf{1}}$.  This
implies that $<Au,u>$ $=\sum_{i=1}^{m}(X_{i}u)^{2}=0$. Since the operator $A$
satisfies the Lie rank condition, from regularity results due to Hörmander
[13], we can infer that $u$ is locally smooth everywhere in $\Omega$. Then we
know that, for a given horizontal curve
$\boldsymbol{\gamma}:[0,1]\rightarrow\Omega$,
$u(\boldsymbol{\gamma}(1))=u(\boldsymbol{\gamma}(0))+\int_{0}^{1}\sum_{i=1}^{m}a_{i}(t)X_{i}u({\boldsymbol{\gamma}(t)})dt$
since $u(\boldsymbol{\gamma}(t))$ satisfies the differential equation
$\dot{u}(\boldsymbol{\gamma}(t))=\sum_{i=1}^{m}a_{i}(t)X_{i}u(\boldsymbol{\gamma}(t))$,
where $a_{i}(t)$ are the essentially bounded functions associated with the
curve $\boldsymbol{\gamma}(t)$ according to (2). Hence,
$u(\boldsymbol{\gamma}(1))-u(\boldsymbol{\gamma}(0))=\int_{0}^{1}\sum_{i=1}^{m}a_{i}(t)X_{i}u({\boldsymbol{\gamma}(t)})dt=0$
because $\sum_{i=1}^{m}(X_{i}u)^{2}=0$.  Note that we require the local
smoothness of $u$ in order to make sense of the term
$\int_{0}^{1}\sum_{i=1}^{m}a_{i}(t)Xu({\boldsymbol{\gamma}(t)})dt$. Since
$\mathcal{V}$ is bracket generating, we can choose $\boldsymbol{\gamma}(t)$
such that $\boldsymbol{\gamma}(0)$ and $\boldsymbol{\gamma}(1)$ are the given
initial and final conditions in $\Omega$. Hence, we have that $u$ is constant
everywhere on $\Omega$. This implies that $A\mathbf{1}=\mathbf{0}$, and hence
$A_{a}^{b}f=\mathbf{0}$ due to the assumption that $a,b$ are uniformly bounded
from below by a positive constant. ∎

##  IV Stabilization with Local Agent Interactions

In Section III, the probability densities that we stabilized were assumed to
be uniformly bounded from below by a positive number. Without this assumption,
the semigroups that were constructed would not be globally asymptotically
stable. In this section, we will introduce a semilinear PDE model for
stabilizing a swarm to probability densities that possibly have supports that
are disconnected.

As in Section III, $\Omega$ will denote an open bounded subset of a manifold,
and we consider a collection of vector fields
$\mathcal{V}=\\{X_{1},...,X_{m}\\}$ satisfying the Chow-Rashevsky condition.
Let $A:=A^{\mathbf{1}}_{\mathbf{1}}$ be the operator defined in Section III,
where $\mathbf{1}$ denotes the function that is equal to $1$ almost everywhere
on $\Omega$. We will also need the spaces
$\mathbf{L}^{2}(\Omega)=L^{2}(\Omega)\times L^{2}(\Omega)$ and
$\mathbf{L}^{\infty}(\Omega)=L^{\infty}(\Omega)\times L^{\infty}(\Omega)$ with
the standard norms inherited from the spaces $L^{2}(\Omega)$ and
$L^{\infty}(\Omega)$.

We will consider the following PDE model,

$\displaystyle(y_{1})_{t}=-Ay_{1}-q_{1}(\mathbf{x},t)y_{1}+q_{2}(\mathbf{x},t)y_{2}$
$\displaystyle~{}in~{}~{}\Omega\times[0,T]$ (12)
$\displaystyle(y_{2})_{t}=q_{1}(\mathbf{x},t)y_{1}-q_{2}(\mathbf{x},t)y_{2}$
$\displaystyle~{}in~{}~{}\Omega\times[0,T]$
$\displaystyle\mathbf{y}(\cdot,0)=\mathbf{y}^{0}$
$\displaystyle~{}in~{}~{}\Omega$ $\displaystyle\mathbf{n}\cdot\nabla y_{1}=0$
$\displaystyle~{}in~{}~{}\partial\Omega\times[0,T],$

where $y_{1}$ and $y_{2}$ are non-negative functions and $q_{i}$ are reaction
parameters. This PDE model is the forward equation of a hybrid switching
diffusion process (HSDP) [62]. In addition to a continuous spatial state
$\mathbf{Z}(t)$, each agent is associated with a discrete state
$Y(t)\in\\{0,1\\}$ at each time $t$. The hybrid switching diffusion process
$(\mathbf{Z}(t),Y(t))$ can be represented as a system of SDEs of the form

$\displaystyle d\mathbf{Z}(t)$ $\displaystyle=$
$\displaystyle\sqrt{2}(1-Y(t))\sum_{i=1}^{m}X_{i}\circ dW_{i}\
+\mathbf{n}(\mathbf{Z}(t))d\psi(t),$ $\displaystyle\mathbf{Z}(0)$
$\displaystyle=$ $\displaystyle~{}\mathbf{Z}_{0}.$ (13)

The PDE (12) is related to the SDE (13), for each $k\in\\{0,1\\}$, through the
relation
$\mathbb{P}(Y(t)=k,\mathbf{Z}(t)\in\Gamma)=\int_{\Gamma}y_{{\color[rgb]{0,0,0}k+1}}(\mathbf{x},t)d\mathbf{x}$
for all $t\in[0,T]$ and all measurable $\Gamma\subset\Omega$. The transitions
of the variable $Y(t)$ from one discrete state to another is determined by two
functions $q_{i}:\Omega\rightarrow[0,\infty]$ in the following way,

$\displaystyle\mathbb{P}(Y(t+h)=1|Y(t)=0)=q_{1}(\mathbf{Z}(t),t)h+o(h)$ (14)
$\displaystyle\mathbb{P}(Y(t+h)=0|Y(t)=1)=q_{2}(\mathbf{Z}(t),t)h+o(h)$ (15)

The state $Y(t)=0$ corresponds to the state in which agents diffuse in space
according to the reflected SDE, and the state $Y(t)=1$ corresponds to a state
in which they are motionless. Therefore, unlike the process considered in
Section III, each agent has two discrete states, between which it jumps
according to the transition rates $q_{i}(\mathbf{x},t)$ (also called reaction
parameters). We will treat the transition rates $q_{i}(\mathbf{x},t)$ as the
control inputs, instead of the velocity and diffusion parameters
$(u_{i},v_{i})$. Since we will allow the control inputs to be functions of the
density of the random variables $(\mathbf{Z}(t),Y(t))$,  i.e., the density of
agents in each state at time $t$, this reaction-based control mechanism
depends on interactions among agents that enable them to estimate these
densities, e.g., via local sensing, wireless communication, or physical
encounters. Due to the density-dependent transition rates, the forward
equation is a semilinear PDE.

We will consider the following problem in this section.

######  Problem IV.1.

Let $y^{d}\in L^{\infty}(\Omega)$ be a target probability density. Construct a
mean-field feedback law $K_{i}:\mathbf{L}^{2}(\Omega)\rightarrow
L^{\infty}(\Omega)$ such that if $u_{i}(\cdot,t)=K_{i}(\mathbf{y}(t))$ for all
$i\in\\{1,2\\}$ and all $t\geq 0$, then the system (12) is globally
asymptotically stable about the equilibrium
$\mathbf{y}^{d}=[\mathbf{0}~{}~{}y^{d}]^{T}$.

Before we address this problem, we make some additional assumptions on the
domain $\Omega$ and the operator $A$. Toward this end, we present the
following definitions.

######  Definition IV.2.

We will say that $\Omega$ is a $C^{1,1}$ domain if each point
$\mathbf{x}\in\partial\Omega$ has a neighborhood $\mathcal{N}$ such that
$\Omega\cap\mathcal{N}$ is represented by the inequality
$x_{n}<\gamma(x_{1},...,x_{n-1})$ in some Cartesian coordinate system for some
function $\gamma:\mathbb{R}^{n-1}\rightarrow\mathbb{R}$ that is at least once
differentiable and has derivatives of order $1$ that are Lipschitz continuous.

######  Definition IV.3.

The domain $\Omega$ will be said to satisfy the chain condition if there
exists a constant $C>0$ such that for every $\mathbf{x},\bar{\mathbf{x}}$
$\in$ $\Omega$ and every positive ${\color[rgb]{0,0,0}j}\in\mathbb{Z}_{+}$,
there exists a sequence of points $\mathbf{x}_{i}$ $\in$ $\Omega$, $0\leq
i\leq{\color[rgb]{0,0,0}j}$, such that $\mathbf{x}_{0}=\mathbf{x}$,
$\mathbf{x}_{\color[rgb]{0,0,0}j}=\bar{\mathbf{x}}$, and
$|\mathbf{x}_{i}-\mathbf{x}_{i+1}|\leq\frac{C}{{\color[rgb]{0,0,0}j}}|\mathbf{x}-\bar{\mathbf{x}}|$
for all $i=0,...,{\color[rgb]{0,0,0}j}-1$. Here $|\cdot|$ denotes the standard
Euclidean norm.

Note that every convex domain satisfies the chain condition.

In this section, we will make some stronger assumptions on the generator and
the domain $\Omega$ than those made in the previous section.

######  Assumption IV.4.

Only one of the following conditions holds:

  1. 1.

If $\Omega\neq M$, then $\Omega$ is a bounded subset of
${\color[rgb]{0,0,0}\mathbb{R}^{n}}$,

$-A=\sum_{i=1}^{{\color[rgb]{0,0,0}n}}\partial_{x_{i}}^{2}=\Delta$ is the
Laplacian, and $\Omega$ is a $C^{1,1}$ domain in the sense of Definition IV.2
and satisfies the chain condition in Definition IV.3.

  2. 2.

The set $\Omega$ is a compact manifold $M$ without a boundary  and
$-A=\sum_{i=1}^{m}X_{i}^{*}X_{i}$.

Given these assumptions, we have the following result due to Gaussian
estimates proved by [16] for the Laplacian $\Delta$, and by [33] for sub-
Laplacians. We will use this result in the subsequent analysis.

######  Theorem IV.5.

Let $(\mathcal{T}(t))_{t\geq 0}$ be the semigroup generated by the operator
$-A$. Let $y^{0}\in L^{2}(\Omega)$ be non-negative. Then there exists a
constant $C>0$ and time $T>0$, independent of $y^{0}$, such that
$\mathcal{T}(t)y^{0}\geq C\|y^{0}\|_{1}$ for all $t\geq T$.

In order to address Problem IV.1, we define the following maps
$F_{i}:L^{2}(\Omega)\rightarrow{\color[rgb]{0,0,0}L^{\infty}(\Omega)}$,
$i\in\\{1,2\\}$,

$(F_{i}(f))(\mathbf{x})=r_{i}(f(\mathbf{x})-y^{d}(\mathbf{x}))$ (16)

for almost every $\mathbf{x}\in\Omega$ and all $f\in L^{2}(\Omega)$, where
${\color[rgb]{0,0,0}r_{i}}:\mathbb{R}\rightarrow[0,c]$ are globally Lipschitz
functions for some positive number $\beta$, such that the functions $r_{1}$
and $r_{2}$ have supports equal to the intervals  $(-\infty,0]$ and
$[0,\infty)$, respectively. Our candidate mean-field feedback law $K_{i}$ for
addressing Problem IV.1 will be
$K_{i}(\mathbf{y})=F_{i}({\color[rgb]{0,0,0}y_{2}})$ for each $i\in\\{1,2\\}$.
Then the resulting closed-loop PDE is given by

$\displaystyle(y_{1})_{t}=-Ay_{1}-F_{1}(y_{2})y_{1}+F_{2}(y_{2})y_{2}$
$\displaystyle~{}~{}in~{}~{}\Omega\times[0,T]$ (17)
$\displaystyle(y_{2})_{t}=F_{1}(y_{2})y_{1}-F_{2}(y_{2})y_{2}$
$\displaystyle~{}~{}in~{}~{}\Omega\times[0,T]$
$\displaystyle\mathbf{y}(\cdot,0)=\mathbf{y}^{0}$
$\displaystyle~{}in~{}~{}\Omega$ $\displaystyle\mathbf{n}\cdot\nabla y_{1}=0$
$\displaystyle~{}in~{}~{}\partial\Omega\times[0,T],$

where the Neumann boundary condition in the last equation is specified only
for the case where the boundary $\partial\Omega$ is nonempty. Since the
transition rates are a functions of the distribution of the random variable,
the relation between the system of SDEs (13) and PDEs (17) is no longer
straightforward. For the choice of control law $F_{i}$, the SDE becomes a
stochastic process of Mckean-Vlasov type [46, 36], and further analysis is
required to establish a rigorous connection between the two systems (12) and
(17). Such an analysis is beyond the scope of this article, and is left for
future work.

Our main goal in this section will be to establish the asymptotic stability of
the PDE (17) given Assumption IV.4. Before we begin the stability analysis of
the above PDE model, we point out that standard approaches to stability
analysis, such as linearization-based approaches or Lyapunov functional
arguments, are not immediately applicable, as we demonstrate in the following
two remarks.

######  Remark IV.6.

(Lack of exponential stability) Consider the linearization of the PDE (17)
about the target equilibrium density
$\mathbf{y}^{d}=[\mathbf{0}~{}~{}y^{d}]^{T}$. It can be verified that the
(Fréchet) derivative of the nonlinear operators $F_{i}$ about $\mathbf{y}^{d}$
is the $\mathbf{0}$ operator. Therefore, the linearization of the PDE about
the equilibrium $\mathbf{y}^{d}$ is:

$\displaystyle(\tilde{y}_{1})_{t}=-A\tilde{y}_{1}$
$\displaystyle~{}~{}in~{}~{}\Omega\times[0,T]$
$\displaystyle(\tilde{y}_{2})_{t}=\mathbf{0}$
$\displaystyle~{}~{}in~{}~{}\Omega\times[0,T]$

Clearly, this PDE is not exponentially stable since the spectrum of its
generator, $\mathbf{A}\mathbf{y}=[Ay_{1}~{}\mathbf{0}]^{T}$, has an infinite
number of eigenvalues at $0$. Hence, the PDE (17) cannot be locally
exponentially stable about the equilibrium $\mathbf{y}^{d}$.

######  Remark IV.7.

(Difficulty in using LaSalle’s principle) Another standard approach to
establish asymptotic stability of dynamical systems is LaSalle’s invariance
principle [35]. However, the application of LaSalle’s invariance principle for
stability analysis of infinite-dimensional dynamical systems, such as the PDE
(17), requires that the trajectories of the system remain in a compact set for
all time. The compactness of trajectories for solutions of parabolic PDEs is
usually inferred from the regularizing effect of the diffusion component of
the dynamics. This is not straightforward to establish for solutions
$\mathbf{y}$ of the PDE (17) due to the fact that the diffusion operator $A$
acts only on the first state $y_{1}$, and therefore it cannot be guaranteed
that the state $y_{2}$ lies in a Sobolev space.

Due to the technical issues pointed out in Remarks IV.6 and IV.7, we will use
an alternative approach to establish asymptotic stability of the PDE (17)
based on the monotonicity properties of the PDE. In order to perform stability
analysis of the PDE (17), we will need a suitable notion of a solution. Toward
this end, we use the following definition.

######  Definition IV.8.

Let $(\mathcal{T}(t))_{t\geq 0}$ be the semigroup generated by the operator
$-A$. We will say that the PDE (17) has a local mild solution if there exist
$T>0$ and $\mathbf{y}\in C([0,T];\mathbf{L}^{2}({\Omega}))$ such that

$\displaystyle y_{1}(\cdot,t)$ $\displaystyle=$
$\displaystyle\mathcal{T}(t)y^{0}_{1}-\int_{0}^{t}\mathcal{T}(t-s)\Big{(}F_{1}(y_{2}(\cdot,s))y_{1}(\cdot,s)\Big{)}ds$
$\displaystyle+\int_{0}^{t}\mathcal{T}(t-s)\Big{(}F_{2}(y_{2}(\cdot,s))y_{2}(\cdot,s)\Big{)}ds,$
$\displaystyle y_{2}(\cdot,t)$ $\displaystyle=$ $\displaystyle
y^{0}_{2}+\int_{0}^{t}F_{1}(y_{2}(\cdot,s))y_{1}(\cdot,s)ds$ (18)
$\displaystyle-\int_{0}^{t}F_{2}(y_{2}(\cdot,s))y_{2}(\cdot,s)ds$

for all $t\in[0,T]$. We will say that the PDE (17) has a unique global
solution if it has a unique local mild solution for every $T>0$.

To establish the existence of solutions of the PDE (17), we will need the
operator
$\mathbf{A}:\mathcal{D}(\mathbf{A}):\rightarrow\mathbf{L}^{2}(\Omega)$,
defined as

$\mathbf{A}\mathbf{y}=\begin{bmatrix}Ay_{1}\\\ \mathbf{0}\end{bmatrix}$

for all $\mathbf{y}\in\mathcal{D}(\mathbf{A})=\mathcal{D}(A)\times
L^{2}(\Omega)$.

Our next goal will be to construct global solutions of the PDE (17). First, we
will show that the solutions of the PDE (17) remain essentially bounded if the
initial condition is essentially bounded. Toward this end, we first establish
this property for a related autonomous linear PDE.

######  Lemma IV.9.

Suppose $\mathbf{y}\in\mathbf{L}^{\infty}(\Omega)$. Let
$\mathbf{a}\in\mathbf{L}^{\infty}(\Omega)$ be non-negative. Consider the
linear bounded operator
$\mathbf{B}:\mathbf{L}^{2}(\Omega)\rightarrow\mathbf{L}^{2}(\Omega)$ defined
by

$\displaystyle({\color[rgb]{0,0,0}\mathbf{B}}\mathbf{y})(\mathbf{x})=\begin{bmatrix}-a_{1}(\mathbf{x})y_{1}(\mathbf{x})+a_{2}(\mathbf{x})y_{2}(\mathbf{x})\\\
a_{1}(\mathbf{x})y_{1}(\mathbf{x})-a_{2}(\mathbf{x})y_{2}(\mathbf{x})\end{bmatrix}$

for almost every $\mathbf{x}\in\Omega$ and all
$\mathbf{y}\in\mathbf{L}^{2}(\Omega)$. Let
$(\mathcal{T}^{\mathbf{C}}(t))_{t\geq 0}$ be the semigroup generated by the
operator $\mathbf{C}=-\mathbf{A}+\mathbf{B}$. Then
$\|\mathcal{T}^{\mathbf{C}}(t)\mathbf{y}^{0}\|_{\infty}\leq
e^{\|\mathbf{a}\|_{\infty}t}\|\mathbf{y}^{0}\|_{\infty}$ for all $t\geq 0$.

###### Proof.

We know that the operator $\mathbf{A}$ generates a semigroup
$(\mathcal{T}^{\mathbf{A}}(t))_{t\geq 0}$ given by

$\mathcal{T}^{\mathbf{A}}(t)=\begin{bmatrix}\mathcal{T}(t)&\mathbf{0}\\\
\mathbf{0}&\mathbf{I}\end{bmatrix}$ (19)

for all $t\geq 0$. Moreover, the semigroup
$(\mathcal{T}^{\mathbf{A}}(t))_{t\geq 0}$ satisfies
$\|\mathcal{T}^{\mathbf{A}}(t)\mathbf{y}^{0}\|_{\infty}\leq\|\mathbf{y}^{0}\|_{\infty}$
for all $\mathbf{y}^{0}\in\mathbf{L}^{\infty}(\Omega)$ and $t\geq 0$
(Corollary III.5). Additionally, we know that the semigroup
$(\mathcal{T}^{\mathbf{B}}(t))_{t\geq 0}$ generated by the bounded operator
$\mathbf{B}$ satisfies the estimate
$\|\mathcal{T}^{\mathbf{B}}(t)\mathbf{y}^{0}\|_{\infty}\leq
e^{\|\mathbf{a}\|_{\infty}t}\|\mathbf{y}^{0}\|_{\infty}$.  Since $\mathbf{B}$
is a bounded operator, and the resolvent of $\mathbf{C}$ has an explicit well-
defined representation [25][p.160], we can conclude that $\lambda-\mathbf{C}$
has a dense range in $\mathbf{L}^{2}(\Omega)$ for all $\lambda$. Then the
result follows from the Lie-Trotter product formula [25][Corollary III.5.8] by
noting that $\mathcal{T}^{\mathbf{C}}(t)=\lim_{N\rightarrow
0}(\mathcal{T}^{\mathbf{A}}(\frac{t}{N})\mathcal{T}^{\mathbf{B}}(\frac{t}{N}))^{N}$,
where the limit holds in the strong operator topology, for all $t\geq 0$. ∎

Now we can show that the $\mathbf{L}^{\infty}-$ estimate proved in Lemma IV.9
can be extended to a class of non-autonomous linear systems that can be
treated as autonomous linear systems over certain intervals of time.

######  Lemma IV.10.

Suppose $\mathbf{y}^{0}\in\mathbf{L}^{\infty}(\Omega)$, $c>0$ and $T>0$. Let
$a_{1},a_{2}\in L^{2}(0,T;L^{2}(\Omega))$ be non-negative and piecewise
constant with respect to time, with $\|a_{1}(t)\|_{\infty}\leq c$ and
$\|a_{2}(t)\|_{\infty}\leq c$ for all $t\in[0,T]$. Then suppose $\mathbf{y}\in
C([0,T];\mathbf{L}^{2}(\Omega))$ is given by

$\displaystyle y_{1}(\cdot,t)$ $\displaystyle=$
$\displaystyle\mathcal{T}(t)y^{0}_{1}-\int_{0}^{t}\mathcal{T}(t-s)\bigg{(}a_{1}(\cdot,t)y_{1}(\cdot,s)\bigg{)}ds$
$\displaystyle+\int_{0}^{t}\mathcal{T}(t-s)\bigg{(}a_{2}(\cdot,t)y_{2}(\cdot,s)\bigg{)}ds$
$\displaystyle y_{2}(\cdot,t)$ $\displaystyle=$ $\displaystyle
y^{0}_{2}+\int_{0}^{t}a_{1}(\cdot,s)y_{1}(\cdot,s)ds-\int_{0}^{t}a_{2}(\cdot,s)y_{2}(\cdot,s)ds$

for all $t\in[0,T]$. Then

${\color[rgb]{0,0,0}\|\mathbf{y}(\cdot,t)\|_{\infty}}\leq
e^{ct}\|\mathbf{y}^{0}\|_{\infty}$ (20)

for all $t\in[0,T]$.

###### Proof.

Let $(t_{i})_{i=0}^{m}$ be a finite sequence of length $m+1\in\mathbb{Z}_{+}$
of strictly increasing time instants, with $t_{0}=0$, such that the functions
$a_{1}$ and $a_{2}$ are constant over the intervals $[t_{i-1},t_{i})$,
$i\in\\{1,...,m\\}$. Then, for each $i\in\\{1,...,m\\}$, consider the bounded
operators
$\mathbf{B}_{i}:\mathbf{L}^{2}(\Omega)\rightarrow\mathbf{L}^{2}(\Omega)$ and
$\mathbf{C}_{i}:\mathcal{D}(\mathbf{A})\rightarrow\mathbf{L}^{2}(\Omega)$

$\displaystyle(\mathbf{B}_{i}\mathbf{y})(\mathbf{x})=\begin{bmatrix}-a_{1}(\mathbf{x},t_{i-1})y_{1}(\mathbf{x})+a_{2}(\mathbf{x},t_{i-1})y_{2}(\mathbf{x})\\\
a_{1}(\mathbf{x},t_{i-1})y_{1}(\mathbf{x})-a_{2}(\mathbf{x},t_{i-1})y_{2}(\mathbf{x})\end{bmatrix}$
(21)

for almost every $\mathbf{x}\in\Omega$ and all
$\mathbf{y}\in\mathbf{L}^{2}(\Omega)$, and
$\mathbf{C}_{i}=\mathbf{A}+\mathbf{B}_{i}$, respectively. Then for each
$i\in\\{1,...,m\\}$, $\mathbf{y}$ is given by

$\mathbf{y}(\cdot,t)=\mathcal{T}^{\mathbf{C}_{i}}(t-t_{i})\mathcal{T}^{\mathbf{C}_{i-1}}(t_{i}-t_{i-1})...\mathcal{T}^{\mathbf{C}_{1}}(t_{1})$
(22)

for all $t\in[t_{i-1},t_{i}]$. Then the result follows from Lemma IV.9. ∎

######  Lemma IV.11.

Suppose $\mathbf{y}^{0}\in\mathbf{L}^{\infty}(\Omega)$, $c>0$, and $T>0$. Let
the functions $a_{1},a_{2}\in L^{2}(0,T;L^{2}(\Omega))$ be non-negative with
$\|a_{1}(t)\|_{\infty}\leq c$ and $\|a_{2}(t)\|_{\infty}\leq c$ for almost
every $t\in[0,T]$, Then  there exists $\mathbf{y}\in
C([0,T];\mathbf{L}^{2}(\Omega))$ given by

$\displaystyle y_{1}(\cdot,t)$ $\displaystyle=$
$\displaystyle\mathcal{T}(t)y^{0}_{1}-\int_{0}^{t}\mathcal{T}(t-s)\Big{(}a_{1}(\cdot,s)y_{1}(\cdot,s)\Big{)}ds$
$\displaystyle+\int_{0}^{t}\mathcal{T}(t-s)\Big{(}a_{2}(\cdot,s)y_{2}(\cdot,s)\Big{)}ds,$
$\displaystyle y_{2}(\cdot,t)$ $\displaystyle=$ $\displaystyle
y^{0}_{2}+\int_{0}^{t}a_{1}(\cdot,s)y_{1}(\cdot,s)ds-\int_{0}^{t}a_{2}(\cdot,s)y_{2}(\cdot,s)ds$

for all $t\geq 0$.  Moreover,

$\|\mathbf{y}(\cdot,t)\|_{\infty}\leq e^{ct}\|\mathbf{y}^{0}\|_{\infty}$ (24)

for all $t\in[0,T]$.

###### Proof.

Given that $a_{1},a_{2}\in L^{2}(0,T;L^{2}(\Omega))$, we know that there
exists a sequence of  piecewise constant (with respect to time) non-negative
functions $(a_{1}^{i})_{i=1}^{\infty},(a_{2}^{i})_{i=1}^{\infty}$ in
$L^{2}(0,T;L^{2}(\Omega))$ such that
$\lim_{i\rightarrow\infty}\|a_{j}^{i}-a_{j}\|_{L^{2}(0,T;L^{2}(\Omega))}=0$,
for $j=1,2$ [55][Proposition 1.36]. Moreover, for each $j\in\\{1,2\\}$, we can
assume that $\|a_{j}^{i}(t)\|_{\infty}\leq c$ for all $t\in[0,T]$ and all
$i\in\mathbb{Z}_{+}$. Consider the corresponding sequence
$(\mathbf{y})_{i=1}^{\infty}$ in $C([0,T];\mathbf{L}^{2}(\Omega))$ defined by

$\displaystyle y^{i}_{1}(\cdot,t)$ $\displaystyle=$
$\displaystyle\mathcal{T}(t)y^{0}_{1}-\int_{0}^{t}\mathcal{T}(t-s)\Big{(}a^{i}_{1}(\cdot,s)y^{i}_{1}(\cdot,s)\Big{)}ds$
$\displaystyle+\int_{0}^{t}\mathcal{T}(t-s)\Big{(}a_{2}(\cdot,s)y^{i}_{2}(\cdot,s)\Big{)}ds,$
$\displaystyle y^{i}_{2}(\cdot,t)$ $\displaystyle=$ $\displaystyle
y^{0}_{2}+\int_{0}^{t}a^{i}_{1}(\cdot,s)y^{i}_{1}(\cdot,s)ds-\int_{0}^{t}a^{i}_{2}(\cdot,s)y^{i}_{2}(\cdot,s)ds$

for each $i\in\mathbb{Z}_{+}$. Let ${\color[rgb]{0,0,0}\mathbf{e}^{i,j}}\in
C([0,T];\mathbf{L}^{2}(\Omega))$ be given by
$\mathbf{e}^{i,j}=\mathbf{y}^{i}-\mathbf{y}^{j}$ for each
$i,j\in\mathbb{Z}_{+}$. Then, from equations (LABEL:eq:vareq1) and
(LABEL:eq:vareqseq), we know that  $\mathbf{e}^{i,j}$ satisfies

$\displaystyle{\color[rgb]{0,0,0}e^{i,j}_{1}}(\cdot,t)$ $\displaystyle=$
$\displaystyle-\int_{0}^{t}\mathcal{T}(t-s)\Big{(}a^{i}_{1}(\cdot,s)y^{i}_{1}(s,\cdot)\Big{)}ds$
$\displaystyle+\int_{0}^{t}\mathcal{T}(t-s)\Big{(}a^{i}_{2}(\cdot,s)y^{i}_{2}(s,\cdot)\Big{)}ds$
$\displaystyle+\int_{0}^{t}\mathcal{T}(t-s)\Big{(}a^{{\color[rgb]{0,0,0}j}}_{1}(\cdot,s)y^{{\color[rgb]{0,0,0}j}}_{1}(\cdot,s)\Big{)}ds$
$\displaystyle-\int_{0}^{t}\mathcal{T}(t-s)\Big{(}a^{{\color[rgb]{0,0,0}j}}_{2}(\cdot,s)y^{{\color[rgb]{0,0,0}j}}_{2}(\cdot,s)\Big{)}ds$
$\displaystyle=$
$\displaystyle-\int_{0}^{t}\mathcal{T}(t-s)\Big{(}(a^{i}_{1}(\cdot,s)-a^{{\color[rgb]{0,0,0}j}}_{1}(\cdot,s))y^{i}_{1}(\cdot,s)\Big{)}ds$
$\displaystyle+\int_{0}^{t}\mathcal{T}(t-s)\Big{(}a^{{\color[rgb]{0,0,0}j}}_{1}(\cdot,s)(y^{{\color[rgb]{0,0,0}j}}_{1}(\cdot,s)-y^{i}_{1}(s,\cdot))\Big{)}ds$
$\displaystyle+\int_{0}^{t}\mathcal{T}(t-s)\Big{(}(a^{i}_{2}(\cdot,s)-a^{{\color[rgb]{0,0,0}j}}_{2}(\cdot,t))y^{i}_{2}(\cdot,s)\Big{)}ds$
$\displaystyle-\int_{0}^{t}\mathcal{T}(t-s)\Big{(}a^{{\color[rgb]{0,0,0}j}}_{2}(\cdot,s)(y^{{\color[rgb]{0,0,0}j}}_{2}(\cdot,s)-y^{i}_{2}(\cdot,s))\Big{)}ds$

for all $t\in[0,T]$. Considering the fact that the semigroup $\mathcal{T}(t)$
is contractive (Proposition III.3), and that $a_{{\color[rgb]{0,0,0}k}}^{i}$
and $y^{i}_{{\color[rgb]{0,0,0}k}}$ are uniformly bounded in
$L^{\infty}((0,T)\times\Omega)$ for $k=1,2$, we can conclude that there exists
a constant $\alpha>0$ such that

$\displaystyle\|e^{i,j}_{1}(\cdot,t)\|_{2}$ $\displaystyle\leq$
$\displaystyle\alpha\|a_{1}^{i}-a^{j}_{1}\|_{\infty}\|y^{0}_{1}\|_{\infty}$
(26)
$\displaystyle+~{}\alpha\|a^{j}_{1}\|_{\infty}\int_{0}^{t}\|e^{i,j}_{1}(s)\|_{2}ds$
$\displaystyle+~{}\alpha\|a_{2}^{i}-a^{j}_{2}\|_{\infty}\|y^{0}_{2}\|_{\infty}$
$\displaystyle+~{}\alpha\|a^{j}_{2}\|_{\infty}\int_{0}^{t}\|e^{i,j}_{2}(s)\|_{2}ds$

for all $t\in[0,T]$. Similarly, we can obtain the estimate

$\displaystyle\|e^{i,j}_{2}(\cdot,t)\|_{2}$ $\displaystyle\leq$
$\displaystyle\alpha\|a_{1}^{i}-a^{j}_{1}\|_{\infty}\|y^{0}_{1}\|_{\infty}$
(27)
$\displaystyle+~{}\alpha\|a^{j}_{1}\|_{\infty}\int_{0}^{t}\|e^{i,j}_{1}(s)\|_{2}ds$
$\displaystyle+~{}\alpha\|a_{2}^{i}-a^{j}_{2}\|_{\infty}\|y^{0}_{2}\|_{\infty}$
$\displaystyle+~{}\alpha\|a^{j}_{2}\|_{\infty}\int_{0}^{t}\|e^{i,j}_{2}(s)\|_{2}ds$

for all $t\in[0,T]$.

Then, by considering the sum
$\|e^{i,j}_{1}(\cdot,t)\|_{2}+\|e^{i,j}_{2}(\cdot,t)\|_{2}$, combining the two
inequalities (26) and (27), and applying the integral form of Gronwall’s
inequality [27], we have that

$\|e^{i,j}_{1}(\cdot,t)\|_{2}+\|e^{i,j}_{2}(\cdot,t)\|_{2}~{}\leq~{}C_{1}e^{C_{2}t}$
(28)

for all $t\in[0,T]$, where $C_{1}$ and $C_{2}$ are constants depending only on
$\beta$ (the upper bound on the functions $r_{i}$) and
$\|\mathbf{y}^{0}\|_{\infty}$. From the inequality (28), we can infer that

$\lim_{i,j\rightarrow\infty}\|\mathbf{e}^{i,j}\|_{C([0,T];\mathbf{L}^{2}(\Omega))}=0$

This implies that the solution of (LABEL:eq:vareq1) is continuous with respect
to $a_{i}$, and by continuity, one can construct a solution $\mathbf{y}$
corresponding to the coefficients $a_{i}$. Considering the estimate (20), we
can conclude that $\mathbf{y}$ satisfies the estimate (24). ∎

From the above lemma, we can conclude the following theorem on global
existence of solutions of the PDE (17).

######  Theorem IV.12.

Suppose $\mathbf{y}^{0}\in\mathbf{L}^{\infty}(\Omega)$. Then the PDE (17) has
a unique global mild solution.

###### Proof.

We use a contraction mapping approach to construct the solution. Consider a
map $\Gamma$ on $C([0,T];\mathbf{L}^{2}(\Omega))$ that is defined by
$\mathbf{v}\mapsto\Gamma(\mathbf{v})\equiv\tilde{\mathbf{v}}$, where
$\tilde{\mathbf{v}}$ is constructed by setting $a_{i}=F_{i}(v_{2})$,
$i\in\\{1,2\\}$, in (LABEL:eq:vareq1):

$\displaystyle\tilde{v}_{1}(\cdot,t)$ $\displaystyle=$
$\displaystyle\mathcal{T}(t)y^{0}_{1}-\int_{0}^{t}\mathcal{T}(t-s)\Big{(}F_{1}(v_{2}(\cdot,s))\tilde{v}_{1}(\cdot,s)\Big{)}ds$
$\displaystyle+~{}\int_{0}^{t}\mathcal{T}(t-s)\Big{(}F_{2}(v_{2}(\cdot,s))\tilde{v}_{2}(\cdot,s)\Big{)}ds,$
$\displaystyle\tilde{v}_{2}(\cdot,t)$ $\displaystyle=$ $\displaystyle
y^{0}_{2}+\int_{0}^{t}F_{1}(v_{2}(\cdot,s))\tilde{v}_{1}(\cdot,s)ds$ (29)
$\displaystyle-\int_{0}^{t}F_{2}(v_{2}(\cdot,s))\tilde{v}_{2}(\cdot,s)ds$

for all $t\in[0,T]$. Then $\mathbf{y}$ is a solution of the PDE (17) if it is
a fixed point of the map $\Gamma$. From Lemma IV.11, any solution of the PDE
(17) must be attained in the set $\mathcal{Y}=\\{\mathbf{v}\in
C([0,T],\mathbf{L}^{2}(\Omega));\|v(t)\|_{\infty}\leq\hat{C}~{}\forall
t\in[0,T]\\}$ for some sufficiently large constant $\hat{C}>0$ that depends on
$\|\mathbf{y}^{0}\|_{\infty}$. Therefore, it suffices to check that the map
$\Gamma$ is a contraction on $\mathcal{Y}$, which guarantees that it has a
unique fixed point.

Let $\mathbf{v},\mathbf{u}\in\mathcal{Y}$. To show that $\Gamma$ is a
contraction on $\mathcal{Y}$, we will derive an upper bound on
$\sup_{t\in[0,T]}\|\Gamma(\mathbf{v})(t)-\Gamma(\mathbf{u})(t)\|_{2}$. In
order to estimate $\|\tilde{v}_{1}(\cdot,t)-\tilde{u}_{1}(\cdot,t)\|_{2}$, we
first compute an upper bound on the corresponding terms in the difference
according to the computation in (29):

$\displaystyle\bigg{\|}-\int_{0}^{t}\mathcal{T}(t-s)\Big{(}F_{1}(v_{2}(\cdot,s))\tilde{v}_{1}(\cdot,s)\Big{)}ds$
$\displaystyle~{}~{}~{}~{}~{}+\int_{0}^{t}\mathcal{T}(t-s)\Big{(}F_{1}(u_{2}(\cdot,s))\tilde{u}_{1}(\cdot,s)\Big{)}ds\bigg{\|}_{2}~{}~{}~{}\leq$
$\displaystyle\bigg{\|}\int_{0}^{t}\mathcal{T}(t-s)\big{(}-F_{1}(v_{2}(\cdot,s))\tilde{v}_{1}(\cdot,s)$
$\displaystyle~{}~{}~{}~{}~{}+~{}F_{1}(u_{2}(\cdot,s))\tilde{v}_{1}(\cdot,s)\big{)}ds\bigg{\|}_{2}~{}~{}~{}+$
$\displaystyle\bigg{\|}\int_{0}^{t}\mathcal{T}(t-s)\big{(}F_{1}(u_{2}(\cdot,s))\tilde{v}_{1}(\cdot,s)-F_{1}(u_{2}(\cdot,s))\tilde{u}_{1}(\cdot,s)\big{)}ds\bigg{\|}_{2}$

Since $r_{i}$ in the definition (16) of $F_{i}$ is a bounded Lipschitz
function (with, say, Lipschitz constant $K>0$), and $\mathcal{T}(t)$ is a
contraction, we can compute the following upper bounds on the left-hand side
of the above inequality:

$\displaystyle\bigg{|}\int_{0}^{t}\|\big{(}-F_{1}(v_{2}(\cdot,s))+F_{2}(u_{2}(\cdot,s)\Big{)}\tilde{v}_{1}(\cdot,s)\|_{2}ds$
$\displaystyle~{}~{}~{}~{}~{}+~{}c\int_{0}^{t}\|\tilde{v}_{1}(\cdot,s)-\tilde{u}_{1}(\cdot,s)\|_{2}ds\bigg{|}$
$\displaystyle\leq~{}~{}K\int_{0}^{t}\|\Big{(}-v_{2}(\cdot,s)+u_{2}(\cdot,s)\Big{)}\tilde{v}_{1}(\cdot,s)\|_{2}ds$
$\displaystyle~{}~{}~{}~{}~{}+~{}c\int_{0}^{t}\|\tilde{v}_{1}(\cdot,s)-\tilde{u}_{1}(\cdot,s)\|_{2}ds$
$\displaystyle\leq~{}~{}\hat{C}K\int_{0}^{t}\|v_{2}(\cdot,s)-u_{2}(\cdot,s)\|_{2}ds$
$\displaystyle~{}~{}~{}~{}~{}+~{}c\int_{0}^{t}\|\tilde{v}_{1}(\cdot,s)-\tilde{u}_{1}(\cdot,s)\|_{2}ds$
$\displaystyle\leq~{}~{}TK\sup_{t\in[0,T]}\|v_{2}(\cdot,t)-u_{2}(\cdot,t)\|_{2}$
$\displaystyle~{}~{}~{}~{}~{}+~{}Tc\sup_{t\in[0,T]}\|\tilde{v}_{1}(\cdot,t)-\tilde{u}_{1}(\cdot,t)\|_{2}.$

Using a similar computation, we can estimate that

$\displaystyle\sup_{t\in[0,T]}(\|\tilde{v}_{1}(\cdot,t)-\tilde{u}_{1}(\cdot,t)\|_{2}+\|\tilde{v}_{2}(\cdot,t)-\tilde{u}_{2}(\cdot,t)\|_{2})$
$\displaystyle\hskip
28.45274pt\leq~{}~{}~{}~{}~{}\tilde{C}T\sup_{t\in[0,T]}\|v_{1}(\cdot,t)-u_{1}(\cdot,t)\|_{2}$
$\displaystyle\hskip
42.67912pt+~{}\tilde{C}T\sup_{t\in[0,T]}\|v_{2}(\cdot,t)-u_{2}(\cdot,t)\|_{2}$
$\displaystyle\hskip
42.67912pt+~{}\tilde{C}T\sup_{t\in[0,T]}\|\tilde{v}_{1}(\cdot,t)-\tilde{u}_{1}(\cdot,t)\|_{2}$
$\displaystyle\hskip
42.67912pt+~{}\tilde{C}T\sup_{t\in[0,T]}\|\tilde{v}_{2}(\cdot,t)-\tilde{u}_{2}(\cdot,t)\|_{2}$

for a sufficiently large constant $\tilde{C}>0$ that depends only on $K$, $c$,
and $\hat{C}$. This implies that

$\displaystyle\sup_{t\in[0,T]}(\|\tilde{v}_{1}(\cdot,t)-\tilde{u}_{1}(\cdot,t)\|_{2}+\|\tilde{v}_{2}(\cdot,t)-\tilde{u}_{2}(\cdot,t)\|_{2})$
$\displaystyle\hskip
21.33955pt\leq~{}~{}~{}~{}\frac{\tilde{C}T}{1-\tilde{C}T}\sup_{t\in[0,T]}\|v_{1}(\cdot,t)-u_{1}(\cdot,t)\|_{2}$
$\displaystyle\hskip
28.45274pt~{}+~{}\frac{\tilde{C}T}{1-\tilde{C}T}\sup_{t\in[0,T]}\|v_{2}(\cdot,t)-u_{2}(\cdot,t)\|_{2}.$

Thus, if $T>0$ is small enough, $\Gamma$ is a contraction map on
$\mathcal{Y}$. This implies that for small enough $T>0$, $\Gamma$ has a fixed
point $\mathbf{y}$ that is a unique local mild solution to the PDE (17). Then,
due to the uniform bound on the solution established in Lemma IV.11, it
follows that the solution can, in fact, be extended to arbitrary $T>0$, and is
therefore global. ∎

Next, our goal will be to prove that $\mathbf{y}^{d}$ is the globally
asymptotically stable equilibrium of the system (17). Toward this end, we
first prove the following preliminary results.

######  Lemma IV.13.

Suppose $y^{0}\in L^{\infty}(\Omega)$ and $T>0$. Let $a\in L^{\infty}(\Omega)$
be non-negative. Consider the multiplication operator
$B:L^{2}(\Omega)\rightarrow L^{2}(\Omega)$ defined by

$\displaystyle(By)(\mathbf{x})=-a(\mathbf{x})y(\mathbf{x})$

for all $\mathbf{x}\in\Omega$ and all $y\in L^{2}(\Omega)$. Let
$(\mathcal{T}^{C}(t))_{t\geq 0}$ be the semigroup generated by the operator
$C=-A+B$. Then $\|\mathcal{T}^{C}(t)y^{0}\|_{\infty}\leq\|y^{0}\|_{\infty}$
for all $t\geq 0$.

###### Proof.

We know that if $(\mathcal{T}(t))_{t\geq 0}$ is the semigroup generated by the
operator $-A$, then from Corollary III.5,
$\|{\color[rgb]{0,0,0}\mathcal{T}(t)}\mathbf{y}^{0}\|_{\infty}\leq\|y^{0}\|_{\infty}$
for all $t\geq 0$. Moreover, $B$ generates the multiplication semigroup
$(e^{-a(\cdot)t})_{t\geq 0}$. Then the result follows from the Lie-Trotter
formula [25]. ∎

######  Lemma IV.14.

Let $T>0$. Let $f,a\in L^{2}(0,T;L^{2}(\Omega))$ be non-negative functions
defined such that $\|f(t)\|_{\infty}$ and $\|a(t)\|_{\infty}$ are bounded by a
constant $C>0$ almost everywhere on $t\in[0,T]$. Suppose $e\in
C([0,T];L^{2}(\Omega))$ is given by

$\displaystyle e(\cdot,t)=$
$\displaystyle-\int_{0}^{t}\mathcal{T}(t-s)\Big{(}a(\cdot,s)e(\cdot,s)\Big{)}ds$
$\displaystyle+\int_{0}^{t}\mathcal{T}(t-s)f(\cdot,s)ds$

for all $t\in[0,T]$. Then $e(\cdot,t)$ is non-negative for all $t\in[0,T]$.

###### Proof.

Since the proof follows a similar line of argument as the proof of Lemma
IV.11, we only sketch an outline of the proof here. As in the proof of Lemma
IV.11, for a given $a\in L^{2}(0,T;L^{2}(\Omega))$ we can construct a sequence
$(a^{i})_{i=1}^{\infty}$ in $L^{2}(0,T;L^{2}(\Omega))$ that is piecewise
constant in time and converges in $L^{2}(0,T;L^{2}(\Omega))$, with
$\|a^{i}(t)\|_{\infty}$ bounded almost everywhere on $[0,T]$ by $C>0$. Let the
sequence $(e_{i})_{i=1}^{\infty}$ in $C([0,T];L^{2}(\Omega))$ be given by

$\displaystyle e_{i}(\cdot,t)=$
$\displaystyle-\int_{0}^{t}\mathcal{T}(t-s)\Big{(}a^{i}(\cdot,s)e_{i}(\cdot,s)\Big{)}ds$
$\displaystyle+~{}\int_{0}^{t}\mathcal{T}(t-s)f(\cdot,s)ds$

for all $t\in[0,T]$.  Since $e_{i}(t)$ is also the solution of the PDE
$\dot{e}_{i}(t)=-Ae_{i}(t)-a^{i}(\cdot,t)e_{i}(t)+f(\cdot,t)$ with initial
condition equal to $0$, and this solution can be constructed using the
positive semigroup $\mathcal{T}^{C}(t)$ in Lemma IV.13, we can conclude that
$(e_{i})_{i=1}^{\infty}$ is non-negative for each $i\in\mathbb{Z}_{+}$. Then,
using the fact that the sequences $(e_{i})_{i=1}^{\infty}$ and
$(a^{i})_{i=1}^{\infty}$ are uniformly bounded in the spaces
$C([0,T];L^{2}(\Omega))$ and $L^{2}(0,T;L^{2}(\Omega))$, respectively, and
applying Gronwall’s lemma, the result follows. ∎

We can use Lemma IV.14 prove the next result, which will enable us to show
later on that the rate of convergence of the solution $\mathbf{y}$ of the PDE
(17) toward $\mathbf{0}$ can be controlled by the rate of convergence of the
solution of a related linear PDE.

######  Theorem IV.15.

(Comparison Principle)
Let $T>0$. Let $y^{0}\in L^{2}(\Omega)$ and $f,g\in L^{2}(0,T;L^{2}(\Omega))$
be non-negative such that $\|f(t)\|_{\infty}$ and $\|g(t)\|_{\infty}$ are
bounded by a constant $C_{1}>0$ almost everywhere on $t\in[0,T]$. Define the
operator $C=-{A}-\|g\|_{\infty}\mathbf{I}$. Let $y(\cdot,t)$ be given by

$\displaystyle y(\cdot,t)=$
$\displaystyle\mathcal{T}(t)y^{0}-\int_{0}^{t}\mathcal{T}(t-s)\Big{(}g(\cdot,s)y(\cdot,s)\Big{)}ds$
$\displaystyle+\int_{0}^{t}\mathcal{T}(t-s)f(\cdot,s)ds$ (30)

for all $t\in[0,T]$. Then $y(\cdot,t)\geq\mathcal{T}^{C}(t)y^{0}$ for all
$t\in[0,T]$, where $(\mathcal{T}^{C}(t))_{t\geq 0}$ is the semigroup generated
by $C$.

###### Proof.

Let $\tilde{y}(\cdot,t)=\mathcal{T}^{C}(t)y^{0}$ for all $t\geq 0$. Then, we
know that $\tilde{y}(\cdot,t)$ satisfies the equation

$\tilde{y}(\cdot,t)=\mathcal{T}(t)y^{0}-\int_{0}^{t}\mathcal{T}(t-s)\|g\|_{\infty}\tilde{y}(\cdot,s)ds$
(31)

for all $t\in[0,T]$. Let $e=y-\tilde{y}$.  Substituting equations (IV.15) and
(31) for $y$ and $\tilde{y}$, respectively, we have that

$\displaystyle e(\cdot,t)$ $\displaystyle=$
$\displaystyle-\int_{0}^{t}\mathcal{T}(t-s)\Big{(}g(\cdot,s)y(\cdot,s)\Big{)}ds$
$\displaystyle+\int_{0}^{t}\mathcal{T}(t-s)f(\cdot,s)ds$
$\displaystyle+\int_{0}^{t}\mathcal{T}(t-s)\|g\|_{\infty}\tilde{y}(\cdot,s)ds$

for all $t\in[0,T]$. By adding and subtracting the term
$\int_{0}^{t}\mathcal{T}(t-s)\Big{(}g(\cdot,s)\tilde{y}(\cdot,s)\Big{)}ds$ in
this expression for $e$, we obtain

$\displaystyle e(\cdot,t)$ $\displaystyle=$
$\displaystyle-\int_{0}^{t}\mathcal{T}(t-s)\Big{(}(g(\cdot,s))e(\cdot,s)\Big{)}ds$
$\displaystyle+\int_{0}^{t}\mathcal{T}(t-s)f(\cdot,s)ds$
$\displaystyle+\int_{0}^{t}\mathcal{T}(t-s)\Big{(}(\|g\|_{\infty}-g(\cdot,s))\tilde{y}(\cdot,s)\Big{)}ds$

for all $t\in[0,T]$. Then the result follows from the non-negativity of $e$,
which is a consequence of Lemma IV.14. ∎

######  Theorem IV.16.

(Positive Lower Bound on Solutions) Let $T>0$. Let $y^{0}\in L^{2}(\Omega)$
and $f,g\in L^{2}(0,T;L^{2}(\Omega))$ be non-negative such that
$\|f(t)\|_{\infty}$ and $\|g(t)\|_{\infty}$ are bounded by a constant
$C_{1}>0$ almost everywhere on $t\in[0,T]$. Let $y(\cdot,t)$ be given by

$\displaystyle y(\cdot,t)$ $\displaystyle=$
$\displaystyle\mathcal{T}(t)y^{0}-\int_{0}^{t}\mathcal{T}(t-s)\Big{(}g(\cdot,s)y(\cdot,s)\Big{)}ds$
(32) $\displaystyle+\int_{0}^{t}\mathcal{T}(t-s)f(\cdot,s)ds$

for all $t\in[0,T]$. Then there exist constants $\tau,\epsilon,\delta>0$,
independent of $y^{0}$ and $T>0$, such that if $\tau+\delta<T$, then
$y(\cdot,t)\geq\epsilon\|y^{0}\|_{1}$ for all $t\geq[\tau,\tau+\delta]$.

###### Proof.

We know from Theorem IV.5 that there exists a constant $k>0$ and time $T>0$,
independent of $y^{0}$, such that $\mathcal{T}(t)y^{0}\geq k\|y^{0}\|_{1}$ for
all $t\geq T$. Let $C=-A-\|g\|_{\infty}\mathbf{I}$. Then the semigroup
$(\mathcal{T}^{C}(t))$ generated by the operator $C$ is given by
$\mathcal{T}(t)=e^{-\|g\|_{\infty}t}\mathcal{T}(t)$ for all $t\geq 0$. The
result then follows from Theorem IV.15. ∎

The following theorem states the fundamental result that the PDE (17)
conserves mass and maintains positivity.

######  Theorem IV.17.

Let ${\color[rgb]{0,0,0}\mathbf{y}^{0}}\in\mathbf{L}^{\infty}(\Omega)$ be non-
negative. Then the unique global mild solution of the PDE (17) is non-
negative, and $\|\mathbf{y}(\cdot,t)\|_{1}=\|\mathbf{y}^{0}\|$ for all $t\geq
0$.

###### Proof.

The conservation of mass is a simple consequence of taking the inner product
of the solution of (17) with a constant function. The positivity property of
solutions follows from [18][Theorem 1] by noting that, if $\lambda>0$ is large
enough, then $G(\mathbf{y})+\lambda\mathbf{y}\geq 0$ for all
$\mathbf{y}\in\mathbf{L}^{2}(\Omega)$ that are non-negative. ∎

We now require some additional notation. For a function $f\in
L^{2}{(\Omega)}$, we define $f_{+}:=\frac{|f|+f}{2}$, the projection of $f$
onto the set of non-negative functions in $L^{2}(\Omega)$, and
$f_{-}:=-\frac{|f|-f}{2}$, the projection of $f$ onto the set of non-positive
functions in $L^{2}(\Omega)$. Given these definitions, we have the following
result on partial monotonicity of solutions of the PDE (17).

######  Proposition IV.18.

(Partial Monotonicity of Solutions) Let
${\color[rgb]{0,0,0}\mathbf{y}^{0}}\in\mathbf{L}^{\infty}(\Omega)$ be
positive. Then, for all $t\geq s\geq 0$, the unique global mild solution of
the PDE (17) satisfies

$\displaystyle(y^{d}-y_{2}(\cdot,t))_{+}~{}\leq~{}(y^{d}-y_{2}(\cdot,s))_{+}$
(33)
$\displaystyle(y^{d}-y_{2}(\cdot,t))_{-}~{}\geq~{}(y^{d}-y_{2}(\cdot,s))_{-}$
(34)

###### Proof.

We will only prove the first inequality (33). Since
$\mathbf{y}^{0}\in\mathbf{L}^{\infty}(\Omega)$, we know that $y_{2}\in
C([0,1];L^{2}(\Omega))$ and $\|y_{2}(t)\|_{\infty}$ is uniformly bounded over
$[0,T]$. Assume that $y^{d}-y^{0}_{2}$ is non-zero and non-negative on a set
$\Omega_{1}\subseteq\Omega$ of positive measure. For the sake of
contradiction, suppose that there exists $t_{2}\in(0,T]$ such that
$y_{2}(\cdot,t_{2})$ is greater than $y^{d}$ on a subset of $\Omega_{1}$ that
has positive Lebesgue measure. Then, due to the fact that $y_{2}\in
C([0,T];L^{2}(\Omega))$, there must exist $t_{1}\in(0,t_{2})$ and a measurable
set $\Omega_{2}\subset\Omega_{1}$ of positive Lebesgue measure, such that for
each $s\in[t_{1},t_{2}]$, $y_{2}(\mathbf{x},s)\geq y^{d}(\mathbf{x})$ for
almost every $\mathbf{x}\in\Omega_{2}$, with $y_{2}(\mathbf{x},t_{2})\neq
y_{2}(\mathbf{x},s)$ for almost every $\mathbf{x}\in\Omega_{2}$ and a subset
of $[t_{1},t_{2}]$ with positive Lebesgue measure. However, we know that

$\displaystyle y_{2}(\cdot,t)$ $\displaystyle=$ $\displaystyle
y_{2}(\cdot,t_{1})+\int_{{\color[rgb]{0,0,0}t_{1}}}^{t}F_{1}(y_{2}(\cdot,\tau))y_{1}(\cdot,\tau)d\tau$
(35)
$\displaystyle-\int_{{\color[rgb]{0,0,0}t_{1}}}^{t}F_{2}(y_{2}(\cdot,\tau)){\color[rgb]{0,0,0}y_{2}(\cdot,\tau)}d\tau$

for all $t\in[t_{1},t_{2}]$ This implies that

$\displaystyle y_{2}(\mathbf{x},t)$ $\displaystyle=$ $\displaystyle
y_{2}(\mathbf{x},t_{1})+\int_{{\color[rgb]{0,0,0}t_{1}}}^{t}F_{1}(y_{2}(\mathbf{x},\tau))y_{1}(\mathbf{x},\tau)d\tau$
$\displaystyle-\int_{{\color[rgb]{0,0,0}t_{1}}}^{t}F_{2}(y_{2}(\mathbf{x},\tau))y_{1}(\mathbf{x},\tau)d\tau$
$\displaystyle=$ $\displaystyle
y_{2}(\mathbf{x},t_{1})-\int_{{\color[rgb]{0,0,0}t_{1}}}^{t}r_{2}(y_{2}(\mathbf{x},\tau)-y^{d}(\mathbf{x})){\color[rgb]{0,0,0}y_{2}(\mathbf{x},\tau)}d\tau$

for almost every $\mathbf{x}\in\Omega_{2}$ and for all $t\in[t_{1},t_{2}]$.
The function $y_{1}$ is non-negative due to Theorem IV.17. Moreover, $r_{2}$
is also non-negative by definition. Thus, we arrive at the contradiction that
$y_{2}(\mathbf{x},t)\leq y_{2}(\mathbf{x},t_{1})$ for almost every
$\mathbf{x}\in\Omega_{1}$ and for all $t\in[t_{1},t_{2}]$. Hence, we must have
that

$y_{2}(\mathbf{x},t)=y_{2}(\mathbf{x},t_{1})+\int_{{\color[rgb]{0,0,0}t_{1}}}^{t}r_{1}(y_{2}(\mathbf{x},\tau)-y^{d}(\mathbf{x}))y_{1}(\mathbf{x},\tau)d\tau$

for almost every $\mathbf{x}\in\Omega_{1}$ and for all $t\in[0,T]$. This
implies that $y_{2}$ is non-decreasing with time, and is less than or equal to
$y^{d}$ almost everywhere on $\Omega_{1}$. This proves the first inequality
(33). Using a similar argument, based on the fact that $r_{1}$ and $r_{2}$ are
non-negative bounded functions, we can arrive at the second inequality (34). ∎

Using the above proposition, we will establish global asymptotic stability of
the system (17) in the $L^{1}$ norm. Towards this end, we first establish
marginal stability of the system about the equilibrium distribution
$\mathbf{y}^{d}$.

######  Theorem IV.19.

($L^{1}$-Lyapunov Stability) Let
$\mathbf{y}^{0}\in\mathbf{L}^{\infty}(\Omega)$ be positive and
$\int_{\Omega}\mathbf{y}^{0}(\mathbf{x})d\mathbf{x}=1$. For every
$\epsilon>0$, if

$\|\mathbf{y}^{0}-\mathbf{y}^{d}\|_{1}\leq\epsilon,$ (36)

then the solution $\mathbf{y}(\cdot,t)$ of the system (17) satisfies

$\|\mathbf{y}(\cdot,t)-\mathbf{y}^{d}\|_{1}\leq 2\epsilon$ (37)

for all $t\geq 0$.

###### Proof.

We know that the solution $\mathbf{y}$ satisfies

$\int_{\Omega}\mathbf{y}(\cdot,t)d\mathbf{x}=\int_{\Omega}y_{1}(\cdot,t)d\mathbf{x}+\int_{\Omega}y_{2}(\cdot,t)d\mathbf{x}=1$

for all $t\in[0,T]$. From Proposition IV.18, we know that
$\|y_{2}(\cdot,t)-y^{d}\|_{1}$ is non-decreasing with time $t$. Hence,
$\|y_{2}(\cdot,t)-y^{d}\|_{1}\leq\epsilon$ for all $t\geq 0$. Then, we have
that

$\int_{\Omega}y_{1}(\mathbf{x},t)d\mathbf{x}+\int_{\Omega}(y_{2}(\mathbf{x},t)-y^{d}(\mathbf{x}))d\mathbf{x}=1-\int_{\Omega}y^{d}(\mathbf{x})d\mathbf{x}$

for all $t\geq 0$. This implies that

$\displaystyle\int_{\Omega}y_{1}(\mathbf{x},t)d\mathbf{x}~{}$
$\displaystyle\leq$
$\displaystyle~{}-\int_{\Omega}(y_{2}(\mathbf{x},t)-y^{d}(\mathbf{x}))d\mathbf{x}$
$\displaystyle\leq$
$\displaystyle~{}\|y_{2}(\cdot,t)-y^{d}\|_{1}~{}\leq~{}\epsilon$

for all $t\geq 0$. This concludes the proof. ∎

######  Proposition IV.20.

Let $\mathbf{y}^{0}\in\mathbf{L}^{\infty}(\Omega)$ be non-negative and
$\|\mathbf{y}^{0}\|_{1}=1$. Then the solution $\mathbf{y}$ of the PDE (17)
satisfies
$\lim_{t\rightarrow\infty}\|({\color[rgb]{0,0,0}y_{2}(\cdot,t)}-y^{d})_{+}\|_{\infty}=0$.

###### Proof.

Suppose that, for the sake of contradiction, this is not true. Then, due to
the partial monotonicity property of the solution $\mathbf{y}$ stated in
Proposition IV.18, there exists a subset $\Omega_{1}\subseteq\Omega$ of
positive measure, and a parameter $\epsilon>0$, such that
${\color[rgb]{0,0,0}y_{2}(\mathbf{x},t)}-y^{d}(\mathbf{x})\geq\epsilon$ for
almost every ${\color[rgb]{0,0,0}\mathbf{x}}\in\Omega_{1}$ and all $t\geq 0$.
However, we know that

$\displaystyle y_{2}(\mathbf{x},t)$ $\displaystyle=$ $\displaystyle
y_{2}(\mathbf{x},{\color[rgb]{0,0,0}0})-\int_{{\color[rgb]{0,0,0}0}}^{t}{\color[rgb]{0,0,0}F_{2}}(y_{2}(\mathbf{x},\tau))y_{2}(\mathbf{x},\tau)d\tau$
$\displaystyle=$ $\displaystyle
y_{2}(\mathbf{x},{\color[rgb]{0,0,0}0})-\int_{{\color[rgb]{0,0,0}0}}^{t}r_{2}(y_{2}(\mathbf{x},\tau)-y^{d}(\mathbf{x}))y_{2}(\mathbf{x},\tau)d\tau$

for almost every $\mathbf{x}\in\Omega_{1}$ and for all $t\geq 0$. We know that
the function $r_{2}$ is non-zero and continuous on the open interval
$t\in(0,\infty)$. Hence, there must exist $\delta>0$ such that

$\displaystyle y_{2}(\mathbf{x},t)~{}$ $\displaystyle\leq$
$\displaystyle~{}y_{2}(\mathbf{x},0)-\int_{0}^{t}\delta
y_{2}(\mathbf{x},\tau)d\tau$ $\displaystyle\leq$
$\displaystyle~{}y_{2}(\mathbf{x},0)-\delta\int_{0}^{t}(y^{d}(\mathbf{x})+\epsilon)d\tau$
(38)

for almost every $\mathbf{x}\in\Omega_{1}$ and for all $t\geq 0$. This leads
to a contradiction. ∎

Finally, we can establish attractivity of the equilibrium point
$\mathbf{y}^{d}\in L^{\infty}(\Omega)$. Towards this end, we first prove in
the lemma below that the density of the agents in the state $Y(t)=0$ (i.e. the
state of motion, given by $y_{1}$, must converge to $0$ eventually.

######  Lemma IV.21.

Let $\mathbf{y}^{0}\in\mathbf{L}^{\infty}(\Omega)$ be non-negative and
$\|\mathbf{y}^{0}\|_{1}=1$. Then
$\lim_{t\rightarrow\infty}\|y_{1}(\cdot,t)\|=0.$ Hence,
$\lim_{t\rightarrow\infty}\|y_{2}(\cdot,t)\|=1.$

###### Proof.

Suppose that, for the sake of contradiction, this is not true. Then there
exists $\epsilon_{1}>0$ and a sequence of increasing time instants
$(t_{i})_{i=1}^{\infty}$ such that $\lim_{i\rightarrow\infty}t_{i}=\infty$ and
$\|{\color[rgb]{0,0,0}y_{1}(\cdot,t_{i})}\|_{1}\geq\epsilon_{1}$ for all
$i\in\mathbb{Z}_{+}$. From Theorem IV.16, we know that this implies that there
exist constants $\tau,\epsilon_{2},\delta>0$ such that
$y_{1}(\cdot,t)\geq\epsilon_{2}\|y_{1}(\cdot,t_{i})\|_{1}\geq\epsilon_{1}\epsilon_{2}$
for all $t{\color[rgb]{0,0,0}~{}\in~{}}[t_{i},t_{i}+\delta]$, for all
$i\in\mathbb{Z}_{+}$. Without loss of generality, we can assume that
$t_{i+1}-t_{i}>\delta$ for all $i\in\mathbb{Z}_{+}$. Let
$\Omega_{1}\subseteq\Omega$ be the subset of largest measure such that
$y^{0}_{2}(\mathbf{x})\geq y^{d}(\mathbf{x})$ for all
$\mathbf{x}\in\Omega_{1}$. Then, from the partial monotonicity property of the
solution $\mathbf{y}$ (Proposition IV.18), we have that, for each
$i\in\mathbb{Z}_{+}$,

$\displaystyle y_{2}(\mathbf{x},t_{i}+\delta)$ $\displaystyle=$ $\displaystyle
y_{2}(\mathbf{x},0)+\int_{0}^{t_{i}+\delta}F_{1}(y_{2}(\mathbf{x},\tau))y_{1}(\mathbf{x},\tau)d\tau$
$\displaystyle\geq$ $\displaystyle y_{2}(\mathbf{x},0)$
$\displaystyle+\sum_{j=1}^{i}\int_{t_{i}}^{t_{i}+\delta}r_{1}(y_{2}(\mathbf{x},\tau)-y^{d}(\mathbf{x}))y_{1}(\mathbf{x},\tau)d\tau$

for almost every $\mathbf{x}\in\Omega_{1}$. This implies that
$\lim_{i\rightarrow\infty}\|({\color[rgb]{0,0,0}y_{2}(\cdot,t_{i})}-y^{d})_{-}\|_{\infty}=0$.
However, we know that $\|\mathbf{y}(\cdot,t)\|_{1}=1$ for all $t\geq 0$. This,
along with the fact that
$\lim_{t\rightarrow\infty}\|({\color[rgb]{0,0,0}y_{2}(\cdot,t)}-y^{d})_{+}\|_{\infty}=0$
( Proposition IV.20) and the assumption that
$\|{\color[rgb]{0,0,0}y_{1}(\cdot,t_{i})}\|_{1}\geq\epsilon_{1}$ for all
$i\in\mathbb{Z}_{+}$, leads to a contradiction. ∎

Using the partial monotonicity property of solutions established in
Proposition IV.18 and the result in Lemma IV.21, we now obtain the following
global asymptotic stability result.

######  Theorem IV.22.

($L^{1}$-Global Asymptotic Stability) Let $\mathbf{y}^{0}\in
L^{\infty}(\Omega)$ be non-negative and $\|\mathbf{y}^{0}\|_{1}=1$. Then
$\lim_{t\rightarrow\infty}\|\mathbf{y}(\cdot,t)-\mathbf{y}^{d}\|_{1}=0$, and
hence the system (17) is globally asymptotically stable about the target
equilibrium distribution $\mathbf{y}^{d}$.

###### Proof.

Let $\Omega_{1}=\\{\mathbf{x}\in\Omega;~{}y^{0}_{2}(\mathbf{x})\geq
y^{d}(\mathbf{x})\\}$. Let $\Omega_{2}=\Omega-\Omega_{1}$. From  Proposition
IV.20, we know that
$\lim_{t\rightarrow\infty}\|y(\cdot,t)|_{\Omega_{1}}-y^{d}|_{\Omega_{1}}\|_{\infty}=0$,
where $\cdot|_{\Omega}$ denotes the restriction operation. This implies that
$\lim_{t\rightarrow\infty}\|y(\cdot,t)|_{\Omega_{1}}-y^{d}|_{\Omega_{1}}\|_{1}=0$.
We also know from  Proposition IV.20 that

$\lim_{t\rightarrow\infty}\int_{\Omega_{1}}\big{(}y_{2}(\mathbf{x},t)-y^{d}(\mathbf{x})\big{)}d\mathbf{x}+\int_{\Omega_{2}}\big{(}y_{2}(\mathbf{x},t)-y^{d}(\mathbf{x})\big{)}d\mathbf{x}=0$
(39)

This implies that

$\lim_{t\rightarrow\infty}\int_{\Omega_{2}}\big{(}y_{2}(\mathbf{x},t)-y^{d}(\mathbf{x})\big{)}d\mathbf{x}=0$
(40)

From Proposition IV.18, we know that $y_{2}(\cdot,t)\leq y^{d}$ almost
everywhere on $\Omega_{2}$ and for all $t\geq 0$. Hence, we conclude that

$\displaystyle\lim_{t\rightarrow\infty}\|y_{2}(\cdot,t)-y^{d}\|_{1}=\lim_{t\rightarrow\infty}\int\limits_{\Omega_{1}}|y_{2}(\mathbf{x},t)-y^{d}(\mathbf{x})|d\mathbf{x}$
$\displaystyle+\int\limits_{\Omega_{2}}|y_{2}(\mathbf{x},t)-y^{d}(\mathbf{x})|d\mathbf{x}$

From this, along with the fact that
$\lim_{t\rightarrow\infty}\|y_{1}(\cdot,t)\|_{1}=0$, we arrive at our result.
∎

##  V SIMULATIONS

In this section, we validate the control laws presented in Sections III and IV
with numerical simulations. The SDEs (III) and (13) were simulated using the
method of Wong-Zakai approximations [58]. The diffusion and reaction parameter
values used in each simulation were chosen with the goal of shortening the
duration of the simulation on a case-by-case basis. Hence, different parameter
values were chosen for each of the examples below. In practice, these
parameters would need to be chosen according to the physical constraints on
the system and the objectives of the user, such as optimizing the rate of
convergence to the target density or controlling the variance of the agent
density around the target density.

###  V-A Density Control without Agent Interactions

In this subsection, we simulate the control approach presented in Section III.
In the following example, we simulate the SDE (III) with the control laws
$v_{i}=0$ and $u_{i}(\mathbf{x})=D/y_{d}(\mathbf{x})$, where $D$ is a
diffusion coefficient. The generator of the process is given by the operator
in (8).

######  Example V.1.

Brockett integrator

In this example, we consider the case where each agent’s motion evolves
according to the Brockett integrator, which has been well-studied in the
control theory literature [14, 5]. The control vector fields for this system
are the following:

$\displaystyle X_{1}(\mathbf{x})=\frac{\partial}{\partial
x_{1}}-x_{2}\frac{\partial}{\partial
x_{3}},~{}~{}~{}~{}X_{2}(\mathbf{x})=\frac{\partial}{\partial
x_{2}}+x_{1}\frac{\partial}{\partial x_{3}}$ (41)

The Lie bracket of the two vector fields is given by

$\displaystyle[X_{1},X_{2}](\mathbf{x})=X_{1}X_{2}-X_{2}X_{1}=2\frac{\partial}{\partial
x_{3}}$ (42)

for all $\mathbf{x}\in{\color[rgb]{0,0,0}\mathbb{R}^{n}}$. Hence, we have that
${\rm
span~{}}\\{X_{1}(\mathbf{x}),X_{2}(\mathbf{x}),[X_{1},X_{2}](\mathbf{x})\\}=T_{\mathbf{x}}\mathbb{R}^{3}$
and therefore, the system is bracket generating. We define the domain
$\Omega=[0,100]^{3}$. The target distribution is given by
$y^{d}=c~{}[\sum_{i=1}^{8}\mathbf{1}_{B_{\mathbf{x}_{i}}}+0.001],$ where $c>0$
is a normalization constant that makes $y^{d}$ a probability density,
$\mathbf{1}_{S}$ denotes the indicator function of a set $S$, and
$B_{\mathbf{x}_{i}}$ denotes a ball of radius $12.5$ centered at
$\mathbf{x}_{i}$, $i=1,...,8$, defined as
$\\{\mathbf{x}_{1},...,\mathbf{x}_{8}\\}=\\{[25~{}25~{}25]^{T},[25~{}25~{}75]^{T},[25~{}75~{}75]^{T},...,[75~{}75~{}75]^{T}\\}.$

The positions of $N_{p}=10,000$ agents are generated from a stochastic
simulation of the SDE (III) and plotted at three times $t$ in Fig. 1. The
figure shows that at time $t=100~{}s$, the distribution of the swarm over the
sphere is close to the target density.  In Fig. 2, it can be seen that the
$L_{1}$ norm of the difference between the current distribution and the target
distribution decreases over time.

(a) $t=0$ s

(b) $t=10$ s

(c) $t=100$ s

Figure 1: (Brockett integrator without agent interactions) Stochastic coverage
by $N_{p}=10,000$ agents (in red) at three times $t$, following the linear
diffusion model (10). Figure 2:  (Brockett integrator without agent
interactions) Time evolution of the $L_{1}$ norm of the difference between the
target distribution and the agent distribution that evolves according to the
linear diffusion model (10).

###  V-B Density Control with Agent Interactions

In this subsection, we simulate the control approach in Section IV. The
functions $r_{i}$ in Eq. (16) are chosen to be

$r_{1}(\mathbf{x})=\begin{cases}-k\mathbf{x}~{}~{}{\rm if}~{}\mathbf{x}<0\\\
0~{}~{}~{}~{}{\rm if}~{}\mathbf{x}=0\end{cases}$ (43)
$r_{2}(\mathbf{x})=\begin{cases}k\mathbf{x}~{}~{}{\rm if}~{}\mathbf{x}>0\\\
0~{}~{}~{}~{}{\rm if}~{}\mathbf{x}=0\end{cases}$ (44)

for all $\mathbf{x}\in\Omega$, where $k$ is a positive scaling constant. Since
we simulate a finite number of agents, instead of the density $y_{1}$, the
agents use the empirical measure
$\frac{1}{N_{p}}\sum_{i=1}^{N_{p}}\delta_{\mathbf{x}_{i}(t)}$ to compute their
transition rates. However, the empirical measure is not absolutely continuous
with respect to the Riemannian volume and does not have a density. Therefore,
the agents use the regularized approximation of the measure
$\frac{1}{N_{p}}\sum_{i=1}^{N_{p}}\delta_{\mathbf{x}_{i}(t)}$, given by

$\tilde{\rho}(\mathbf{x},t)=c(\epsilon)\frac{1}{N_{p}}\sum_{i=1}^{N_{p}}K_{\epsilon}(\mathbf{x},\mathbf{x}_{i}(t))$
(45)

for all $\mathbf{x}\in M$, where the kernel function $K_{\epsilon}$ is chosen
such that $\lim_{\epsilon\rightarrow
0}c(\epsilon)K_{\epsilon}(\cdot,\mathbf{y})=\delta_{\mathbf{y}}$ for each
$\mathbf{y}\in\Omega$, and the function $c(\epsilon)$ is a normalization
parameter defined such that
$c(\epsilon)\int_{M}K_{\epsilon}(\mathbf{x},\mathbf{y})d\mathbf{x}=1$. For
each of the examples below, we will specify the kernel function used. The
positions of each agent are generated according to the SDE (13). The
transition rates $q_{i}(\mathbf{x},t)$ of the agents are defined as

$q_{i}(\mathbf{x},t)=r_{i}\Big{(}c(\epsilon)\frac{1}{N_{p}}\sum_{i=1}^{N_{p}}K_{\epsilon}(\mathbf{x},\mathbf{x}_{i}(t))-y^{d}(\mathbf{x})\Big{)}.$
(46)

######  Example V.2.

Brockett integrator

In this example, each agent moves according to Eq. (13) with the control
vector fields as defined in Example V.1. The kernel function is given by:

$K_{\epsilon}(\mathbf{x},\mathbf{y})=\begin{cases}\exp{\frac{-1}{1-(|\mathbf{x}-\mathbf{y}|/\epsilon)^{2}}}~{}{\rm~{}if~{}}|\mathbf{x}-\mathbf{y}|<\epsilon\\\
0~{}~{}{\rm otherwise}\end{cases}$ (47)

for all $\mathbf{x},\mathbf{y}\in\mathbb{R}^{3}$. The target density is set to
$y^{d}=c\sum_{i=1}^{8}\mathbf{1}_{B_{\mathbf{x}_{i}}}$, similar to the $y^{d}$
defined in Example V.1, where $c>0$ is a normalization constant. Note that in
this example, the probability density is allowed to take value equal to $0$ in
certain regions of the the domain, unlike in Example V.1. The reaction
constant in Eqs. (43)-(44) is set to $k=500$, and the parameter of the kernel
in (47) is defined as $\epsilon=5$.

(a) $t=0$ s

(b) $t=10$ s

(c) $t=100$ s

Figure 3: (Brockett integrator with agent interactions) Stochastic coverage of
$\mathbb{R}^{3}$ by $N=1,000$ agents at three times $t$, following the
semilinear PDE model (17). Blue agents are in the motion state; red agents are
in the motionless state. Figure 4:  (Brockett integrator with agent
interactions) Time evolution of the $L_{1}$ norm of the difference between the
target distribution and the agent distribution in the motionless state that
evolves according to the semilinear PDE model (17).

The positions of $N_{p}=1,000$ agents are generated from a stochastic
simulation of the SDE (13) and plotted in Fig. 3 at three times $t$. As can be
seen in this figure, at time  $t=100~{}s$ the swarm is uniformly distributed
over the sets $B_{\mathbf{x}_{i}}$ according to the target density.  Figure 4
shows that the $L_{1}$ norm of the difference between the current distribution
of agents in the motionless state and the target distribution decreases over
time. At  $t=100~{}s$, there are only $40$ agents in the state of motion. In
contrast, when using the control approach without agent interactions (Section
III), all agents are constantly in motion.  The $L_{1}$ norm at time
$t=100~{}s$ is smaller for the case with agent interactions (Fig. 4) than for
the case without agent interactions (Fig. 2).

While the interacting-agent control approach, unlike the control approach
without interactions, enables agents to stop moving once the target density is
reached (and therefore stop unnecessarily expending energy), the time until
the interacting agents converge to the target density was found to be
sensitive to the reaction constant $k$. Lower values of $k$, e.g. $k=10$,
resulted in a slower rate of agent transitions to the motionless state. On the
other hand, if the value of $k$ was chosen too large, some of the agents
prematurely transitioned to the motionless state in regions close to their
initial positions. The performance of the interacting-agent control law was
also affected by the parameter $\epsilon$. If $\epsilon$ was taken to be too
small, for example $\epsilon=0.1$, the agents did not converge to any
distribution, but instead remained in a state of motion. This can be
attributed to the fact that given a fixed value of $N_{p}$, the sum
$c(\epsilon)\frac{1}{N_{p}}\sum_{i=1}^{N_{p}}K_{\epsilon}(\mathbf{x},\mathbf{x}_{i}(t))$
becomes a less accurate approximation of the density $y(\mathbf{x},t)$ as the
value of $\epsilon$ is decreased. On the other hand, if $\epsilon$ is taken to
be too large, then the agent density converges to a regularized approximation
of the target density, rather than the target density itself. Due to space
limitations, we do not include numerical results on the effects of these
parameters here.

######  Example V.3.

Underactuated system on the sphere

In this example, we consider a system on the $2$-dimensional sphere embedded
in $\mathbb{R}^{3}$ given by
$S^{2}=\\{\mathbf{x}\in\mathbb{R}^{3};~{}\mathbf{x}^{T}\mathbf{x}=1\\}$.  We
define the following matrices $\mathbf{B}_{i}$, $i=1,2,3$:

$\displaystyle\mathbf{B}_{1}=\begin{bmatrix}0&-1&0\\\ 1&0&0\\\
0&0&0\end{bmatrix},~{}\mathbf{B}_{2}=\begin{bmatrix}0&0&1\\\ 0&0&0\\\
-1&0&0\end{bmatrix},~{}\mathbf{B}_{3}=\begin{bmatrix}0&0&0\\\ 0&0&-1\\\
0&1&0\end{bmatrix}.$

Each of the matrices defines a vector field $\tilde{X}_{i}$ on $S^{2}$ given
by

$(\tilde{X}_{i}f)(\mathbf{x})=\lim_{t\rightarrow
0}\frac{f(e^{t\mathbf{B}_{i}}\mathbf{x})-f(\mathbf{x})}{t}$ (48)

for all $\mathbf{x}\in S^{2}$ and all functions $f\in C^{\infty}(S^{2})$. We
assume that each agent can control its motion along the vector fields
${\tilde{X}_{1},\tilde{X}_{2}}$. Note that in this case, the system is
underactuated. This is because ${\rm
span}\\{\tilde{X}_{1}(\mathbf{x}),\tilde{X}_{2}(\mathbf{x}),[\tilde{X}_{1}(\mathbf{x}),\tilde{X}_{2}(\mathbf{x})]\\}=T_{\mathbf{x}}S^{2}$,
where it can be verified that
$[\tilde{X}_{1},\tilde{X}_{2}](\mathbf{x})=\tilde{X}_{3}(\mathbf{x})$ for all
$\mathbf{x}\in S^{2}$.

The kernel function is defined as

$K_{\epsilon}(\mathbf{x},\mathbf{y})=\begin{cases}\exp{\frac{-1}{1-({\rm
acos}(\mathbf{x}^{T}\mathbf{y})/\epsilon)^{2}}}~{}~{}{\rm if}~{}{\rm
acos}(\mathbf{x}^{T}\mathbf{y})<\epsilon\\\ 0~{}~{}{\rm otherwise}\end{cases}$
(49)

for all $\mathbf{x},\mathbf{y}\in S^{2}$. The target density
$y^{d}:S^{2}\rightarrow\mathbb{R}_{\geq 0}$ (with respect to the Haar measure)
is given by

$y^{d}(\mathbf{x})=\begin{cases}c~{}~{}{\rm if}~{}x^{2}_{i}\geq 0.75~{}{\rm
for}~{}i\in\\{1,2,3\\}\\\ 0~{}~{}{\rm otherwise}\end{cases}$ (50)

for all $\mathbf{x},\mathbf{y}\in S^{2}$, where $c$ is a normalization
parameter chosen such that this function integrates to $1$. We set
$\epsilon=0.1$.

The positions of $N_{p}=1,000$ agents are generated from a stochastic
simulation of the SDE (13) and plotted in Fig. 5 at three times $t$. The
target density $y^{d}$ is depicted on the surface of the sphere using a color
density plot. Blue regions are assigned a low target density of agents, while
yellow regions are assigned a high target density. The agent positions are
superimposed on the density plot to enable comparison between the actual and
target densities. Figure 5 shows that at time $t=100~{}s$, the distribution of
the swarm over the sphere is close to the target density. As for the case of
the Brockett integrator in Example V.2, only a small fraction of the swarm
($62$ agents) is in the state of motion once the swarm has converged closely
to the target density ($t=100~{}s$).

(a) $t=0$ s

(b) $t=10$ s

(c) $t=100$ s

Figure 5: ( Underactuated system on $S^{2}$ with agent interactions)
Stochastic coverage of $S^{2}$ by $N=1,000$ agents at three times $t$,
following the semilinear PDE model (17). White agents are in the motion state;
red agents are in the motionless state.

##  VI CONCLUSION

In this article, we have generalized our diffusion-based multi-agent coverage
approach to the case where the agents have nonholonomic dynamics. We
established exponential stability of the resulting Kolmogorov forward
equation, whose generator is a hypoelliptic operator. In addition, we
constructed a hybrid switching diffusion process of mean-field type such that
the probability density of the random variable that represents the
distribution of a swarm can be stabilized to a target density that is not
necessarily positive everywhere on the domain. One possible direction for
future work is to investigate the tradeoffs between control laws with and
without agent interaction.  Another is to incorporate pairwise interactions
between agents that model collision avoidance maneuvers, which would require
the inclusion of corresponding interaction terms in the PDE model. One could
also investigate the convergence of the $N$-agent system of hybrid switching
diffusion processes to the solution of the semilinear PDE.

## References

  * [1] Behçet Açıkmeşe and David S Bayard.  Markov chain approach to probabilistic guidance for swarms of autonomous agents.  Asian Journal of Control, 17(4):1105–1124, 2015.
  * [2] Andrei Agrachev, Davide Barilari, and Ugo Boscain.  A comprehensive introduction to sub-Riemannian geometry, volume 181\.  Cambridge University Press, 2019.
  * [3] Andrei Agrachev, Ugo Boscain, Jean-Paul Gauthier, and Francesco Rossi.  The intrinsic hypoelliptic Laplacian and its heat kernel on unimodular lie groups.  Journal of Functional Analysis, 256(8):2621–2655, 2009.
  * [4] Andrei A Agrachev and Yuri Sachkov.  Control theory from the geometric viewpoint, volume 87.  Springer Science & Business Media, 2013.
  * [5] Alessandro Astolfi.  Discontinuous control of the Brockett integrator.  European Journal of Control, 4(1):49–63, 1998.
  * [6] Fabrice Baudoin and Michel Bonnefont.  Log-Sobolev inequalities for subelliptic operators satisfying a generalized curvature dimension inequality.  Journal of Functional Analysis, 262(6):2646–2676, 2012.
  * [7] Jean-David Benamou and Yann Brenier.  A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem.  Numerische Mathematik, 84(3):375–393, 2000.
  * [8] Spring Berman, Ádám Halász, M Ani Hsieh, and Vijay Kumar.  Optimized stochastic policies for task allocation in swarms of robots.  IEEE Transactions on Robotics, 25(4):927–937, 2009.
  * [9] Subhrajit Bhattacharya, Robert Ghrist, and Vijay Kumar.  Multi-robot coverage and exploration on Riemannian manifolds with boundaries.  The International Journal of Robotics Research, 33(1):113–137, 2014\.
  * [10] Shiba Biswal, Karthik Elamvazhuthi, and Spring Berman.  Stabilization of nonlinear discrete-time systems to target measures using stochastic feedback laws.  IEEE Transactions on Automatic Control, 66(5):1957–1972, 2020.
  * [11] Shiba Biswal, Karthik Elamvazhuthi, and Spring Berman.  Stabilization of multi-agent systems to target distributions using local interactions.  IFAC-PapersOnLine, 54(9):60–65, 2021.
  * [12] Benoît Bonnet and Francesco Rossi.  The Pontryagin maximum principle in the Wasserstein space.  Calculus of Variations and Partial Differential Equations, 58(1):11, 2019.
  * [13] Marco Bramanti.  An invitation to hypoelliptic operators and Hörmander’s vector fields.  Springer, 2014.
  * [14] Roger W Brockett.  Asymptotic stability and feedback stabilization.  R. W. Brockett, R. S. Millman and H. J. Sussmann, eds., Differential geometric control theory, 27(1):181–191, 1983.
  * [15] Francesco Bullo and Andrew D Lewis.  Geometric control of mechanical systems: modeling, analysis, and design for simple mechanical control systems, volume 49.  Springer, 2019.
  * [16] Mourad Choulli and Laurent Kayser.  Gaussian lower bound for the Neumann Green function of a general parabolic operator.  Positivity, 19(3):625–646, 2015.
  * [17] Jorge Cortes, Sonia Martinez, Timur Karatas, and Francesco Bullo.  Coverage control for mobile sensing networks.  IEEE Transactions on robotics and Automation, 20(2):243–255, 2004\.
  * [18] Michel Duprez and Antoine Perasso.  Criterion of positivity for semilinear problems with applications in biology.  Positivity, 21(4):1383–1392, 2017.
  * [19] Karthik Elamvazhuthi, Chase Adams, and Spring Berman.  Coverage and field estimation on bounded domains by diffusive swarms.  In IEEE Conference on Decision and Control (CDC), pages 2867–2874, 2016.
  * [20] Karthik Elamvazhuthi and Spring Berman.  Nonlinear generalizations of diffusion-based coverage by robotic swarms.  In 2018 IEEE Conference on Decision and Control (CDC), pages 1341–1346. IEEE, 2018.
  * [21] Karthik Elamvazhuthi and Spring Berman.  Mean-field models in swarm robotics: A survey.  Bioinspiration & Biomimetics, 15(1):015001, 2019.
  * [22] Karthik Elamvazhuthi, Shiba Biswal, and Spring Berman.  Mean-field stabilization of robotic swarms to probability distributions with disconnected supports.  In 2018 Annual American Control Conference (ACC), pages 885–892. IEEE, 2018.
  * [23] Karthik Elamvazhuthi, Hendrik Kuiper, and Spring Berman.  PDE-based optimization for stochastic mapping and coverage strategies using robotic ensembles.  Automatica, 95:356–367, 2018.
  * [24] Karthik Elamvazhuthi, Hendrik Kuiper, Matthias Kawski, and Spring Berman.  Bilinear controllability of a class of advection–diffusion–reaction systems.  IEEE Transactions on Automatic Control, 64(6):2282–2297, 2018.
  * [25] Klaus-Jochen Engel and Rainer Nagel.  One-Parameter Semigroups for Linear Evolution Equations, volume 194\.  Springer Science & Business Media, 2000.
  * [26] Utku Eren and Behçet Açıkmeşe.  Velocity field generation for density control of swarms using heat equation and smoothing kernels.  IFAC-PapersOnLine, 50(1):9405–9411, 2017.
  * [27] Lawrence C Evans.  Partial differential equations.  Graduate Studies in Mathematics, 19, 1998.
  * [28] Massimo Fornasier and Francesco Solombrino.  Mean-field optimal control.  ESAIM: Control, Optimisation and Calculus of Variations, 20(4):1123–1152, 2014.
  * [29] Nicola Garofalo and Duy-Minh Nhieu.  Isoperimetric and Sobolev inequalities for Carnot-Carathéodory spaces and the existence of minimal surfaces.  Communications on Pure and Applied Mathematics, 49(10):1081–1144, 1996.
  * [30] Thanh-Trung Han and Shuzhi Sam Ge.  Styled-velocity flocking of autonomous vehicles: A systematic design.  IEEE Transactions on Automatic Control, 60(8), 2015.
  * [31] Minyi Huang, Roland P Malhamé, Peter E Caines, et al.  Large population stochastic dynamic games: closed-loop mckean-vlasov systems and the nash certainty equivalence principle.  Communications in Information & Systems, 6(3):221–252, 2006.
  * [32] David Jerison and Carlos E. Kenig.  The functional calculus for the Laplacian on Lipschitz domains.  Journées Equations aux Dérivées Partielles, pages 1–10, 1989.
  * [33] David S Jerison and Antonio Sáchez-CAlle.  Estimates for the heat kernel for a sum of squares of vector fields.  Indiana University mathematics journal, 35(4):835–854, 1986.
  * [34] Ioannis Karatzas and Steven E Shreve.  Brownian motion.  In Brownian Motion and Stochastic Calculus, pages 47–127. Springer, 1998.
  * [35] Hassan K Khalil and Jessy W Grizzle.  Nonlinear systems, volume 3.  Prentice hall Upper Saddle River, NJ, 2002.
  * [36] Vassili N Kolokoltsov.  Nonlinear Markov processes and kinetic equations, volume 182.  Cambridge University Press, 2010.
  * [37] Vishaal Krishnan and Sonia Martínez.  Distributed optimal transport for the deployment of swarms.  In 2018 IEEE Conference on Decision and Control (CDC), pages 4583–4588. IEEE, 2018.
  * [38] Jean-Michel Lasry and Pierre-Louis Lions.  Mean field games.  Japanese journal of mathematics, 2(1):229–260, 2007.
  * [39] John M Lee.  Smooth manifolds.  In Introduction to Smooth Manifolds, pages 1–31. Springer, 2013\.
  * [40] John M Lee.  Introduction to Riemannian manifolds.  Springer, 2018.
  * [41] Taeyoung Lee.  Bayesian attitude estimation with the matrix fisher distribution on so (3).  IEEE Transactions on Automatic Control, 63(10):3377–3392, 2018\.
  * [42] Peter Li, Shing Tung Yau, et al.  On the parabolic kernel of the schrödinger operator.  Acta Mathematica, 156:153–201, 1986.
  * [43] Alessandra Lunardi.  Analytic Semigroups and Optimal Regularity in Parabolic Problems.  Springer Science & Business Media, 2012.
  * [44] T William Mather and M Ani Hsieh.  Distributed robot ensemble control for deployment to multiple sites.  Robotics: Science and Systems VII, 2011.
  * [45] Siddharth Mayya, Sean Wilson, and Magnus Egerstedt.  Closed-loop task allocation in robot swarms using inter-robot encounters.  Swarm Intelligence, 13(2):115–143, 2019.
  * [46] Henry P McKean.  Propagation of chaos for a class of non-linear parabolic equations.  Stochastic Differential Equations (Lecture Series in Differential Equations, Session 7, Catholic Univ., 1967), pages 41–57, 1967\.
  * [47] Alexandre R Mesquita and Joao P Hespanha.  Jump control of probability densities with applications to autonomous vehicle motion.  IEEE Transactions on Automatic Control, 57(10):2588–2598, 2012\.
  * [48] Alexandre R Mesquita, João P Hespanha, and Karl Åström.  Optimotaxis: A stochastic multi-agent optimization procedure with point measurements.  In International workshop on hybrid systems: Computation and control, pages 358–371. Springer, 2008.
  * [49] Duy-Minh Nhieu.  The Neumann problem for sub-Laplacians on carnot groups and the extension theorem for Sobolev spaces.  Annali di Matematica Pura ed Applicata, 180(1):1–25, 2001.
  * [50] El-Maati Ouhabaz.  Analysis of Heat Equations on Domains (LMS-31).  Princeton University Press, 2009.
  * [51] Benedetto Piccoli, Francesco Rossi, and Emmanuel Trélat.  Control to flocking of the kinetic Cucker–Smale model.  SIAM Journal on Mathematical Analysis, 47(6):4685–4719, 2015.
  * [52] Andrey Pilipenko.  An Introduction to Stochastic Differential Equations with Reflection, volume 1.  Universitätsverlag Potsdam, 2014.
  * [53] Hannes Risken.  Fokker-Planck equation.  In The Fokker-Planck Equation, pages 63–95. Springer, 1996.
  * [54] Francesco Rossi.  Large time behavior for the heat equation on Carnot groups.  Nonlinear Differential Equations and Applications NoDEA, 20(3):1393–1407, 2013.
  * [55] Tomáš Roubíček.  Nonlinear partial differential equations with applications, volume 153.  Springer Science & Business Media, 2013.
  * [56] Alain Sarlette and Rodolphe Sepulchre.  Consensus optimization on manifolds.  SIAM Journal on Control and Optimization, 48(1):56–76, 2009.
  * [57] Hiroshi Tanaka.  Stochastic differential equations with reflecting boundary condition in convex regions.  Hiroshima Mathematical Journal, 9(1):163–177, 1979.
  * [58] Krystyna Twardowska.  Wong-Zakai approximations for stochastic differential equations.  Acta Applicandae Mathematica, 43(3):317–359, 1996.
  * [59] Gabriele Valentini.  Achieving consensus in robot swarms.  2017\.
  * [60] Nicholas T Varopoulos, Laurent Saloff-Coste, and Thierry Coulhon.  Analysis and geometry on groups, volume 100.  Cambridge university press, 2008.
  * [61] Sean Wilson, Theodore P Pavlic, Ganesh P Kumar, Aurélie Buffin, Stephen C Pratt, and Spring Berman.  Design of ant-inspired stochastic control policies for collective transport by robotic swarms.  Swarm Intelligence, 8(4):303–327, 2014.
  * [62] George Yin and Chao Zhu.  Hybrid Switching Diffusions: Properties and Applications, volume 63.  Springer New York, 2010.

|  Karthik Elamvazhuthi received the B.Tech. degree in mechatronics
engineering from the Manipal Institute of Technology, Manipal, India in 2011.
He received the M.S. and Ph.D. degrees in mechanical engineering from Arizona
State University, Tempe, AZ, in 2014 and 2019, respectively. He was a CAM
Assistant Adjunct Professor at the Department of Mathematics, University of
California, Los Angeles from 2019–2022. Currently, he is a postdoctoral
scholar at the Department of Mechanical Engineering, University of California,
Riverside. His research interests include modeling and control of robotic
swarms using methods from partial differential equations and stochastic
processes.
---|---
|  Spring Berman (M’07) received the M.S.E. and Ph.D. degrees in Mechanical
Engineering and Applied Mechanics from the University of Pennsylvania in 2008
and 2010, respectively. From 2010 to 2012, she was a Postdoctoral Fellow in
Computer Science at Harvard University. She was then an Assistant Professor of
Mechanical and Aerospace Engineering at Arizona State University, where she
has been an Associate Professor since 2018. Her research focuses on the
analysis of behaviors in biological and engineered collectives and the
synthesis of control strategies for robotic swarms.
---|---