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

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.

## Authors

• 6 publications
• 2 publications
• 8 publications
06/21/2022

### Parameter estimation for a linear parabolic SPDE model in two space dimensions with a small noise

We study parameter estimation for a linear parabolic second-order stocha...
03/08/2021

### Efficient numerical approximation of a non-regular Fokker–Planck equation associated with first-passage time distributions

In neuroscience, the distribution of a decision time is modelled by mean...
06/29/2020

### Using Reinforcement Learning to Herd a Robotic Swarm to a Target Distribution

In this paper, we present a reinforcement learning approach to designing...
11/14/2017

### Weak convergence of Galerkin approximations for fractional elliptic stochastic PDEs with spatial white noise

The numerical approximation of the solution u to a stochastic partial di...
03/01/2022

### Nonlocal diffusion models with consistent local and fractional limits

For some spatially nonlocal diffusion models with a finite range of nonl...
04/19/2022

### On the Dynamics of Inference and Learning

Statistical Inference is the process of determining a probability distri...
09/25/2017

### Bayesian Filtering for ODEs with Bounded Derivatives

Recently there has been increasing interest in probabilistic solvers for...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 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 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 agents deployed on the one-dimensional domain . The position of each agent, indexed by , evolves according to a stochastic process , where

denotes time. Since the random variables that correspond to the dynamics of each agent are independent and identically distributed, we can drop the subscript

and define the problem in terms of a single stochastic process, . The deterministic motion of each agent is defined by a vector field , where . This motion is perturbed by the Wiener process , 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 , the stochastic process satisfies the following SDE [28]:

 dZ(t) = v(Z,t)dt+dW(t)+dψ(t), (1) Z(0) = Z0,

where is called the reflecting function, a stochastic process that constrains to the domain .

###### Problem II.1.

Given and such that , determine if there exists a feedback control law such that the process (1) satisfies for each Borel subset .

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

 yt=yxx−(vy)x in  (0,1)×[0,T] y(⋅,0)=y0 in  (0,1) (yx−vy)(0,⋅)=(yx−vy)(1,⋅)=0 in  [0,T]

The solution of this equation represents the probability density of a single agent occupying position at time , 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 for all and all measurable . Problem II.1 can be reframed in terms of equation (II) as a PDE controllability problem as follows:

###### Problem II.2.

Given , , and such that , determine whether there exists a space- and time-dependent parameter such that the solution of the PDE (II) satisfies .

## Iii Preliminaries and Notation

We define as the Hilbert space of square-integrable real-valued functions over the unit interval, . The Hilbertian structure of is induced by the standard inner product, , given by:

 ⟨p,q⟩2=∫10p(x)q(x)dx (3)

for each . The norm on the space is defined as

 ∥p∥2=⟨p,p⟩1/22 (4)

for all . For a function and a given constant , we write to imply that for almost every .

We define the Sobolev space . We equip this space with the usual Sobolev norm , given by

 ∥p∥H1=(∥p∥22+∥px∥22)1/2

for each . We will also need the space , which will be equipped with the norm given by

 ∥p∥H2=(∥p∥2H1+∥pxx∥22)1/2

for each . We define as the space of essentially bounded measurable functions on . The space is equipped with the norm

 ∥z∥∞=ess supx∈(0,1)|z(x)|, (5)

where denotes the essential supremum attained by its argument over the interval . For a given , will refer to the set of all functions such that

 ∫10|p(x)|2a(x)dx<∞. (6)

The space is a Hilbert space with respect to the weighted inner product , given by:

 ⟨p,q⟩a=∫10p(x)q(x)a(x)dx (7)

for each . We also define the Sobolev spaces

 W1,∞(0,1)={z∈L∞(0,1):zx∈L∞(0,1)}, (8)
 (9)

The space consists of all continuous functions for which

 ∥u∥C([0,T];L2(0,1)) := max0≤t≤T∥u(t)∥2 < ∞.

We will need an appropriate notion of a solution of the PDE (II). Toward this end, let be a closed linear operator on a Hilbert space with domain . For a given time , a mild solution of the ODE

 ˙u(t)=Au(t);u(0)=u0∈L2(0,1) (10)

is a function such that for each . Under appropriate conditions satisfied by , the mild solution is given by a semigroup of linear operators, , that are generated by the operator . That is, the solution of the above ODE is given by for each .

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 be a closed linear operator with domain for each . For a certain time interval , a piecewise constant family of operators is given by a map, , for which there exists a partition of such that for each and for each . Then a mild solution of the ODE

 ˙u(t)=A(t)u(t);u(0)=u0∈X (11)

is a function such that

 u(t)=u0+∑i∈Z+Ai∫min{t,ai+1}min{t,ai}u(s)ds (12)

for each . 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 and patching these solutions together to obtain the solution . Note that the mild solution is defined with respect to an operator or collection of operators ; 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 and such that , for some strictly positive constant and . Suppose for all and almost every in the PDE, (II), then a unique mild solution

of the PDE exists and satisfies the estimate

 ∥y(⋅,t)−f∥2≤Me−λt∥y0−f∥2 (13)

for some positive constants, , and all .

The above result is a straightforward corollary of a result well-known in mathematics and physics literature; namely, that if is the gradient of a potential function , then is the unique stationary distribution of the stochastic process (1). From this, if one observes that

, then it follows that any weakly differentiable function with a strictly positive lower bound can be written as an exponential distribution,

. Hence, any function, , 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, . 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 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 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 with the zero-flux boundary condition, is a self-adjoint, negative semi-definite, closed and densely defined operator on the weighted space . is isomorphic to the space due to the upper and lower bounds on . Hence, the spectrum of the generator,

, lies on the real line and is bounded above by the principal eigenvalue.

is an eigenvalue that has a positive eigenvector,

. This implies that 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, . From this the result on expoenential stability of the steady state solution follows. See also [3] for more details to establish the spectral gap, , using Poincare inequalities.

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

###### Theorem IV.2.

Consider the PDE,

 yt=ayxx in  (0,1)×[0,T] y(⋅,0)=y0 in  (0,1) yx(0,⋅)=yx(1,⋅)=0 in  [0,T]

Let be such that for some strictly positive constant . Let . Then a unique mild solution of the above PDE exists. Additionally, if there exists a positive constant, , such that , then for each .

###### Proof.

See appendix. ∎

Now, we make some definitions which will be used subsequently. Given such that for some strictly positive constant we define the operator, , by

 Aau=(au)xx (15)

for each .

###### Corollary IV.3.

Let be such that for some strictly positive constant . Consider the PDE,

 yt=(ay)xx in  (0,1)×[0,T] y(⋅,0)=y0 in  (0,1) (ay)x(0,⋅)=(ay)x(1,⋅)=0 in  [0,T]

Let . Then a unique mild solution of the above PDE exists. Additionally, if there exists a positive constant, , such that , then for each , for some strictly positive constant .

###### 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, , transforms the PDE (LABEL:eq:cllp1), to the PDE,

 ut=auxx in  (0,1)×[0,T] u(⋅,0)=u0 in  (0,1) ux(0,⋅)=ux(1,⋅)=0 in  [0,T]

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

###### Lemma IV.4.

Let be such that . Additionally, assume , where such that for some constant , and . Then the mild solution, , of the PDE, (LABEL:eq:cllp1), satisfies for each . Moreover, the following estimates hold:

 ∥y(⋅,t)−f∥2 ≤M0e−λt (18) ∥Aay(⋅,t)∥2 ≤M1e−λt (19)

for some strictly positive constants , and .

###### Proof.

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

For the second estimate, we consider the space where is the graph norm given by

 ∥z∥g=∥z∥2+∥Aaz∥2 (20)

for each . generates a semigroup of operators, , on the space, . Let be the semigroup of operators generated on the space, by . is the restriction of the operator, , for each , on the space . for each . It follows that and are similar for each . Hence, the spectrum and growth bounds of and are the same. Hence, it follows that we have

 ∥y(⋅,t)−f∥g≤M1e−λt (21)

for some and for all . This implies the estimate, (19), since . ∎

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 and be such that for some positive constant . Suppose for all and almost every in the PDE, (II), If there exists a positive constant, , such that , then the unique mild solution of the PDE satisfies the estimate for each , and for some positive constant, .

Moreover, we have that , whenever and is the operator defined in equation (15) and is the operator given by

 Bfu=uxx−(fxfu)x (22)

for each .

###### Lemma IV.6.

Consider the heat equation with Neumann boundary condition, that is, in (II). Let be such that . Let be the unique mild solution. Then for each there exists a positive constant, , such that .

###### Proof.

The solution of the PDE (II), can be represented using the Neumann heat kernel, . That is, there exists a measurable map, such that the mild solution, , is related to by the following relation:

 y(x,t)=∫10K(t,x,z)y0(x)dz (23)

for each and almost every . From [18][Corollary 2.1], we know that the Neumann heat kernel, , satisfies the following lower bound:

 K(t,x,z)≥1(4πt)1/2exp(−(x−z)24t) (24)

for each and almost every . From this the lower bound on follows. ∎

###### Lemma IV.7.

Let be such that for some strictly positive constant, . Suppose such that for some strictly positive constant , and . Let be the final time. Define the vector field by

 v(⋅,t)=yxy−αm(ay)xy (25)

with in (II) whenever and . Here, we define if .

Then there exists an such that and the mild solution of the PDE, (II), satisfies .

###### Proof.

Substituting whenever , in (II), it can be seen that the solution of the PDE, (II) exists over each time interval for each . This is true because this solution can be constructed from mild solutions of the closed-loop PDE

 ~yt=αm(a~y)xx in  (0,1)×[0,1m2) ~y(⋅,0)=~y0=y(⋅,m−1∑n=11n2) in  (0,1) (a~y)x(0,⋅)=(a~y)x(1,⋅)=0 in  [0,1m2)

and we get the relation for each and each . Then it follows from lemma IV.4 that

 |y(⋅,m∑n=11n2)−f∥L2(Ω)≤M0e−αλ∑mn=1nn2 =M0e−αλ∑mn=11n)

for each , for some strictly positive constants and independent of . Since the summation is diverging we have if the solution is defined over the interval . By continuity of on , it follows that and can be extended to a unique mild solution defined over the time interval .

It is additionally required to prove that . As will be shown further ahead, these results will follow if is chosen to be large enough. More specifically, it will be established that if is large enough we can get uniform bounds on , and as is varied over the interval, .

First, we derive bounds on the term, . Due to the lower bound on the initial condition , and corollary IV.3, it follows that, there exists a positive constant such that

 y(⋅,t)≥d (27)

for all . This gives us the uniform upper bound on the term .

Next, we consider the term, . We note that . Hence, we can apply the estimates in lemma IV.4 to get

 ∥yx(⋅,m∑i=11n2)∥H1≤~Me−αλ∑mi=11n (28)

for some strictly positive constants, . Here, we have implicitly used the fact that is twice weakly differentiable and the equivalence between the norm, and the norm, . From this it follows that is uniformly bounded on the interval, , due to the continuous embedding, [8][Theorem 8.8].

Lastly, we need to bound the term, . As in the estimates for in the above arguments, from lemma IV.4, we have the estimate

 ∥αm(ay)x(⋅,m∑i=11n2∥H1≤αm~M1e−αλ∑mi=11n. (29)

for some strictly positive constant, . The right-hand side in the estimate,(29), is not uniformly bounded for arbitrary due to dependence on . However, we note that where is the Euler-Mascheroni constant [14][section 1.5]. Therefore, by setting the right-side becomes uniformly bounded for all .

From the estimates, (27)-(29) it follows that if is large enough, then . This concludes the proof. ∎

###### Remark IV.8.

In the above lemma, the choice, is definitely not unique. In fact any control law of the form, for numerous other values of and 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 with a continuous function, such that .

From the above lemma the following corollary follows.

###### Corollary IV.9.

Let be such that for some strictly positive constant, . Suppose such that , and , for some strictly positive constant, . Let be the final time. Then there exists such that the mild solution, , of the PDE, (II), satisfies .

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 . 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 . Both issues can cause to tend to as . Firstly, with the same control law as in (25), needs to have a strict lower bound, as assumed in lemma IV.7. Otherwise the term in the denominator, , causes blow up in near . Secondly, for general the numerator terms ( and )) in also cause blow up near because it might be true that . These issues can be remedied, as shown in the following proof, by modifying the control in eq. (25) appropriately.

###### Theorem IV.10.

Let be such that and . Suppose . Let be the final time. Then there exists , such that the unique mild solution,, of the PDE, (II), satisfies .

###### Proof.

Set , in (II), for each where is small enough. Then the PDE, (II), is the heat equation with Neumann boundary condition. From lemma IV.6, it follows that the solution, , satisfies for some strictly positive constant . Then for each let . Semigroups generated by elliptic operators are analytic. Hence, from regularizing properties of analytic semigroups, [19][Theorem 2.1.1], it follows that where is the operator given by

 Bfu=uxx−(fxfu)x (30)

for each . From theorem IV.2 we know this implies . Then the result follows from corollary IV.9. ∎

###### Remark IV.11.

(The case when ) 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 is only once differentiable then it is no longer true that is a subset of , in general. However, even if we have that . From this, and bounds on , it follows that and hence, using the product rule, it follows that .

However, needs to be defined using its weak formulation in this case to make sense of the term . We avoid this issue for now and leave the more general case when for future work.

As pointed out in the last remark, using the approach in this paper, the requirement that , 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 is bounded is relaxed to admit square-integrable vector fields, then it is sufficient to establish the bounds on and , which immediately follows from estimates such as those in 28 and 29. Such estimates on the 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 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, , defined by for each . Using the product rule for functions in Sobolev spaces, we can represent the operator, , as for each .

The validity of the product rule of differentiation used above, that is whenever , can be seen by constructing an approximating sequence in converging to in . Then it can be shown the integral for each . Then taking the limit gives us the validity of the product rule.

We define the bilinear form, , corresponding to this operator by,

 b(u,ϕ)=⟨aux,ϕx⟩2+⟨axux,ϕ⟩2 (31)

for each . is related to , by

 ⟨~Aau,ϕ⟩2=−b(u,ϕ) (32)

for all and all . Using the above representation of , from [23][Corollary 4.3] it follows that since is an accretive, closed and continuous bilinear form on , the associated operator, generates a positivity preserving semigroup, , on such that for each in (LABEL:eq:varsys1), the unique mild solution of the PDE can be represented by . Note that the above mentioned properties of the bilinear form have been established in [23]. By positivity preserving we mean that if then for each . Additionally, we note that has a eigenvalue at . The function , defined by for almost every , is an eigenvector corresponding to this eigenvalue, . Hence, if for some strictly positive parameter, . Then , for each . Since is positivity preserving, it follows that