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

arXiv:1703.07306v2 [] 24 Mar 2017

Karthik Elamvazhuthi, Hendrik Kuiper, and Spring Berman
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.

I. INTRODUCTION
In recent years, there has been much work on the modeling
and control of swarms of homogeneous agents using meanfield 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
This work was supported by National Science Foundation (NSF) Award
CMMI-1436960 and by ONR Young Investigator Award N00014-16-1-2605.
Karthik Elamvazhuthi and Spring Berman are with the School
for Engineering of Matter, Transport and Energy, Arizona
State University, Tempe, AZ, 85281 USA {karthikevaz,
Spring.Berman}@asu.edu. Hendrik Kuiper is with the School of
Mathematical and Statistical Sciences, Arizona State University, Tempe,
AZ, 85281 USA {kuiper@asu.edu}

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 Rn 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 FokkerPlanck 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 onedimensional domain [0, 1]. The position of each agent, indexed by i ∈ {1, 2, ..., n}, evolves according to a stochastic
process Zi (t) ∈ [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) ∈ [0, 1]. The deterministic
motion of each agent is defined by a vector field v(x, t) ∈ R,
where x ∈ [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]:
dZ(t) =
Z(0) =

v(Z, t)dt + dW (t) + dψ(t),
Z0 ,

(1)

where dψ(t) ∈ 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] → R+ such that
R1
f (x)dx = 1, determine if there exists a feedback control
0
law v : [0, 1] × [0,R T ] → R such that the process (1) satisfies
P(Z(T ) ∈ Γ) = Γ f (x)dx for each Borel subset Γ ⊂ [0, 1].
The Kolmogorov forward equation corresponding to the
SDE (1) is given by [25]:

+
Problem II.2. Given T >R 0, y0 : [0, 1]
R 1 → R , and
1
+
f : [0, 1] → R such that 0 y0 (x)dx = 0 f (x)dx = 1,
determine whether there exists a space- and time-dependent
parameter v : [0, 1] × [0, T ] → R such that the solution y of
the PDE (2) satisfies y(·, T ) = f .

III. PRELIMINARIES AND NOTATION
We define L2 (0, 1) as the Hilbert space of squareintegrable real-valued functions over the unit interval,
(0, 1) ⊂ R. The Hilbertian structure of L2 (0, 1) is induced by
the standard inner product, h·, ·i2 : L2 (0, 1) × L2 (0, 1) → R,
given by:
Z
1

p(x)q(x)dx

hp, qi2 =

(3)

0

for each p, q ∈ L2 (0, 1). The norm k·k2 on the space L2 (0, 1)
is defined as
1/2
kpk2 = hp, pi2
(4)
for all p ∈ L2 (0, 1). For a function r ∈ L2 (0, 1) and a given
constant c, we write r ≥ c to imply that r(x) ≥ c for almost
every x ∈ (0, 1).

We define the Sobolev space H 1 (0, 1) = z ∈ L2 (0, 1) :
zx ∈ L2 (0, 1) . We equip this space with the usual Sobolev
norm k · kH 1 , given by
1/2

kpkH 1 = kpk22 + kpx k22
for each p ∈ H 1 (0, 1). We will also need the space
H 2 (0, 1) = z ∈ H 1 (0, 1) : zxx ∈ L2 (0, 1) , which will
be equipped with the norm k · kH 2 given by
1/2

kpkH 2 = kpk2H 1 + kpxx k22
for each p ∈ H 1 (0, 1). We define L∞ (0, 1) as the space
of essentially bounded measurable functions on (0, 1). The
space L∞ (0, 1) is equipped with the norm
kzk∞ = ess supx∈(0,1)|z(x)|,

(5)

where ess supx∈(0,1)(·) denotes the essential supremum
attained by its argument over the interval (0, 1). For a given
a ∈ L∞ (0, 1), L2a (0, 1) will refer to the set of all functions
p such that
Z
1

|p(x)|2 a(x)dx < ∞.

(6)

0

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 space L2a (0, 1) is a Hilbert space with respect to the
weighted inner product h·, ·ia : L2a (0, 1) × L2a (0, 1) → R,
given by:
Z
1

(2)

hp, qia =

p(x)q(x)a(x)dx

(7)

0

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

for each f, g ∈ L2 (0, 1). We also define the Sobolev spaces

W 1,∞ (0, 1) = z ∈ L∞ (0, 1) : zx ∈ L∞ (0, 1) , (8)

W 2,∞ (0, 1) = z ∈ W 1,∞ (0, 1) : zxx ∈ L∞ (0, 1) . (9)

The space C([0, T ]; L2 (0, 1)) consists of all continuous
functions u : [0, T ] → L2 (0, 1) for which
kukC([0,T ];L2(0,1)) :=

max ku(t)k2 < ∞.

0≤t≤T

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

v(x, t) = fx (x)/f (x) for all t ∈ [0, ∞) and almost every
x ∈ (0, 1) in the PDE, (2), then a unique mild solution
y ∈ C([0, ∞); L2 (0, 1)) of the PDE exists and satisfies the
estimate

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

ky(·, t) − f k2 ≤ M e−λt ky0 − f k2

(10)

2

is a function
u ∈ C([0, T ]; L (0, 1)) such that u(t) =
Rt
u0 + A 0 u(s)ds for each t ∈ [0, T ]. Under appropriate
conditions satisfied by A, the mild solution is given by a
semigroup of linear operators, (T (t))t≥0 , that are generated
by the operator A. That is, the solution of the above ODE
is given by u(t) = T (t)u0 for each t ∈ [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 Ai be a closed linear operator with domain D(Ai ) for
each i ∈ Z+ . For a certain time interval [0, T ], a piecewise
constant family of operators is given by a map, t 7→ A(t),
for which there exists a partition of [0, T ] = ∪i∈Z+ [ai , ai+1 )
such that ai ≤ ai+1 for each i ∈ Z+ and A(t) = Ai for
each t ∈ [ai , ai+1 ). Then a mild solution of the ODE
u̇(t) = A(t)u(t); u(0) = u0 ∈ X
is a function u ∈ C([0, T ]; X) such that
Z min{t,ai+1 }
X
u(t) = u0 +
Ai
u(s)ds
i∈Z+

(11)

(12)

min{t,ai }

for each t ∈ [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
[ai , ai+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 y0 ∈ L2 (0, 1) and f ∈ W 1,∞ (0, 1)
such that y0 ≥ R0, f ≥ k > 0 Rfor some strictly positive
1
1
constant k and 0 y0 (x)dx = 0 f (x)dx = 1. Suppose

(13)

for some positive constants, M ≥ 1, λ and all t ∈ [0, ∞).
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 φ, then eφ
is the unique stationary distribution of the stochastic process
(1). From this, if one observes that fx /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, (2), 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 (·)xx −(· ffx )x with
the zero-flux boundary condition, is a self-adjoint, negative
semi-definite, closed and densely defined operator on the
weighted space L2f (0, 1). L2f (0, 1) is isomorphic to the space
L2 (0, 1) due to the upper and lower bounds on f . Hence, the
spectrum of the generator, (·)xx − ( ffx )x , lies on the real line
and is bounded above by the principal eigenvalue. 0 is an
eigenvalue that has a positive eigenvector, f . This implies
that 0 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 f 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
y(·, 0) = y0
yx (0, ·) = yx (1, ·) = 0

in (0, 1) × [0, T ]
in (0, 1)
in [0, T ]
(14)

Let a ∈ W 2,∞ (0, 1) be such that a ≥ e for some strictly
positive constant e. Let y0 ∈ L2 (0, 1). Then a unique mild
solution y ∈ C([0, T ]; L2 (0, 1)) of the above PDE exists.

Additionally, if there exists a positive constant, c > 0, such
that y0 ≥ c, then y(·, t) ≥ c for each t ∈ [0, ∞).
Proof. See appendix.
Now, we make some definitions which will be used
subsequently. Given a ∈ W 2,∞ (0, 1) such that a ≥ e for
some strictly positive constant e we define the operator,
Aa : D(Aa ) → L2 (0, 1), by
Aa u = (au)xx

(15)
2

for each u ∈ D(Aa ) = {w ∈ H (0, 1); (awx )(0) =
(aw)x (1) = 0}.
Corollary IV.3. Let a ∈ W 2,∞ (0, 1) be such that a ≥ e for
some strictly positive constant e. Consider the PDE,
yt = (ay)xx
y(·, 0) = y0

in (0, 1) × [0, T ]
in (0, 1)

(ay)x (0, ·) = (ay)x (1, ·) = 0

in [0, T ]

estimate, the assumption of the regularity of the initial
condition is not required.
For the second estimate, we consider the space (D(A)a , k·
kg ) where k · kg is the graph norm given by
kzkg = kzk2 + kAa zk2

(20)

for each z ∈ D(Aa ). Aa generates a semigroup of operators,
(T1 (t))t≥0 , on the space, (D(Aa ), k · kg ). Let (T (t))t≥0 be
the semigroup of operators generated on the space, L2 (0, 1)
by Aa . T1 (t) is the restriction of the operator, T (t), for each
t ∈ [0, ∞), on the space (D(Aa ), k · kg ). T1 (t) = (I −
Aa )−1 T (t)(I − Aa ) for each t ∈ [0, ∞). It follows that
T1 (t) and T (t) are similar for each t ∈ [0, ∞). Hence, the
spectrum and growth bounds of (T1 (t))t≥0 and (T (t))t≥0
are the same. Hence, it follows that we have
ky(·, t) − f kg ≤ M1 e−λt

(21)

(16)

for some M1 > 0 and for all t ∈ [0, ∞). This implies the
estimate, (19), since Aa f = 0.

Let y0 ∈ L (Ω). Then a unique mild solution y ∈
C([0, T ]; L2 (0, 1)) of the above PDE exists. Additionally, if
there exists a positive constant, c > 0, such that y0 ≥ c, then
y(·, t) ≥ d for each t ∈ [0, ∞), for some strictly positive
constant d.

Controllability of the PDE, (2), 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.

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

Theorem IV.5. Let y0 ∈ L2 (0, 1) and f ∈ W 2,∞ (0, 1) be
such that f ≥ k > 0 for some positive constant k. Suppose
v(x, t) = fx (x)/f (x) for all t ∈ [0, ∞) and almost every
x ∈ (0, 1) in the PDE, (2), If there exists a positive constant,
c > 0, such that y0 ≥ c, then the unique mild solution of the
PDE satisfies the estimate y(·, t) ≥ d for each t ∈ [0, ∞),
and for some positive constant, d.
Moreover, we have that D(Bf ) = D(Aa ), whenever a =
1/f and Aa is the operator defined in equation (15) and
Bf : D(Bf ) → L2 (0, 1) is the operator given by

2

ut = auxx
u(·, 0) = u0

in (0, 1) × [0, T ]
in (0, 1)

ux (0, ·) = ux (1, ·) = 0

in [0, T ]
(17)

In the following lemma, we will establish some estimates
on the rate of convergence of the solution, y, of the PDE
(16), assuming the initial condition is regular enough.
Lemma IV.4. Let y0 ∈ D(Aa ) be such that y0 ≥ 0.
Additionally, assume a = 1/f , where f ∈ W 2,∞
R 1 (0, 1) such
that f ≥ c > 0 for some constant c, and 0 f (x)dx =
R1
2
0 y0 (x)dx. Then the mild solution, y ∈ C([0, ∞); L (0, 1)),
of the PDE, (16), satisfies y(·, t) ∈ D(Aa ) for each t ∈
[0, ∞). Moreover, the following estimates hold:
ky(·, t) − f k2

≤ M0 e−λt

(18)

kAa y(·, t)k2

≤ M1 e−λt

(19)

for some strictly positive constants M0 , M1 and λ.
Proof. The first estimate is just a restatement of
[13][Theorem IV.4], where it was shown that f is the
eigenvector of Aa corresponding to simple principal
eigenvalue 0, and the rest of spectrum is in the left-half
complex plane, left to some negative number −λ. For this

Bf u = uxx − (

fx
u)x
f

for each u ∈ D(Bf ) = {w ∈ H 2 (0, 1); (wx −
(wx − ffx w)(1) = 0}.

(22)
fx
f w)(0)

=

Lemma IV.6. Consider the heat equation with Neumann
boundary condition, that is, v ≡ 0 in (2). Let y0 ∈ L2 (0, 1)
be such that y0 ≥ 0. Let y ∈ C([0, T ]; L2 (0, 1)) be the
unique mild solution. Then for each t ∈ (0, ∞) there exists
a positive constant, ct > 0, such that y(·, t) ≥ ct .
Proof. The solution y of the PDE (2), can be represented
using the Neumann heat kernel, K. That is, there exists a
measurable map, K : (0, ∞)× [0, 1]2 → [0, ∞) such that the
mild solution, y, is related to K by the following relation:
Z 1
y(x, t) =
K(t, x, z)y0 (x)dz
(23)
0

for each t ∈ (0, ∞) and almost every x ∈ (0, 1). From
[18][Corollary 2.1], we know that the Neumann heat kernel,

K, satisfies the following lower bound:
K(t, x, z) ≥

−(x − z)2
1
exp(
)
1/2
4t
(4πt)

it follows that, there exists a positive constant d > 0 such
that
(24)
y(·, t) ≥ d

for each t > 0 and almost every x, z ∈ (0, 1). From this the
lower bound on y(·, t) follows.
Lemma IV.7. Let y0 ∈ D(Aa ) be such that y0 ≥ c for some
strictly positive constant, c. Suppose f ∈ W 2,∞ (0, 1) such
that Rf ≥ s for some
s, a = 1/f
R 1 strictly positive constant
P∞ 1
1
and 0 f (x)dx = 0 y0 (x)dx. Let T =
n=1 n2 be the
final time. Define the vector field v by
yx
(ay)x
− αm
(25)
y
y
Pm
Pm−1
with a = 1/f in (2) wheneverPt ∈ [ n=1 n12 , n=1 n12 )
1
and m ∈ Z+ . Here, we define m
n=1 n2 = 0 if m = 0.
Then there exists an α > 0 such that v ∈
L∞ (0, T ; L∞(0, T )) and the mild solution y of the PDE,
(2), satisfies y(·, T ) = f .
v(·, t) =

x
whenever t ∈
Proof. Substituting v(·, t) = yyx − αm (ay)
y
Pm−1 1 Pm 1
solution
[ n=1 n2 , n=1 n2 ), in (2), it can be seen that the
Pm
of the PDE, (2) exists over each time interval [0, n=1 n12 )
for each m ∈ Z+ . This is true because this solution can be
constructed from mild solutions of the closed-loop PDE

ỹt = αm(aỹ)xx
Pm−1
ỹ(·, 0) = ỹ0 = y(·, n=1

in (0, 1) × [0,
1
n2 )

(aỹ)x (0, ·) = (aỹ)x (1, ·) = 0

1
)
m2

in (0, 1)
1
in [0, 2 )
m
(26)

Pm−1

and we get the relation y(·, n=1 n12 +j) = ỹ(0, j) for each
j ∈ [0, m12 ) and each m ∈ Z+ . Then it follows from lemma
IV.4 that
m
Pm
X
n
1
) − f kL2 (Ω) ≤ M0 e−αλ n=1 n2
|y(·,
2
n
n=1

= M0 e−αλ

Pm

1
n=1 n )

for each m ∈ Z+ , for some strictly positive constants
Pm M0
and λ independent of m. Since the summation n=1 n1 is
diverging we have y(·, T ) = f if the solution is defined over
the interval [0, T ]. By continuity of y on [0, T ), it follows
that y ∈ C([0, T ); L2 (0, 1)) and can be extended to a unique
mild solution y ∈ C([0, T ]; L2 (0, 1)) defined over the time
interval [0, T ].
It is additionally required to prove that v
∈
L∞ (0, T ; L∞(0, 1)). As will be shown further ahead, these
results will follow if α > 0 is chosen to be large
enough. More specifically, it will be established that if α
is large enough we can get uniform bounds on 1/y(·, t),
αm(ay)x (·, t) and yx (·, t) as t is varied over the interval,
[0, T ).
First, we derive bounds on the term, 1/y(·, t). Due to the
lower bound on the initial condition y0 , and corollary IV.3,

(27)

for all t ∈ [0, T ). This gives us the uniform upper bound
1/d on the term 1/y(·, t).
Next, we consider the term, yx (·, t). We note that y0 ∈
D(Aa ). Hence, we can apply the estimates in lemma IV.4 to
get
m
X
Pm 1
1
(28)
)kH 1 ≤ M̃ e−αλ i=1 n
kyx (·,
2
n
i=1
for some strictly positive constants, M̃ . Here, we have implicitly used the fact that a is twice weakly differentiable and
the equivalence between the norm, k·kg and the norm, k·kH 2 .
From this it follows that kyx (·, t)k∞ is uniformly bounded
on the interval, [0, T ), due to the continuous embedding,
H 1 (0, 1) ֒→ L∞ (0, 1) [8][Theorem 8.8].
Lastly, we need to bound the term, αm(ay)x (·, t). As in
the estimates for yx (·, t) in the above arguments, from lemma
IV.4, we have the estimate
m
X
Pm 1
1
kαm(ay)x (·,
(29)
kH 1 ≤ αmM̃1 e−αλ i=1 n .
2
n
i=1
for some strictly positive constant, M̃1 . The right-hand side
in the estimate,(29), is not uniformly bounded for arbitrary
α > 0 due to dependence
on m. However, we note that
Pm
limm→∞ −ln m + i=1 n1 = γ where γ > 0 is the EulerMascheroni constant [14][section 1.5]. Therefore, by setting
α ≥ 1/λ the right-side becomes uniformly bounded for all
m ∈ Z+ .
From the estimates, (27)-(29) it follows that if α > 0 is
large enough, then v ∈ L∞ (0, T ; L∞(0, 1)). This concludes
the proof.
Remark IV.8. In the above lemma, the choice, v(·, t) =
(ay)x
yx
is definitely not unique. In fact any control
y − αm y
x
for numerous other
law of the form, v(·, t) = yyx −αmβ (ay)
y
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 replaceR the parameter m
T
with a continuous function, m(t) such that 0 m(τ )dτ = ∞.
From the above lemma the following corollary follows.

Corollary IV.9. Let y0 ∈ D(Aa ) be such that y0 ≥ c for
some strictly positive constant, c.RSuppose f ∈ RW 2,∞ (0, 1)
1
1
such that f ≥ s, a = 1/f and 0 f (x)dx = 0 y0 (x)dx,
for some strictly positive constant, s. Let T > 0 be the final
time. Then there exists v ∈ L∞ (0, T ; L∞(0, T )) such that
the mild solution, y, of the PDE, (2), satisfies y(·, 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 y0 . 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
y0 ∈ L2 (0, 1). Both issues can cause v to tend to ∞ as
t → 0. Firstly, with the same control law as in (25), y0
needs to have a strict lower bound, as assumed in lemma
IV.7. Otherwise the term in the denominator, y, causes blow
x
near t = 0. Secondly, for
up in v(, t) = yyx − αm (ay)
y
2
general y0 ∈ L (0, 1) the numerator terms (yx and (ay)x ))
x
in v(, t) = yyx − αm (ay)
also cause blow up near t = 0
y
because it might be true that y0 ∈
/ D(Aa ). These issues can
be remedied, as shown in the following proof, by modifying
the control in eq. (25) appropriately.
IV.10. Let y0 ∈ L2 (0, 1) be such that y0 ≥ 0 and
RTheorem
1
f (x)dx = 1. Suppose f ∈ W 2,∞ (0, 1). Let T > 0 be the
0
final time. Then there exists v ∈ L∞ (0, T ; L∞ (0, 1)), such
that the unique mild solution,y, of the PDE, (2), satisfies
y(T ) = f .
Proof. Set v(·, t) = 0, in (2), for each t ∈ [0, ǫ/2] where
ǫ ∈ (0, T ) is small enough. Then the PDE, (2), is the heat
equation with Neumann boundary condition. From lemma
IV.6, it follows that the solution, y, satisfies y(·, ǫ) ≥ c for
some strictly positive constant c. Then for each t ∈ (ǫ/2, ǫ]
let v(·, t) = ffx . Semigroups generated by elliptic operators
are analytic. Hence, from regularizing properties of analytic
semigroups, [19][Theorem 2.1.1], it follows that y(·, ǫ) ∈
D(Bf ) where Bf : D(Bf ) → L2 (0, 1) is the operator given
by
fx
Bf u = uxx − ( u)x
(30)
f
for each u ∈ D(Bf ) = {w ∈ H 2 (0, 1); (wx − ffx w)(0) =
(wx − ffx w)(1) = 0}. From theorem IV.2 we know this
implies y(·, ǫ) ∈ D(Aa ). Then the result follows from
corollary IV.9.
Remark IV.11. (The case when f ∈ W 1,∞ (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 D(Aa ) is a subset of H 2 (0, 1),
in general. However, even if u ∈
/ H 2 (0, 1) we have that
k(ay)xx k2 < ∞. From this, and bounds on kayk2, it follows
that k(ay)x k∞ < ∞ and hence, using the product rule, it
follows that kyx k∞ < ∞.
However, Bf needs to be defined using its weak formulation in this case to make sense of the term ( ffx u)x . We avoid
this issue for now and leave the more general case when
f ∈ W 1,∞ (0, 1) for future work.

As pointed out in the last remark, using the approach
in this paper, the requirement that f ∈ W 2,∞ (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 squareintegrable vector fields, then it is sufficient to establish the
bounds on k∇y(·, t)k2 and k∇(ay)(·, t)k2 , which immediately follows from estimates such as those in 28 and 29. Such
estimates on the L2 norm of the control are not dimension
dependent.
We would also like to point out that in proving controllability properities of the system (2), we have taken advantage
of the fact that diffusion enables infinite speed of propagation
of the solution y of the PDE (2) (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
A. Proof of Theorem IV.2 (See page 3)
Theorem IV.2. Consider the PDE,
yt = ayxx
y(·, 0) = y0
yx (0, ·) = yx (1, ·) = 0

in (0, 1) × [0, T ]
in (0, 1)
in [0, T ]
(14)

Let a ∈ W 2,∞ (0, 1) be such that a ≥ e for some strictly
positive constant e. Let y0 ∈ L2 (0, 1). Then a unique mild
solution y ∈ C([0, T ]; L2 (0, 1)) of the above PDE exists.
Additionally, if there exists a positive constant, c > 0, such
that y0 ≥ c, then y(·, t) ≥ c for each t ∈ [0, ∞).
Proof. We define the operator, Ãa : D(Ãa ) → L2 (0, 1),
defined by Ãa u = a(x)uxx for each u ∈ D(Ãa ) = {w ∈
H 2 (0, 1); wx (0) = wx (1) = 0}. Using the product rule for
functions in Sobolev spaces, we can represent the operator,
Ã, as Ãu = (aux )x − ax ux for each u ∈ D(Ãa ).

The validity of the product rule of differentiation used
above, that is (pq)x = px q + pqx whenever p, q ∈ H 1 (0, 1),
can be seen by constructing an approximating sequence pnx in
1
C ∞ (0, 1) converging
be shown
R 1 nto p in H (0, 1). ThenRit1 can
the integral − 0 p (s)q(s)φx (s)ds = 0 (pnx (s)q(s) +
pn (s)qx (s))φ(s)ds for each φ ∈ C ∞ (0, 1). Then taking the
limit n → ∞ gives us the validity of the product rule.
We define the bilinear form, b : H 1 (0, 1)× H 1 (0, 1) → R,
corresponding to this operator by,
b(u, φ) = aux , φx

2

+ ax u x , φ

2

(31)

for each u, φ ∈ H 1 (0, 1). Aa is related to b, by
hÃa u, φi2 = −b(u, φ)

(32)

for all u ∈ D(Ãa ) and all φ ∈ H 1 (0, 1). Using the above
representation of 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, Ãa generates a
positivity preserving semigroup, (S(t))t≥0 , on L2 (0, 1) such
that for each y0 ∈ L2 (0, 1) in (14), the unique mild solution
of the PDE can be represented by y(·, t) = S(t)y0 . Note
that the above mentioned properties of the bilinear form
have been established in [23]. By positivity preserving we
mean that if y0 ≥ 0 then S(t)y0 ≥ 0 for each t ∈ [0, ∞).
Additionally, we note that Ãa has a eigenvalue at 0. The
function 1, defined by 1(x) = 1 for almost every x ∈ (0, 1),
is an eigenvector corresponding to this eigenvalue, 0. Hence,
if y0 ≥ c > 0 for some strictly positive parameter, c. Then
S(t)y0 = S(t)(c1 + y0 − c1) = c1 + S(t)(y0 − c1), for
each t ∈ [0, ∞). Since (S(t))t≥0 is positivity preserving, it
follows that S(t)(y0 − c1) ≥ 0 and hence c is a lower bound
on the solution S(t)y0 if c is a lower bound on the initial
condition, y0 .
B. Proof of Theorem IV.5 (See page 4)
Theorem IV.5. Let y0 ∈ L2 (0, 1) and f ∈ W 2,∞ (0, 1) be
such that f ≥ k > 0 for some positive constant k. Suppose
v(x, t) = fx (x)/f (x) for all t ∈ [0, ∞) and almost every
x ∈ (0, 1) in the PDE, (2), If there exists a positive constant,
c > 0, such that y0 ≥ c, then the unique mild solution of the
PDE satisfies the estimate y(·, t) ≥ d for each t ∈ [0, ∞),
and for some positive constant, d.
Moreover, we have that D(Bf ) = D(Aa ), whenever a =
1/f and Aa is the operator defined in equation (15) and
Bf : D(Bf ) → L2 (0, 1) is the operator given by
Bf u = uxx − (

fx
u)x
f

for each u ∈ D(Bf ) = {w ∈ H 2 (0, 1); (wx −
(wx − ffx w)(1) = 0}.

(22)
fx
f w)(0)

=

Proof. We only prove that D(Bf ) = D(Aa ). The proof for
the other statements follows the same line of arguments as
in theorem IV.2.
The result, D(Bf ) = D(Aa ), follows from the quotient rule applied at the boundary. That is, (au)x (0) =
f ux −ufx
(0) = 0 implies (ux − u ffx )(0) = 0. The quotient
f2

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.

R EFERENCES
[1] Luigi Ambrosio, Giuseppe Savaré, and Lorenzo Zambotti. Existence
and stability for fokker–planck equations with log-concave reference
measure. Probability theory and related fields, 145(3-4):517–564,
2009.
[2] Mario Annunziato and Alfio Borzı̀. Optimal control of probability
density functions of stochastic processes. Mathematical Modelling
and Analysis, 15(4):393–407, 2010.
[3] Dominique Bakry, Ivan Gentil, and Michel Ledoux. Analysis and
geometry of Markov diffusion operators, volume 348. Springer Science
& Business Media, 2013.
[4] 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.
[5] Alain Bensoussan, Jens Frehse, Phillip Yam, et al. Mean field games
and mean field type control theory, volume 101. Springer, 2013.
[6] A Blaquiere.
Controllability of a fokker-planck equation, the
schrödinger system, and a related stochastic optimal control (revised
version). Dynamics and Control, 2(3):235–253, 1992.
[7] Mattia Bongini, Massimo Fornasier, Francesco Rossi, and Francesco
Solombrino. Mean-field pontryagin maximum principle. arXiv
preprint arXiv:1504.02236, 2015.
[8] Haim Brezis. Functional Analysis, Sobolev Spaces and Partial
Differential Equations. Springer Science & Business Media, 2010.
[9] Roger Brockett. Notes on the control of the liouville equation. In
Control of Partial Differential Equations, pages 101–129. Springer,
2012.
[10] Yongxin Chen, Tryphon T Georgiou, and Michele Pavon. Optimal
steering of a linear stochastic system to a final probability distribution,
part i. IEEE Transactions on Automatic Control, 61(5):1158–1169,
2016.
[11] Yongxin Chen, Tryphon T Georgiou, and Michele Pavon. Optimal
steering of a linear stochastic system to a final probability distribution,
part ii. IEEE Transactions on Automatic Control, 61(5):1170–1180,
2016.
[12] Paolo Dai Pra. A stochastic control approach to reciprocal diffusion
processes. Applied mathematics and Optimization, 23(1):313–329,
1991.
[13] Karthik Elamvazhuthi, Chase Adams, and Spring Berman. Coverage
and field estimation on bounded domains by diffusive swarms. In
Decision and Control (CDC), 2016 IEEE 55th Conference on, pages
2867–2874. IEEE, 2016.
[14] Steven R Finch. Mathematical constants. Cambridge university press,
2003.
[15] Greg Foderaro, Silvia Ferrari, and Thomas A Wettergren. Distributed
optimal control for multi-agent trajectory optimization. Automatica,
50(1):149–154, 2014.
[16] Massimo Fornasier, Benedetto Piccoli, and Francesco Rossi.
Mean-field sparse optimal control.
Phil. Trans. R. Soc. A,
372(2028):20130400, 2014.
[17] Diogo A Gomes et al. Mean field games models: a brief survey.
Dynamic Games and Applications, 4(2):110–154, 2014.
[18] Peter Li and Shing Tung Yau. On the parabolic kernel of the
schrödinger operator. Acta Mathematica, 156(1):153–201, 1986.
[19] Alessandra Lunardi. Analytic semigroups and optimal regularity in
parabolic problems. Springer Science & Business Media, 2012.
[20] Peter A Markowich and Cédric Villani. On the trend to equilibrium
for the fokker-planck equation: an interplay between physics and
functional analysis. Mat. Contemp, 19:1–29, 2000.
[21] Alexandre R Mesquita and João P Hespanha. Jump control of
probability densities with applications to autonomous vehicle motion.
IEEE Transactions on Automatic Control, 57(10):2588–2598, 2012.
[22] Dejan Milutinovic and Pedro Lima. Modeling and optimal centralized
control of a large-size robotic population. IEEE Transactions on
Robotics, 22(6):1280–1285, 2006.
[23] El-Maati Ouhabaz. Analysis of heat equations on domains.(LMS-31).
Princeton University Press, 2009.

[24] Benedetto Piccoli, Francesco Rossi, and Emmanuel Trélat. Control
to flocking of the kinetic Cucker–Smale model. SIAM Journal on
Mathematical Analysis, 47(6):4685–4719, 2015.
[25] Andrey Pilipenko. An Introduction to Stochastic Differential Equations
with Reflection, volume 1. Universitätsverlag Potsdam, 2014.
[26] Alessio Porretta. On the planning problem for the mean field games
system. Dynamic Games and Applications, 4(2):231–256, 2014.
[27] Amanda Prorok, Nikolaus Corell, and Alcherio Martinoli. Multilevel spatial modeling for stochastic distributed robotic systems. The
International Journal of Robotics Research, page 0278364911399521,
2011.
[28] Hiroshi Tanaka et al. Stochastic differential equations with reflecting
boundary condition in convex regions. Hiroshima Mathematical
Journal, 9(1):163–177, 1979.
[29] Cédric Villani. Optimal transport: old and new, volume 338. Springer
Science & Business Media, 2008.

