Controllability to Equilibria of the 1-D Fokker-Planck Equation with Zero-Flux Boundary Condition

TL;DR

This paper demonstrates that under specific regularity and boundary conditions, the probability distribution of a robotic swarm governed by a 1-D Fokker-Planck equation can be exactly controlled from any initial to a desired target distribution using bounded drift fields.

Contribution

It introduces a control method for the 1-D Fokker-Planck equation with zero-flux boundary conditions, allowing exact distribution steering with bounded controls, which was not previously established.

Findings

Any initial distribution can be transported to a target distribution under certain regularity conditions.

Exact control is possible with bounded drift vector fields.

The proof employs classical linear semigroup theory.

Abstract

We consider the problem of controlling the spatiotemporal probability distribution of a robotic swarm that evolves according to a reflected diffusion process, using the space- and time-dependent drift vector field parameter as the control variable. In contrast to previous work on control of the Fokker-Planck equation, a zero-flux boundary condition is imposed on the partial differential equation that governs the swarm probability distribution, and only bounded vector fields are considered to be admissible as control parameters. Under these constraints, we show that any initial probability distribution can be transported to a target probability distribution under certain assumptions on the regularity of the target distribution. In particular, we show that if the target distribution is (essentially) bounded, has bounded first-order and second-order partial derivatives, and is bounded from below by a strictly positive constant, then this distribution can be reached exactly using a drift vector field that is bounded in space and time. Our proof is constructive and based on classical linear semigroup theoretic concepts.

Equations

$$dZ(t)$$

$$Z(0)$$

$$y_t=y_{xx}-(vy{)}_x$$

$$⟨p,q{⟩}_2={\int }_0^{1}p(x)q(x)dx$$

$$\Vert p{\Vert }_2=⟨p,p{⟩}_2^{1/2}$$

$$\|p\|_{H^{1}}=\Big{(}\|p\|^{2}_{2}+\|p_{x}\|^{2}_{2}\Big{)}^{1/2}$$

$$\|p\|_{H^{2}}=\Big{(}\|p\|^{2}_{H^{1}}+\|p_{xx}\|^{2}_{2}\Big{)}^{1/2}$$

$$\Vert z{\Vert }_{\infty }=esssu{p}_{x\in (0,1)}|z(x)|,$$

$${\int }_0^1|p(x){|}^{2}a(x)dx<\infty .$$

$$⟨p,q{⟩}_a={\int }_0^{1}p(x)q(x)a(x)dx$$

$$W^{1,\infty}(0,1)=\big{\{}z\in L^{\infty}(0,1):z_{x}\in L^{\infty}(0,1)\big{\}},$$

$$W^{2,\infty}(0,1)=\big{\{}z\in W^{1,\infty}(0,1):z_{xx}\in L^{\infty}(0,1)\big{\}}.$$

$$\Vert u{\Vert }_{C([0,T];L^2(0,1))}:=\underset{0\le t\le T}{\max }\Vert u(t){\Vert }_2<\infty .$$

$$\dot{u}(t)=Au(t);u(0)=u_0\in L^2(0,1)$$

$$\dot{u}(t)=A(t)u(t);u(0)=u_0\in X$$

$$u(t)=u_0+\sum _{i\in {\mathbb{Z}}_{+}}{A}_i{\int }_{min\{t,a_i\}}^{min\{t,a_{i+1}\}}u(s)ds$$

$$\Vert y(\cdot ,t)-f{\Vert }_2\le M{e}^{-\lambda t}\Vert y_0-f{\Vert }_2$$

$$y_t=a{y}_{xx}$$

$$A_{a}u=(au{)}_{xx}$$

$$y_t=(ay{)}_{xx}$$

$$u_t=a{u}_{xx}$$

$$\Vert y(\cdot ,t)-f{\Vert }_2$$

$$\Vert A_{a}y(\cdot ,t){\Vert }_2$$

$$\Vert z{\Vert }_g=\Vert z{\Vert }_2+\Vert A_{a}z{\Vert }_2$$

$$\Vert y(\cdot ,t)-f{\Vert }_g\le M_1{e}^{-\lambda t}$$

$$B_{f}u=u_{xx}-(\frac{f_x}{f}u{)}_x$$

$$y(x,t)={\int }_0^{1}K(t,x,z)y_0(x)dz$$

$$K(t,x,z)\ge \frac{1}{(4\pi t{)}^{1/2}}exp(\frac{-(x-z{)}^2}{4t})$$

$$v(\cdot ,t)=\frac{y_x}{y}-\alpha m\frac{(ay{)}_x}{y}$$

$${\tilde{y}}_t=\alpha m(a\tilde{y}{)}_{xx}$$

$$|y(\cdot ,\sum _{n=1}^m\frac{1}{n^2})-f{\Vert }_{L^2(\Omega )}\le M_0{e}^{-\alpha \lambda {\sum }_{n=1}^m\frac{n}{n^2}}$$

Full text

I INTRODUCTION

In recent years, there has been much work on the modeling and control of swarms of homogeneous agents using mean-field models. These mean-field models are typically defined by a system of partial differential equations (PDEs) that describe how an initial probability measure is pushed forward under the action of an ordinary differential equation (ODE) or a stochastic differential equation (SDE). In this context, spatial information on the agent positions is modeled using probability measures, and the mean-field control problem is to design parameters of the system of PDEs so that these probability measures evolve in a desirable manner.

This perspective has led to a number of works on optimal control of PDEs with the goal of optimizing swarm behavior. In [22], the authors use a maximum principle for control of infinite-dimensional systems to design optimal switching parameters that achieve target swarm densities. The work in [2] addresses the problem of optimal control of the Fokker-Planck equation to ensure that its solution tracks a predefined time-dependent reference density. A similar approach is used for mean-field games and mean-field type controls in [5, 17], which consider controlled versions of Vlasov-Mckean type SDEs whose coefficients are coupled to the distribution of the stochastic process. While most prior work has considered the case where there is noise in the agent dynamics, there has also been some recent work on the noise-less case where each agent evolves according to an ODE [7, 16] rather than an SDE.

The literature on control of mean-field models has mostly focused on the synthesis of optimal control strategies. However, there has been some work on questions of stabilization and controllability of such models. For example, the design of output feedback laws for designing globally stable invariant distributions of stochastic processes was considered in [21]. Control of swarm protocols governed by the kinetic Cucker-Smale model was addressed in [24] to produce flocking behavior. The Benamou-Brenier fluid dynamic formulation of optimal transport problems [4] can also be interpreted as a problem of optimal control and controllability of the continuity equation. See also the work by Brockett [9] on related problems on the control of Liouville equations. Closer to the work in this paper are the studies [6, 12, 26] on controllability of Fokker-Planck equations. In [6] and [12], it was proven that solutions of the Fokker-Planck equation evolving on $\mathbb{R}^{n}$ can be controlled to a large class of target distributions. This result has been extended to the case where the agent dynamics are governed by general linear SDEs and the initial and target distributions are Gaussian [10, 11]. Poretta [26] considered the problem of controllability of the Fokker-Planck equation when the solutions evolve on a torus, along with the well-posedness of an associated mean-field game problem.

In this work, we consider controllability of the Fokker-Planck equation on a bounded one-dimensional domain with zero-flux boundary condition. This problem is of practical importance for swarm robotic applications in which the agents’ spatiotemporal behavior can be modeled by the Fokker-Planck equation and the agents are constrained to evolve in a bounded domain [13, 15, 22, 27]. Specifically, we prove (see Theorem IV.10) that when the target probability density is sufficiently regular, it can be reached from any square-integrable initial probability distribution in a given finite time using bounded control inputs. Our arguments are entirely based on linear operator semigroup theoretic concepts and make a straightforward exploitation of the fact that the exponential convergence rate to the target equilibrium can be increased arbitrarily using an appropriate density-feedback law. This approach differs from those in [6] and [12], which are based on probabilistic and stochastic control theoretic concepts, and the approach in [26], which uses observability inequality type arguments. Although it might be possible to adapt these methods for our scenario, the constraint of the vector field being bounded is not imposed in these works, which is more relevant for practical scenarios of interest to us. On the other hand, unlike the works [12, 6, 26], we do not address any issue regarding the optimality of the control laws. We note that although we restrict our analysis in this paper to the case of 1-D domains for the sake of simplicity, our approach has natural extensions to the more practical multi-dimensional case, which will be the subject of our future work.

II PROBLEM FORMULATION

Consider a swarm of $n$ agents deployed on the one-dimensional domain $[0,1]$. The position of each agent, indexed by $i\in\{1,2,...,n\}$, evolves according to a stochastic process $Z_{i}(t)\in[0,1]$, where $t$ denotes time. Since the random variables that correspond to the dynamics of each agent are independent and identically distributed, we can drop the subscript $i$ and define the problem in terms of a single stochastic process, $Z(t)\in[0,1]$. The deterministic motion of each agent is defined by a vector field $v(x,t)\in\mathbb{R}$, where $x\in[0,1]$. This motion is perturbed by the Wiener process $W(t)$, which models noise. This process can be a model for stochasticity arising from inherent sensor and actuator noise. Alternatively, noise could be actively programmed into the agents’ motion to implement more exploratory agent behaviors and to take advantage of the smoothening effect of the process on the agents’ probability densities. Given the parameter $v(x,t)$, the stochastic process $Z(t)$ satisfies the following SDE [28]:

[TABLE]

where $d\psi(t)\in\mathbb{R}$ is called the reflecting function, a stochastic process that constrains $Z(t)$ to the domain $[0,1]$.

Problem II.1**.**

Given $T>0$ and $f:[0,1]\rightarrow\mathbb{R}^{+}$ such that $\int_{0}^{1}f(x)dx=1$, determine if there exists a feedback control law $v:[0,1]\times[0,T]\rightarrow\mathbb{R}$ such that the process (1) satisfies $\mathbb{P}(Z(T)\in\Gamma)=\int_{\Gamma}f(x)dx$ for each Borel subset $\Gamma\subset[0,1]$.

The Kolmogorov forward equation corresponding to the SDE (1) is given by [25]:

[TABLE]

The solution $y(x,t)$ of this equation represents the probability density of a single agent occupying position $x\in[0,1]$ at time $t$, or alternatively, the density of a population of agents at this position and time. The PDE (II) is related to the SDE (1) by the relation $\mathbb{P}(Z(t)\in\Gamma)=\int_{\Gamma}y(x,t)dx$ for all $t\in[0,T]$ and all measurable $\Gamma$. Problem II.1 can be reframed in terms of equation (II) as a PDE controllability problem as follows:

Problem II.2**.**

Given $T>0$, $y_{0}:[0,1]\rightarrow\mathbb{R}^{+}$, and $f:[0,1]\rightarrow\mathbb{R}^{+}$ such that $\int_{0}^{1}y_{0}(x)dx=\int_{0}^{1}f(x)dx=1$, determine whether there exists a space- and time-dependent parameter $v:[0,1]\times[0,T]\rightarrow\mathbb{R}$ such that the solution $y$ of the PDE (II) satisfies $y(\cdot,T)=f$.

III PRELIMINARIES AND NOTATION

We define $L^{2}(0,1)$ as the Hilbert space of square-integrable real-valued functions over the unit interval, $(0,1)\subset\mathbb{R}$. The Hilbertian structure of $L^{2}(0,1)$ is induced by the standard inner product, $\langle\cdot,\cdot\rangle_{2}:L^{2}(0,1)\times L^{2}(0,1)\rightarrow\mathbb{R}$, given by:

[TABLE]

for each $p,q\in L^{2}(0,1)$. The norm $\|\cdot\|_{2}$ on the space $L^{2}(0,1)$ is defined as

[TABLE]

for all $p\in L^{2}(0,1)$. For a function $r\in L^{2}(0,1)$ and a given constant $c$, we write $r\geq c$ to imply that $r(x)\geq c$ for almost every $x\in(0,1)$.

We define the Sobolev space H^{1}(0,1)=\big{\{}z\in L^{2}(0,1):z_{x}\in L^{2}(0,1)\big{\}}. We equip this space with the usual Sobolev norm $\|\cdot\|_{H^{1}}$, given by

[TABLE]

for each $p\in H^{1}(0,1)$. We will also need the space H^{2}(0,1)=\big{\{}z\in H^{1}(0,1):z_{xx}\in L^{2}(0,1)\big{\}}, which will be equipped with the norm $\|\cdot\|_{H^{2}}$ given by

[TABLE]

for each $p\in H^{1}(0,1)$. We define $L^{\infty}(0,1)$ as the space of essentially bounded measurable functions on $(0,1)$. The space $L^{\infty}(0,1)$ is equipped with the norm

[TABLE]

where $ess~{}sup_{x\in(0,1)}(\cdot)$ denotes the essential supremum attained by its argument over the interval $(0,1)$. For a given $a\in L^{\infty}(0,1)$, $L^{2}_{a}(0,1)$ will refer to the set of all functions $p$ such that

[TABLE]

The space $L^{2}_{a}(0,1)$ is a Hilbert space with respect to the weighted inner product $\langle\cdot,\cdot\rangle_{a}:L^{2}_{a}(0,1)\times L^{2}_{a}(0,1)\rightarrow\mathbb{R}$, given by:

[TABLE]

for each $f,g\in L^{2}(0,1)$. We also define the Sobolev spaces

[TABLE]

[TABLE]

The space $C([0,T];L^{2}(0,1))$ consists of all continuous functions $u:[0,T]\rightarrow L^{2}(0,1)$ for which

[TABLE]

We will need an appropriate notion of a solution of the PDE (II). Toward this end, let $A$ be a closed linear operator on a Hilbert space $X$ with domain $\mathcal{D}(A)$. For a given time $T>0$, a mild solution of the ODE

[TABLE]

is a function $u\in C([0,T];L^{2}(0,1))$ 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 is given by a semigroup of linear operators, $(T(t))_{t\geq 0}$, that are generated by the operator $A$. That is, the solution of the above ODE is given by $u(t)=T(t)u_{0}$ for each $t\in[0,T]$.

The differential equations that we analyze in this paper will be non-autonomous in general. Hence, we must adapt the notion of a mild solution to these types of equations. Let $A_{i}$ be a closed linear operator with domain $\mathcal{D}(A_{i})$ for each $i\in\mathbb{Z}_{+}$. For a certain time interval $[0,T]$, a piecewise constant family of operators is given by a map, $t\mapsto A(t)$, for which there exists a partition of $[0,T]=\cup_{i\in\mathbb{Z}_{+}}[a_{i},a_{i+1})$ such that $a_{i}\leq a_{i+1}$ for each $i\in\mathbb{Z}_{+}$ and $A(t)=A_{i}$ for each $t\in[a_{i},a_{i+1})$. Then a mild solution of the ODE

[TABLE]

is a function $u\in C([0,T];X)$ such that

[TABLE]

for each $t\in[0,T]$. There is in fact a more general notion of mild solutions that arise from two-parameter semigroups of operators generated by time-varying linear operators. However, the definition in (12) will be sufficient for our purposes, since one can construct solutions of the ODE (11) by treating it as an autonomous system in each time interval $[a_{i},a_{i+1})$ and patching these solutions together to obtain the solution $u$. Note that the mild solution is defined with respect to an operator $A$ or collection of operators $A(t)$; when we refer to such a solution, the associated operator(s) will be clear from the context.

IV CONTROLLABILITY ANALYSIS

In this section, we prove our main result, Theorem IV.10, as a solution to Problem II.2. However, we first note the following result on exponential stabilizability of equilibrium distributions which is also useful from the point of view of multi-agent control problems and gives an ’approximate’ result to the controllability problem II.2. Particularly, this theorem gives a candidate time-independent vector field if one desires convergence to a given target distribution at an exponential rate. This is similar to our previous work [13], where we used a spatially inhomogenous diffusion coefficient to stabilize desired probability densities. However, one needs to assume higher regularity on the target distributions in the following scenario than the work in [13].

Theorem IV.1**.**

Let $y_{0}\in L^{2}(0,1)$ and $f\in W^{1,\infty}(0,1)$ such that $y_{0}\geq 0$, $f\geq k>0$ for some strictly positive constant $k$ and $\int_{0}^{1}y_{0}(x)dx=\int_{0}^{1}f(x)dx=1$. Suppose $v(x,t)=f_{x}(x)/f(x)$ for all $t\in[0,\infty)$ and almost every $x\in(0,1)$ in the PDE, (II), then a unique mild solution $y\in C([0,\infty);L^{2}(0,1))$ of the PDE exists and satisfies the estimate

[TABLE]

for some positive constants, $M\geq 1$, $\lambda$ and all $t\in[0,\infty)$.

The above result is a straightforward corollary of a result well-known in mathematics and physics literature; namely, that if $v$ is the gradient of a potential function $\phi$, then $e^{\phi}$ is the unique stationary distribution of the stochastic process (1). From this, if one observes that $f_{x}/f=(ln~{}f)_{x}$, then it follows that any weakly differentiable function with a strictly positive lower bound can be written as an exponential distribution, $f=e^{(ln~{}f)}$. Hence, any function, $f$, that is at least once weakly differentiable can be then designed to be the equilibrium distribution of the Fokker-Planck equation by choosing the vector field to be the gradient of the potential function, $lnf$. See remarks in [20], [1] and [3]. There are multiple ways to establish the stability result. One such approach is using optimal transport methods. Using this approach the above result follows when certain convexity conditions on $f$ are satisfied and the PDE, (II), is framed as a gradient flow for an appropriate energy functional on the 2-Wasserstein space [29]. This is shown in [20] and [1]. For more general $f$ as stated in the above theorem, the result follows using classical functional analytic methods, which might not sufficient when the domain of the PDE is unbounded or infinite dimensional as in the works [20] and [1] respectively. Particularly, the operator $(\cdot)_{xx}-(\cdot~{}\frac{f_{x}}{f})_{x}$ with the zero-flux boundary condition, is a self-adjoint, negative semi-definite, closed and densely defined operator on the weighted space $L^{2}_{f}(0,1)$. $L^{2}_{f}(0,1)$ is isomorphic to the space $L^{2}(0,1)$ due to the upper and lower bounds on $f$. Hence, the spectrum of the generator, $(\cdot)_{xx}-(\frac{f_{x}}{f})_{x}$, lies on the real line and is bounded above by the principal eigenvalue. [math] is an eigenvalue that has a positive eigenvector, $f$. This implies that [math] is that principal eigenvalue. Principal eigenvalues of elliptic operators are simple and therefore the the rest of the spectrum lies to the left of some negative number, $-\lambda$. From this the result on expoenential stability of the steady state solution $f$ follows. See also [3] for more details to establish the spectral gap, $\lambda$, using Poincare inequalities.

Next, we collect some preliminary results that will be needed to prove our main theorem.

Theorem IV.2**.**

Consider the PDE,

[TABLE]

Let $a\in W^{2,\infty}(0,1)$ be such that $a\geq e$ for some strictly positive constant $e$. Let $y_{0}\in L^{2}(0,1)$. Then a unique mild solution $y\in C([0,T];L^{2}(0,1))$ of the above PDE exists. Additionally, if there exists a positive constant, $c>0$, such that $y_{0}\geq c$, then $y(\cdot,t)\geq c$ for each $t\in[0,\infty)$.

Proof.

See appendix. ∎

Now, we make some definitions which will be used subsequently. Given $a\in W^{2,\infty}(0,1)$ such that $a\geq e$ for some strictly positive constant $e$ we define the operator, $A_{a}:\mathcal{D}(A_{a})\rightarrow L^{2}(0,1)$, by

[TABLE]

for each $u\in{D}(A_{a})=\{w\in H^{2}(0,1);(aw_{x})(0)=(aw)_{x}(1)=0\}$.

Corollary IV.3**.**

Let $a\in W^{2,\infty}(0,1)$ be such that $a\geq e$ for some strictly positive constant $e$. Consider the PDE,

[TABLE]

Let $y_{0}\in L^{2}({\Omega})$. Then a unique mild solution $y\in C([0,T];L^{2}(0,1))$ of the above PDE exists. Additionally, if there exists a positive constant, $c>0$, such that $y_{0}\geq c$, then $y(\cdot,t)\geq d$ for each $t\in[0,\infty)$, for some strictly positive constant $d$.

Proof.

Both the existence of mild solutions of the PDE, (LABEL:eq:cllp1), and the lower bound on solutions follow from theorem IV.2 by noting that the coordinate transformation, $u=ay$, transforms the PDE (LABEL:eq:cllp1), to the PDE,

[TABLE]

∎

In the following lemma, we will establish some estimates on the rate of convergence of the solution, $y$, of the PDE (LABEL:eq:cllp1), assuming the initial condition is regular enough.

Lemma IV.4**.**

Let $y_{0}\in\mathcal{D}(A_{a})$ be such that $y_{0}\geq 0$. Additionally, assume $a=1/f$ , where $f\in W^{2,\infty}(0,1)$ such that $f\geq c>0$ for some constant $c$, and $\int_{0}^{1}f(x)dx=\int_{0}^{1}y_{0}(x)dx$. Then the mild solution, $y\in C([0,\infty);L^{2}(0,1))$, of the PDE, (LABEL:eq:cllp1), satisfies $y(\cdot,t)\in\mathcal{D}(A_{a})$ for each $t\in[0,\infty)$. Moreover, the following estimates hold:

[TABLE]

for some strictly positive constants $M_{0}$, $M_{1}$ and $\lambda$.

Proof.

The first estimate is just a restatement of [13][Theorem IV.4], where it was shown that $f$ is the eigenvector of $A_{a}$ corresponding to simple principal eigenvalue [math], and the rest of spectrum is in the left-half complex plane, left to some negative number $-\lambda$. For this estimate, the assumption of the regularity of the initial condition is not required.

For the second estimate, we consider the space $(\mathcal{D}(A)_{a},\|\cdot\|_{g})$ where $\|\cdot\|_{g}$ is the graph norm given by

[TABLE]

for each $z\in\mathcal{D}(A_{a})$. $A_{a}$ generates a semigroup of operators, $(T_{1}(t))_{t\geq 0}$, on the space, $(\mathcal{D}(A_{a}),\|\cdot\|_{g})$. Let $(T(t))_{t\geq 0}$ be the semigroup of operators generated on the space, $L^{2}(0,1)$ by $A_{a}$. $T_{1}(t)$ is the restriction of the operator, $T(t)$, for each $t\in[0,\infty)$, on the space $(\mathcal{D}(A_{a}),\|\cdot\|_{g})$. $T_{1}(t)=(I-A_{a})^{-1}T(t)(I-A_{a})$ for each $t\in[0,\infty)$. It follows that $T_{1}(t)$ and $T(t)$ are similar for each $t\in[0,\infty)$. Hence, the spectrum and growth bounds of $(T_{1}(t))_{t\geq 0}$ and $(T(t))_{t\geq 0}$ are the same. Hence, it follows that we have

[TABLE]

for some $M_{1}>0$ and for all $t\in[0,\infty)$. This implies the estimate, (19), since $A_{a}f=0$. ∎

Controllability of the PDE, (II), will initially be established with some assumptions on regularity conditions and lower bounds satisfied by the initial conditions (lemma IV.7). The next two results will help us relax these assumptions further ahead in the main theorem IV.10.

Theorem IV.5**.**

Let $y_{0}\in L^{2}(0,1)$ and $f\in W^{2,\infty}(0,1)$ be such that $f\geq k>0$ for some positive constant $k$. Suppose $v(x,t)=f_{x}(x)/f(x)$ for all $t\in[0,\infty)$ and almost every $x\in(0,1)$ in the PDE, (II), If there exists a positive constant, $c>0$, such that $y_{0}\geq c$, then the unique mild solution of the PDE satisfies the estimate $y(\cdot,t)\geq d$ for each $t\in[0,\infty)$, and for some positive constant, $d$.

Moreover, we have that $\mathcal{D}(B_{f})=\mathcal{D}(A_{a})$, whenever $a=1/f$ and $A_{a}$ is the operator defined in equation (15) and $B_{f}:\mathcal{D}(B_{f})\rightarrow L^{2}(0,1)$ is the operator given by

[TABLE]

for each $u\in{D}(B_{f})=\{w\in H^{2}(0,1);(w_{x}-\frac{f_{x}}{f}w)(0)=(w_{x}-\frac{f_{x}}{f}w)(1)=0\}$.

Lemma IV.6**.**

Consider the heat equation with Neumann boundary condition, that is, $v\equiv 0$ in (II). Let $y_{0}\in L^{2}(0,1)$ be such that $y_{0}\geq 0$. Let $y\in C([0,T];L^{2}(0,1))$ be the unique mild solution. Then for each $t\in(0,\infty)$ there exists a positive constant, $c_{t}>0$, such that $y(\cdot,t)\geq c_{t}$.

Proof.

The solution $y$ of the PDE (II), can be represented using the Neumann heat kernel, $K$. That is, there exists a measurable map, $K:(0,\infty)\times[0,1]^{2}\rightarrow[0,\infty)$ such that the mild solution, $y$, is related to $K$ by the following relation:

[TABLE]

for each $t\in(0,\infty)$ and almost every $x\in(0,1)$. From [18][Corollary 2.1], we know that the Neumann heat kernel, $K$, satisfies the following lower bound:

[TABLE]

for each $t>0$ and almost every $x,z\in(0,1)$. From this the lower bound on $y(\cdot,t)$ follows. ∎

Lemma IV.7**.**

Let $y_{0}\in\mathcal{D}(A_{a})$ be such that $y_{0}\geq c$ for some strictly positive constant, $c$. Suppose $f\in W^{2,\infty}(0,1)$ such that $f\geq s$ for some strictly positive constant $s$, $a=1/f$ and $\int_{0}^{1}f(x)dx=\int_{0}^{1}y_{0}(x)dx$. Let $T=\sum_{n=1}^{\infty}\frac{1}{n^{2}}$ be the final time. Define the vector field $v$ by

[TABLE]

with $a=1/f$ in (II) whenever $t\in[\sum_{n=1}^{m-1}\frac{1}{n^{2}},\sum_{n=1}^{m}\frac{1}{n^{2}})$ and $m\in\mathbb{Z}_{+}$. Here, we define $\sum_{n=1}^{m}\frac{1}{n^{2}}=0$ if $m=0$.

Then there exists an $\alpha>0$ such that $v\in L^{\infty}(0,T;L^{\infty}(0,T))$ and the mild solution $y$ of the PDE, (II), satisfies $y(\cdot,T)=f$.

Proof.

Substituting $v(\cdot,t)=\frac{y_{x}}{y}-\alpha m\frac{(ay)_{x}}{y}$ whenever $t\in[\sum_{n=1}^{m-1}\frac{1}{n^{2}},\sum_{n=1}^{m}\frac{1}{n^{2}})$, in (II), it can be seen that the solution of the PDE, (II) exists over each time interval $[0,\sum_{n=1}^{m}\frac{1}{n^{2}})$ for each $m\in\mathbb{Z}_{+}$. This is true because this solution can be constructed from mild solutions of the closed-loop PDE

[TABLE]

and we get the relation $y(\cdot,\sum_{n=1}^{m-1}\frac{1}{n^{2}}+j)=\tilde{y}(0,j)$ for each $j\in[0,\frac{1}{m^{2}})$ and each $m\in\mathbb{Z}_{+}$. Then it follows from lemma IV.4 that

[TABLE]

for each $m\in\mathbb{Z}_{+}$, for some strictly positive constants $M_{0}$ and $\lambda$ independent of $m$. Since the summation $\sum_{n=1}^{m}\frac{1}{n}$ is diverging we have $y(\cdot,T)=f$ if the solution is defined over the interval $[0,T]$. By continuity of $y$ on $[0,T)$, it follows that $y\in C([0,T);L^{2}(0,1))$ and can be extended to a unique mild solution $y\in C([0,T];L^{2}(0,1))$ defined over the time interval $[0,T]$.

It is additionally required to prove that $v\in L^{\infty}(0,T;L^{\infty}(0,1))$. As will be shown further ahead, these results will follow if $\alpha>0$ is chosen to be large enough. More specifically, it will be established that if $\alpha$ is large enough we can get uniform bounds on $1/y(\cdot,t)$, $\alpha m(ay)_{x}(\cdot,t)$ and $y_{x}(\cdot,t)$ as $t$ is varied over the interval, $[0,T)$.

First, we derive bounds on the term, $1/y(\cdot,t)$. Due to the lower bound on the initial condition $y_{0}$, and corollary IV.3, it follows that, there exists a positive constant $d>0$ such that

[TABLE]

for all $t\in[0,T)$. This gives us the uniform upper bound $1/d$ on the term $1/y(\cdot,t)$.

Next, we consider the term, $y_{x}(\cdot,t)$. We note that $y_{0}\in\mathcal{D}(A_{a})$. Hence, we can apply the estimates in lemma IV.4 to get

[TABLE]

for some strictly positive constants, $\tilde{M}$. Here, we have implicitly used the fact that $a$ is twice weakly differentiable and the equivalence between the norm, $\|\cdot\|_{g}$ and the norm, $\|\cdot\|_{H^{2}}$. From this it follows that $\|y_{x}(\cdot,t)\|_{\infty}$ is uniformly bounded on the interval, $[0,T)$, due to the continuous embedding, $H^{1}(0,1)\hookrightarrow L^{\infty}(0,1)$ [8][Theorem 8.8].

Lastly, we need to bound the term, $\alpha m(ay)_{x}(\cdot,t)$. As in the estimates for $y_{x}(\cdot,t)$ in the above arguments, from lemma IV.4, we have the estimate

[TABLE]

for some strictly positive constant, $\tilde{M}_{1}$. The right-hand side in the estimate,(29), is not uniformly bounded for arbitrary $\alpha>0$ due to dependence on $m$. However, we note that $\lim_{m\rightarrow\infty}-ln~{}m+\sum_{i=1}^{m}\frac{1}{n}=\gamma$ where $\gamma>0$ is the Euler-Mascheroni constant [14][section 1.5]. Therefore, by setting $\alpha\geq 1/\lambda$ the right-side becomes uniformly bounded for all $m\in\mathbb{Z}_{+}$.

From the estimates, (27)-(29) it follows that if $\alpha>0$ is large enough, then $v\in L^{\infty}(0,T;L^{\infty}(0,1))$. This concludes the proof. ∎

Remark IV.8**.**

In the above lemma, the choice, $v(\cdot,t)=\frac{y_{x}}{y}-\alpha m\frac{(ay)_{x}}{y}$ is definitely not unique. In fact any control law of the form, $v(\cdot,t)=\frac{y_{x}}{y}-\alpha m^{\beta}\frac{(ay)_{x}}{y}$ for numerous other values of $\beta$ and $\alpha$ will also achieve the desired objective due to the fact that an exponential function of a variable grows faster than a polynomial function, as the variable tends to infinity. Additionally, we could also replace the parameter $m$ with a continuous function, $m(t)$ such that $\int_{0}^{T}m(\tau)d\tau=\infty$.

From the above lemma the following corollary follows.

Corollary IV.9**.**

Let $y_{0}\in\mathcal{D}(A_{a})$ be such that $y_{0}\geq c$ for some strictly positive constant, $c$. Suppose $f\in W^{2,\infty}(0,1)$ such that $f\geq s$, $a=1/f$ and $\int_{0}^{1}f(x)dx=\int_{0}^{1}y_{0}(x)dx$, for some strictly positive constant, $s$. Let $T>0$ be the final time. Then there exists $v\in L^{\infty}(0,T;L^{\infty}(0,T))$ such that the mild solution, $y$, of the PDE, (II), satisfies $y(\cdot,T)=f$.

The above corollary follows from lemma IV.7 using a straightforward scaling argument.

Now, we are ready to state and prove our main theorem, where we relax the assumptions made in the previous corollary on the initial condition $y_{0}$. A few comments are due before we prove this theorem. There are two main problems with extending corollary IV.9 with the same control as defined in eq. (25), for general positive initial conditions $y_{0}\in L^{2}(0,1)$. Both issues can cause $v$ to tend to $\infty$ as $t\rightarrow 0$. Firstly, with the same control law as in (25), $y_{0}$ needs to have a strict lower bound, as assumed in lemma IV.7. Otherwise the term in the denominator, $y$, causes blow up in $v(⋅,t)=\frac{y_{x}}{y}-\alpha m\frac{(ay)_{x}}{y}$ near $t=0$. Secondly, for general $y_{0}\in L^{2}(0,1)$ the numerator terms ($y_{x}$ and $(ay)_{x}$)) in $v(⋅,t)=\frac{y_{x}}{y}-\alpha m\frac{(ay)_{x}}{y}$ also cause blow up near $t=0$ because it might be true that $y_{0}\notin\mathcal{D}(A_{a})$. These issues can be remedied, as shown in the following proof, by modifying the control in eq. (25) appropriately.

Theorem IV.10**.**

Let $y_{0}\in L^{2}(0,1)$ be such that $y_{0}\geq 0$ and $\int_{0}^{1}f(x)dx=1$. Suppose $f\in W^{2,\infty}(0,1)$. Let $T>0$ be the final time. Then there exists $v\in L^{\infty}(0,T;L^{\infty}(0,1))$, such that the unique mild solution,$y$, of the PDE, (II), satisfies $y(T)=f$.

Proof.

Set $v(\cdot,t)=0$, in (II), for each $t\in[0,\epsilon/2]$ where $\epsilon\in(0,T)$ is small enough. Then the PDE, (II), is the heat equation with Neumann boundary condition. From lemma IV.6, it follows that the solution, $y$, satisfies $y(\cdot,\epsilon)\geq c$ for some strictly positive constant $c$. Then for each $t\in(\epsilon/2,\epsilon]$ let $v(\cdot,t)=\frac{f_{x}}{f}$. Semigroups generated by elliptic operators are analytic. Hence, from regularizing properties of analytic semigroups, [19][Theorem 2.1.1], it follows that $y(\cdot,\epsilon)\in\mathcal{D}(B_{f})$ where $B_{f}:\mathcal{D}(B_{f})\rightarrow L^{2}(0,1)$ is the operator given by

[TABLE]

for each $u\in{D}(B_{f})=\{w\in H^{2}(0,1);(w_{x}-\frac{f_{x}}{f}w)(0)=(w_{x}-\frac{f_{x}}{f}w)(1)=0\}$. From theorem IV.2 we know this implies $y(\cdot,\epsilon)\in\mathcal{D}(A_{a})$. Then the result follows from corollary IV.9. ∎

Remark IV.11**.**

(The case when $f\in W^{1,\infty}(0,1)$) Comparing theorem IV.10 and theorem IV.1, it is apparent that there is a gap in the result we have obtained. That is, while theorem IV.1 states that any strictly positive, at least once differentiable function can be reached asymptotically, theorem IV.10 requires the target densities to be at least twice differentiable in order to be reachable in finite time. However, this assumption can be relaxed by modifying the argument in lemma IV.7. Particularly, if $f$ is only once differentiable then it is no longer true that $\mathcal{D}(A_{a})$ is a subset of $H^{2}(0,1)$, in general. However, even if $u\notin H^{2}(0,1)$ we have that $\|(ay)_{xx}\|_{2}<\infty$. From this, and bounds on $\|ay\|_{2}$, it follows that $\|(ay)_{x}\|_{\infty}<\infty$ and hence, using the product rule, it follows that $\|y_{x}\|_{\infty}<\infty$.

However, $B_{f}$ needs to be defined using its weak formulation in this case to make sense of the term $(\frac{f_{x}}{f}u)_{x}$. We avoid this issue for now and leave the more general case when $f\in W^{1,\infty}(0,1)$ for future work.

As pointed out in the last remark, using the approach in this paper, the requirement that $f\in W^{2,\infty}(0,1)$, can be relaxed. On the other hand, it is not clear how much regularity needs to be assumed when extending the technique in this paper to the case when the diffusion process evolves on a higher dimensional Euclidean space since embedding results depend on the dimension of the domain. However, if the constraint that $v$ is bounded is relaxed to admit square-integrable vector fields, then it is sufficient to establish the bounds on $\|\nabla y(\cdot,t)\|_{2}$ and $\|\nabla(ay)(\cdot,t)\|_{2}$, which immediately follows from estimates such as those in 28 and 29. Such estimates on the $L^{2}$ norm of the control are not dimension dependent.

We would also like to point out that in proving controllability properities of the system (II), we have taken advantage of the fact that diffusion enables infinite speed of propagation of the solution $y$ of the PDE (II) (lemma IV.6). Hence, the control laws constructed in this paper might need further modification if implemented in practice, since robots have limitations on their speed of movement. One possibility is to introduce ’virtual particles’ that do propagate at infinite speeds and hence avoiding the division-by-zero in the control law (25). Another possibility is to have a more realistic model of noise in the system.

V CONCLUSION

In this paper, we proved controllability properties of the Fokker-Planck equation with zero-flux boundary condition. In contrast to previous work, we established controllability with bounded control inputs. Our approach to establishing controllability using spectral properties of the elliptic operators under consideration is also novel. In our opinion, this provides a simpler approach to conclude controllability than methods in previous similar works. Future work will focus on extending the arguments in this paper to the case where the diffusion process evolves on higher-dimensional domains.

VI APPENDIX

VI-A Proof of Theorem IV.2 (See page IV.2)

See IV.2

Proof.

We define the operator, $\tilde{A}_{a}:\mathcal{D}(\tilde{A}_{a})\rightarrow L^{2}(0,1)$, defined by $\tilde{A}_{a}u=a(x)u_{xx}$ for each $u\in{D}(\tilde{A}_{a})=\{w\in H^{2}(0,1);w_{x}(0)=w_{x}(1)=0\}$. Using the product rule for functions in Sobolev spaces, we can represent the operator, $\tilde{A}$, as $\tilde{A}u=(au_{x})_{x}-a_{x}u_{x}$ for each $u\in\mathcal{D}(\tilde{A}_{a})$.

The validity of the product rule of differentiation used above, that is $(pq)_{x}=p_{x}q+pq_{x}$ whenever $p,q\in H^{1}(0,1)$, can be seen by constructing an approximating sequence $p^{n}_{x}$ in $C^{\infty}(0,1)$ converging to $p$ in $H^{1}(0,1)$. Then it can be shown the integral $-\int_{0}^{1}p^{n}(s)q(s)\phi_{x}(s)ds=\int_{0}^{1}(p^{n}_{x}(s)q(s)+p^{n}(s)q_{x}(s))\phi(s)ds$ for each $\phi\in C^{\infty}(0,1)$. Then taking the limit $n\rightarrow\infty$ gives us the validity of the product rule.

We define the bilinear form, $b:H^{1}(0,1)\times H^{1}(0,1)\rightarrow\mathbb{R}$, corresponding to this operator by,

[TABLE]

for each $u,\phi\in H^{1}(0,1)$. $A_{a}$ is related to $b$, by

[TABLE]

for all $u\in\mathcal{D}(\tilde{A}_{a})$ and all $\phi\in H^{1}(0,1)$. Using the above representation of $\tilde{A}_{a}$, from [23][Corollary 4.3] it follows that since $b$ is an accretive, closed and continuous bilinear form on $H^{1}(0,1)$, the associated operator, $\tilde{A}_{a}$ generates a positivity preserving semigroup, $(S(t))_{t\geq 0}$, on $L^{2}(0,1)$ such that for each $y_{0}\in L^{2}(0,1)$ in (LABEL:eq:varsys1), the unique mild solution of the PDE can be represented by $y(\cdot,t)=S(t)y_{0}$. Note that the above mentioned properties of the bilinear form have been established in [23]. By positivity preserving we mean that if $y_{0}\geq 0$ then $S(t)y_{0}\geq 0$ for each $t\in[0,\infty)$. Additionally, we note that $\tilde{A}_{a}$ has a eigenvalue at [math]. The function $\mathbf{1}$, defined by $\mathbf{1}(x)=1$ for almost every $x\in(0,1)$, is an eigenvector corresponding to this eigenvalue, [math]. Hence, if $y_{0}\geq c>0$ for some strictly positive parameter, $c$. Then $S(t)y_{0}=S(t)(c\mathbf{1}+y_{0}-c\mathbf{1})=c\mathbf{1}+S(t)(y_{0}-c\mathbf{1})$, for each $t\in[0,\infty)$. Since $(S(t))_{t\geq 0}$ is positivity preserving, it follows that $S(t)(y_{0}-c\mathbf{1})\geq 0$ and hence $c$ is a lower bound on the solution $S(t)y_{0}$ if $c$ is a lower bound on the initial condition, $y_{0}$. ∎

VI-B Proof of Theorem IV.5 (See page IV.5)

See IV.5

Proof.

We only prove that $\mathcal{D}(B_{f})=\mathcal{D}(A_{a})$. The proof for the other statements follows the same line of arguments as in theorem IV.2.

The result, $\mathcal{D}(B_{f})=\mathcal{D}(A_{a})$, follows from the quotient rule applied at the boundary. That is, (au)_{x}(0)=\big{(}\frac{fu_{x}-uf_{x}}{f^{2}}\big{)}(0)=0 implies $(u_{x}-u\frac{f_{x}}{f})(0)=0$. The quotient rule for functions in Sobolev spaces can be seen to be true by an argument similar to the one made in theorem IV.2 for verification of the product rule.

∎