Complete analytic solution to Brownian unicycle dynamics
Agostino Martinelli
October 8, 2018
Abstract
This paper derives a complete analytical solution for the probability distribution of the
configuration of a non-holonomic vehicle that moves in two spatial dimensions by satisfying
the unicycle kinematic constraints and in presence of Brownian noises. In contrast to pre-
vious solutions, the one here derived holds even in the case of arbitrary linear and angular
speed. This solution is obtained by deriving the analytical expression of any-order moment
of the probability distribution. To the best of our knowledge, an analytical expression for
any-order moment that holds even in the case of arbitrary linear and angular speed, has
never been derived before. To compute these moments, a direct integration of the Langevin
equation is carried out and each moment is expressed as a multiple integral of the determin-
istic motion (i.e., the known motion that would result in absence of noise). For the special
case when the ratio between the linear and angular speed is constant, the multiple integrals
can be easily solved and expressed as the real or the imaginary part of suitable analytic
functions. As an application of the derived analytical results, the paper investigates the
diffusivity of the considered Brownian motion for constant and for arbitrary time-dependent
linear and angular speed.
Keywords: Brownian motion, Stochastic dynamics, Fokker-Plank equation without detailed
balance, Stochastic processes, Localization.
1
arXiv:1501.03300v1  [cs.RO]  14 Jan 2015
Contents
1
Introduction
3
2
Stochastic motion model
4
3
Computation of the probability distribution
5
3.1
First and second-order statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
3.2
Computation of any-order moment . . . . . . . . . . . . . . . . . . . . . . . . . .
7
4
Diffusivity for constant and for arbitrary time-dependent linear and angular
speed
9
4.1
Constant ratio between linear and angular speed . . . . . . . . . . . . . . . . . .
9
4.2
Arbitrary time-dependent ratio between linear and angular speed . . . . . . . . .
13
5
Conclusion
14
A Computation of ⟨u(s)p w(s)q⟩
15
B Computation of
D
u(s)p w(s)q ˜θrE
18
C Computation of

D(s)4
=

u(s)2 w(s)2
20
2
1
Introduction
In recent years a great interest has been devoted to investigate the statistical properties of the
motion of active (or self-propelled) particles. These particles differ from passive particles since
they move under the action of an internal force. Typical examples of these particles are Brownian
motors [9, 10, 23], biological and artificial microswimmers [8, 12, 13, 25, 26, 27, 28], macroscopic
animals [1, 4, 17] and even pedestrians [29].
In the last decade, this investigation has also been addressed in the framework of mobile
robotics, in particular by analyzing the 2D cases of the unicyle, cart, and car [30, 31], and the
3D case of flexible needle steering [20, 22]. In these works, analytical solutions are obtained by
directly solving the corresponding Fokker-Planck equations by using the Fourier transform. An-
other method proposed for obtaining a closed-form solution is the use of exponential coordinates
[14, 21]. In all these works, the benefit of having closed-form solutions is clearly illustrated by
introducing motion planning methods based on them.
In this paper we want to deal with wheeled robots, which represent a special case of active
particles moving in a two dimensional space. Their configuration is characterized by a 3D-vector
(2 components characterize their position and one their orientation). On the other hand, in most
of cases, a wheeled robot cannot move along any direction since it must satisfy the so-called
non-holonomic constraints [11]. Our stochastic motion model (presented in section 2) differs
from the one considered in [14, 21, 30, 31] since we assume independent errors acting on the
linear and angular components. The model adopted in [14, 21, 30, 31] refers to the case of a
differential drive system where the errors acting on each wheel are independent.
Starting from our stochastic motion model, we compute the analytical expression of any-
order moment of the probability distribution (section 3). Specifically, each moment is expressed
as a multiple integral of the deterministic motion performed by the mobile robot (i.e., the known
motion that would result in absence of noise). In other words, the expression is a multiple integral
on two time-dependent functions, which describe the time behaviour of the deterministic linear
and angular speed. This allows obtaining any order moment for any deterministic motion, i.e.,
when the mobile robot is propelled by arbitrary time-dependent linear and angular speed. For
the special case when the ratio between these two speeds is constant, the multiple integrals can be
expressed as the real or the imaginary part of suitable analytic functions. In section 4 we show the
power of the derived analytical results by investigating the diffusivity of the considered Brownian
motion for constant and for arbitrary time-dependent linear and angular speed. Conclusions are
provided in section 5.
The interest in deriving the statistical properties of the motion of a wheeled robot comes
from the possibility of improving its localization (both in precision and speed). The localization
is a fundamental problem in mobile robotics which must be solved to autonomously and safely
navigate. It is a non linear estimation problem. In most of cases, the localization is carried out
by using nonlinear filters (e.g., Extended and Unscented Kalman Filters, Particle filters, etc). On
the other hand, non linear filters are strictly connected with the solution of partial differential
equations (see [3, 5, 6, 7] and the next section 2). Very easily speaking, for a dynamic stochastic
process described by stochastic differential equations, it is possible to introduce a probability
distribution which provides the probability that the process takes a given set of values. This
probability distribution must satisfy a partial differential equation. In this sense, a non linear
filter could be considered as a numerical solution of this partial differential equation. In the
specific case of localization, the dynamic process is the motion of the robot together with the
noisy data delivered by its sensors.
We already computed the statistics up to the second order for a similar motion model [15,
16] (in the expressions in [15] there are some typos, corrected in [16]). Here, we extend the
3
computation in [15] by including any-order statistics. We also show that the expressions hold for
a much more general model (the one described in section 2) and not only for a simple specific
odometry error model (as shown in [15]). As a result, this paper provides a complete analytical
solution to a Fokker-Planck equation for which the detailed balance condition is not verified1.
Note that the goal of this paper is to analytically compute the moments of the probability
distribution up to any order for the configuration of a non-holonomic vehicle. As it will be seen,
this is a very hard task from a computational point of view. These analytical results could be
very useful to improve the localization which is, as previously mentioned, a numerical solution
of a partial differential equation.
2
Stochastic motion model
We consider a mobile robot which moves in a 2D−environment. The configuration of the mobile
robot is characterized by its position and orientation. In Cartesian coordinates, we denote them
by r ≡[x, y]T and θ, respectively. We assume that the motion of the mobile robot satisfies the
unicycle model [11], which is the most general nonholonomic model for 2D−motions, where the
shift of the mobile robot can only occur along one direction ([cos θ, sin θ]T ):


˙x =
v cos θ
˙y =
v sin θ
˙θ =
ω
(1)
The quantities v = v(t) and ω = ω(t) are the linear and angular speed, respectively. We assume
that these functions are known and we call the motion that results from them the deterministic
motion. Now we want to introduce a stochastic model that generalizes equation (1) by accounting
a Gaussian white noise. We assume that also the noise satisfies the same nonholonomic constraint.
The stochastic differential equations are:


dx(t) =
cos θ(t) [v(t)dt + α(t)dwr(t)]
dy(t) =
sin θ(t) [v(t)dt + α(t)dwr(t)]
dθ(t) =
ω(t)dt + β(t)dwθ(t)
(2)
where [dwr(t), dwθ(t)]T is a standard Wiener process of dimension two [18].
The functions
α(t) and β(t) are modeled according to the following physical requirement (diffusion in the
overdamped regime). We require that, when the mobile robot moves during the infinitesimal
time interval dt, the shift is a random Gaussian variable, whose variance increases linearly with
the traveled distance. This is obtained by setting α(t) =
p
Kr|v(t)|, with Kr a positive parameter
which characterizes our system (interaction robot-environment). Similarly, we require that the
rotation accomplished during the same time interval dt is a random Gaussian variable, whose
variance increases linearly with the traveled distance. Hence, we set β(t) =
p
Kθ|v(t)|, with Kθ
another positive parameter which characterizes our system (interaction robot-environment). The
equations in (2) are the Langevin equations in the overdamped limit. From now on, we assume
that the mobile robot can only move ahead2, i.e., v(t) ≥0. We remind the reader that v(t) is a
deterministic and known function of time. The curve length ds = v(t)dt is the curve length of the
deterministic motion, which is independent of the Wiener process [dwr(t), dwθ(t)]T . Hence, we
1Note that, deriving any-order moment, corresponds to analytically derive the probability density distribution
[19].
2Note that this constraint does not limit the robot motion: for instance, a go and back motion is easily obtained
by accomplishing 180deg rotations
4
are allowed to use the deterministic curve length ds = v(t)dt instead of the time t. By denoting
the ratio between the angular and the linear speed by µ, equations in (2) read:


dx(s) =
cos θ(s) [ds + K
1
2r dwr(s)]
dy(s) =
sin θ(s) [ds + K
1
2r dwr(s)]
dθ(s) =
µ(s)ds + K
1
2
θ dwθ(s)
(3)
The associated Smoluchowski equation is:
∂p
∂s = −∇· (D1p) + ∇· (D2∇p) = 0
(4)
where:
• p = p(x, y, θ; s(t)) is the probability density for the mobile robot at the configuration
(x, y, θ) and at time t;
• D1 = [cos θ, sin θ, µ]T is the drift vector;
• D2 = 1
2


Kr cos2 θ
Kr cos θ sin θ
0
Kr cos θ sin θ
Kr sin2 θ
0
0
0
Kθ

is the diffusion tensor.
Our goal is to obtain the probability density p(x, y, θ; s). This will be done in the next two
sections.
3
Computation of the probability distribution
The probability density p(x, y, θ; s) satisfies the Smoluchowski partial differential equation in
(4), which is a special case of the Fokker-Planck equation. Since the detailed balance condition
is not satisfied [24], we follow a different procedure to compute p(x, y, θ; s). Specifically, we
use the Langevin equation in (3) to compute the moments, up to any order, of the probability
distribution. First of all, we remark that the third equation in (3) is independent of the first two,
is independent of θ and is linear in dwθ. As a result, the probability density only in terms of
θ, i.e., the probability Pθ(θ; s) ≡
R
dx
R
dy p(x, y, θ; s), is a Gaussian distribution with mean
value θ0 +
R s
0 ds′µ(s′) and variance Kθs, where θ0 is the initial orientation. This same result
could also be obtained by integrating (4) in x and y and by using the divergence theorem. In
other words, we have the following distribution for the orientation at a given value of s:
θ(s) = N(θ(s), Kθs)
(5)
where N(·, ·) denotes the normal distribution with mean value and variance the first and the
second argument; θ(s) ≡θ0 +
R s
0 ds′ µ(s′).
Note that, according to the unicycle model, a
trajectory is completely characterized by its starting point and by the orientation vs the curve
length, i.e., by the function θ(s). In the following, we will call the deterministic function θ(s),
the deterministic trajectory.
5
3.1
First and second-order statistics
Let us consider the first equation in (3). We compute the expression of x(s) by a direct integra-
tion. We divide the interval (0, s) in N equal segments, δs ≡
s
N . We have:
x(s) = lim
N→∞
N
X
j=1
(δs + ϵj) cos θj
(6)
where ϵj is a random Gaussian variable satisfying ⟨ϵj⟩= 0, ⟨ϵj ϵk⟩= δjkKrδs for j, k = 1, · · · , N
(δjk is the Kronecker delta) and θj ≡θ(jδs). On the other hand, we have:
θj = θ(jδs) +
j
X
m=1
δθm ≡θj + ∆θj
(7)
where, according to our stochastic model in (3), δθj is a random Gaussian variable satisfy-
ing ⟨δθj⟩= 0, ⟨δθj ϵk⟩= 0 and ⟨δθj δθk⟩= δjkKθδs for j, k = 1, · · · , N.
As a result,
∆θj ≡Pj
m=1 δθm is also a random Gaussian variable satisfying ⟨∆θj⟩= 0, ⟨∆θj ∆θk⟩=
 Kθjδs
if j ≤k
Kθkδs
if j > k .
Starting from (6) and (7) and by remarking that ⟨cos ∆θj⟩= e−
kθjδs
2
and ⟨sin ∆θj⟩= 0, it is easy to obtain the mean value of x(s). We have:
⟨x(s)⟩= lim
N→∞
N
X
j=1
(δs + ⟨ϵj⟩) ⟨cos θj⟩=
(8)
= lim
N→∞
N
X
j=1
δs cos θj e−
kθjδs
2
=
Z s
0
ds′ cos θ(s′) e−
Kθs′
2
Similarly, it is possible to obtain the mean value of y(s):
⟨y(s)⟩=
Z s
0
ds′ sin θ(s′) e−
Kθs′
2
(9)
Finally, in a similar manner, but with some more computation, we obtain all the second-order
moments:

x(s)2
=
Z s
0
ds′
Z s−s′
0
ds′′e−
Kθs′′
2
{(1 + χc(s′))
(10)
cos[θ(s′ + s′′) −θ(s′)] −χs(s′) sin[θ(s′ + s′′) −θ(s′)]
	
+
+Kr
2

s +
Z s
0
ds′ χc(s′)


y(s)2
=
Z s
0
ds′
Z s−s′
0
ds′′e−
Kθs′′
2
{(1 −χc(s′))
(11)
cos[θ(s′ + s′′) −θ(s′)] + χs(s′) sin[θ(s′ + s′′) −θ(s′)]
	
+
+Kr
2

s −
Z s
0
ds′ χc(s′)

⟨x(s) y(s)⟩=
Z s
0
ds′
Z s−s′
0
ds′′e−
Kθs′′
2
{χs(s′)
(12)
6
cos[θ(s′ + s′′) −θ(s′)] + χc(s′) sin[θ(s′ + s′′) −θ(s′)]
	
+
+Kr
2
Z s
0
ds′ χs(s′)
σxθ(s) ≡⟨x(s) θ(s)⟩−⟨x(s)⟩θ(s) = 2Kθ
∂⟨y(s)⟩
∂Kθ
(13)
σyθ(s) ≡⟨y(s) θ(s)⟩−⟨y(s)⟩θ(s) = −2Kθ
∂⟨x(s)⟩
∂Kθ
(14)
where χc(s′) ≡cos[2θ(s′)] e−2Kθs′ and χs(s′) ≡sin[2θ(s′)] e−2Kθs′.
Obtaining the expression of higher-order moments will demand more tricky computation and
will be dealt separately, in the next subsection. Here we conclude by considering the quantity
D(s)2 ≡x(s)2 + y(s)2, which provides the time-evolution of the square of the distance of the
mobile robot from its initial position. From (10) and (11) we obtain a simple expression for its
mean value:

D(s)2
=
(15)
Krs + 2
Z s
0
ds′
Z s−s′
0
ds′′e−
Kθs′′
2
cos[θ(s′ + s′′) −θ(s′)]
3.2
Computation of any-order moment
We introduce the following two complex random quantities:
u(s) ≡lim
N→∞
N
X
j=1
(δs + ϵj)eiθj; w(s) ≡lim
N→∞
N
X
j=1
(δs + ϵj)e−iθj
(16)
From (6) and (16) it is immediate to realize that x(s) = u(s)+w(s)
2
. Similarly we have y(s) =
u(s)−w(s)
2i
. Hence, in order to compute any-order moment which involves x(s) and y(s) it suffices
to compute ⟨u(s)p w(s)q⟩for any p, q ∈N. The computation of this quantity requires several
tricky steps, which are provided in appendix A. The key is to separate all the independent random
quantities in order to compute their mean values. This is obtained by arranging all the sums in
a suitable manner (see appendix A for all the details). We have:
⟨u(s)p w(s)q⟩=
⌊p+q
2
⌋
X
n=0
Kn
r
2n
X
l=0

p
2n −l
q
l

(17)
X
m
2n −l
m
 l
m

m!(2n −l −m −1)!!(l −m −1)!! sm
ei(p−q)θ0
2n −l −m
2

!
l −m
2

!(p −2n + l)!(q −l)!
X
c
Z s
0
ds1
Z s
s1
ds2 · · ·
Z s
sβ−1
dsβexp
(β−1
X
b=0
[i(p −q + Φb)·
·[θ(sb+1) −θ(sb)] −(p −q + Φb)2(sb+1 −sb)Kθ
2

where:
7
• ⌊p+q
2 ⌋is the largest integer not greater than p+q
2 ;
• the second sum on l (i.e., P2n
l=0) is restricted to the values of l for which 2n −l ≤p and
l ≤q (or, equivalently, we are using the convention that
x
y


= 0 when y > x);
• the sum on m (i.e., P
m) goes from 0 to the minimum between l and 2n −l and it is
restricted to the integers m with the same parity of l (hence, both 2n−l−m
2
and l−m
2
are
integers);
• the symbols "!" and "!!" denote the factorial and the double factorial, respectively [2] (note
that 0! = 0!! = (−1)!! = 1);
• β = p + q −n −m is the dimension of the remaining multiple integral (note that a multiple
integral of dimension m has already been computed and provided the term sm);
• the sum on c, i.e., P
c, is the sum over all the vectors c of dimension β, whose entries
are −2, −1, 1 and 2: specifically, each vector c has l−m
2
entries equal to 2, q −l equal
to 1, p −2n + l equal to −1 and 2n−l−m
2
equal to −2 (note that the sum P
c consists of
 β
q−l

 β−q+l
p−2n+l

n−m
l−m
2


addends);
• Φb ≡Pb
a=1 ca (note that Φβ = q −p);
• s0 ≡0 and θ(s0) = θ0 (note that s0 is not a variable of integration).
In order to complete the derivation of the statistics for our problem, we need to compute any-
order moment which also involves the orientation θ. We provide the formula for the quantity
D
u(s)p w(s)q ˜θ(s)rE
, where ˜θ(s) ≡θ(s) −θ(s). The details of this computation are provided in
appendix B. We have, for any p, q, r ∈N:
D
u(s)p w(s)q ˜θ(s)rE
=
⌊p+q
2
⌋
X
n=0
Kn
r
2n
X
l=0

p
2n −l
q
l

(18)
X
m
2n −l
m
 l
m

m!(2n −l −m −1)!!(l −m −1)!! sm ei(p−q)θ0
2n −l −m
2

!
l −m
2

!(p −2n + l)!(q −l)!
X
c
X
γ
r!(γβ+1 −1)!!K
r
2
θ
Qβ
b=0 γb+1!
Z s
0
ds1
Z s
s1
ds2 · · ·
Z s
sβ−1
dsβ(s −sβ)
γβ+1
2
exp
(
i
β−1
X
b=0
(p −q + Φb)[θ(sb+1) −θ(sb)]
)
·
·
β−1
Y
b=0

(sb+1 −sb)
γb+1
2 exp

−(p −q + Φb)2Kθ(sb+1 −sb)
2

⌊
γb+1
2
⌋
X
a=0
γb+1!

i(p −q + Φb)
p
Kθ(sb+1 −sb)
γb+1−2a
a!(γb+1 −2a)!2a





8
where the sum over γ, i.e. P
γ, is the sum over all the vectors γ = [γ1, γ2, · · · , γβ, γβ+1], where γb
(b = 1, · · · , β) are positive integers and γβ+1 is a positive integer with even parity. Additionally,
they satisfy the constraint: Pβ+1
b=1 γb = r.
4
Diffusivity for constant and for arbitrary time-dependent
linear and angular speed
In this section we want to illustrate the power of the formulas derived in the previous section,
which hold for any time-dependent linear and angular speed, to compute the statistics in several
specific cases. In particular, we compare the results obtained by using our formulas with the
results that it is possible to obtain by running Monte Carlo simulations. We generate many trials
of random motion via Monte Carlo simulations according to the motion model in (3). The initial
configuration is the same for all the trials and it is x(0) = y(0) = θ(0) = 0. The final curve
length is also the same and it is s = 1. We consider several values for the two parameters Kθ and
Kr and we set µ(s) constant in 4.1 and variable in 4.2. We will see that, as the number of trials
increases, the statistical properties of the motion directly obtained from the trials approach the
ones obtained by using our formulas. In particular, we will refer to the motion diffusivity and we
compute the moments up to the forth order. We start by computing for each trial, the square of
the final distance from the origin, i.e. the distance at s = 1:
D(1)2 = x(1)2 + y(1)2
(19)
Once we have the previous quantity for many trials, we compute its mean value and its variance.
On the other hand, our formulas allow us to directly obtain both the mean value and the
variance. Specifically, equation (15) provides the mean value. Regarding the variance we have:
σD(s)2 ≡

D(s)4
−

D(s)22. Additionally, D(s)4 = x(s)4 + 2x(s)2y(s)2 + y(s)4 = u(s)2w(s)2,
where we used x(s) = u(s)+w(s)
2
and y(s) = u(s)−w(s)
2i
. The quantity

u(s)2w(s)2
can be easily
computed by using equation (17) with p = q = 2 (see appendix C for the details). Its analytical
expression includes the following six terms: n = 0, l = 0, m = 0 (C000); n = 1, l = 0, m = 0
(C100); n = 1, l = 1, m = 1 (C111); n = 1, l = 2, m = 0 (C120); n = 2, l = 2, m = 0 (C220) and
n = 2, l = 2, m = 2 (C222). We have:

D(s)4
= C000 + C100 + C111 + C120 + C220 + C222
(20)
See appendix C for the analytical expression of the previous terms. In the next two subsections,
we compare the results for the mean value and the variance of D(1)2 obtained by running many
trials and by using our analytic expressions. In 4.1 we refer to a constant ratio between the linear
and the angular speed while in 4.2 we refer to a generic varying ratio. Note that this analysis is
a test for the 4th order statistics of the considered Brownian motion.
4.1
Constant ratio between linear and angular speed
When µ(s) = µ0 we have that θ(s) is linear in s:
θ(s) = θ0 + µ0s
(21)
and, the corresponding active trajectory, is a circumference with radius
1
µ0 (note that when the
angular speed vanishes, the radius becomes infinite and the active trajectory becomes a straight
line). The computation of the integrals in (17) and (18) is immediate. As a result, the expression
9
of

D(s)2
can be easily provided in terms of µ0, s, Kr and Kθ. Let us introduce the following
complex quantity:
z ≡−Kθ
2 + iµ0
(22)
By a direct analytical computation of the double integral in (15) we easily obtain:

D(s)2
= Krs + 2ℜ
ezs −1 −zs
z2

(23)
where the symbol ℜ{·} denotes the real part of a given complex quantity.
We do not provide here the analytical expression for

D(s)4
. We only remark that all the
integrals appearing in the expressions in appendix C can be computed in a similar manner and
are expressed as the real part of suitable analytic functions of the following complex quantities:
z, −3Kθ
2
+ iµ0 and −2Kθ + 2iµ0.
Figure 1: Trajectories generated in the case of constant ratio between linear and angular speed.
The blue line is the trajectory generated without noise (i.e, Kθ = Kr = 0) and the other lines
are five trajectories obtained by setting Kθ = Kr = 0.01m.
In the following, we illustrate some numerical results obtained by setting µ(s) = µ0 = 5.
We consider two cases of Brownian noises, the former is characterized by Kθ = Kr = 0.01m
and the latter by Kθ = Kr = 1m. Figure 1 displays the trajectory generated without noise
(blue line) together with five trajectories obtained by setting Kθ = Kr = 0.01m.
For the
trajectory without noise, the final D(1)2 is equal to 0.0573 m2. We use the expression in (23)
to compute the mean value and a similar expression, which depends on the three mentioned
complex quantities z, −3Kθ
2
+ iµ0 and −2Kθ + 2iµ0, to compute the variance. We obtain the
10
Kθ = Kr (m)
trials

D(1)2
(m2)
σ2
D(1)2 (m4)
0.01
103
0.0688
0.0011
0.01
104
0.0664
0.0012
0.01
105
0.0679
0.0012
1
103
1.1189
1.3736
1
104
1.1124
1.2767
1
105
1.1126
1.3035
Table 1: Constant ratio between linear and angular speed. Values of

D(1)2
and σ2
D(1)2 obtained
by running 103, 104 and 105 trials.
following values:

D(1)2
= 0.0680 m2 and σ2
D(1)2 = 0.0012 m4 when Kθ = Kr = 0.01m and

D(1)2
= 1.1130 m2 and σ2
D(1)2 = 1.3052 m4 when Kθ = Kr = 1m.
We ran 105 trials and we computed from them the mean value and the variance. Figure 2
displays the values of D(1)2 obtained in the first 104 trials together with the value of D(1)2 for
the trajectory without noise (red line), the value of

D(1)2
obtained by averaging D(1)2 on
the 104 trials (green line) and the value of

D(1)2
obtained by our formulas (black line). The
image on the left refers to the case Kθ = Kr = 0.01m and the image on the right to the case
Kθ = Kr = 1m.
Figure 2: Constant ratio between linear and angular speed. Values of D(1)2 for 104 trials together
with the value of D(1)2 for the trajectory without noise (red line), the value of

D(1)2
obtained
by averaging D(1)2 on the trials (green line) and the value of

D(1)2
obtained by our formulas
(black line). Kθ = Kr = 0.01m and Kθ = Kr = 1m in the left and right image, respectively.
Table 1 reports the values of

D(1)2
and σ2
D(1)2 obtained by running 103, 104 and 105 trials.
It is possible to see that, as the number of trials increases, the results converge to the values
obtained by using our formulas.
11
Figure 3: Trajectories generated in the case of varying ratio between linear and angular speed.
The blue line is the trajectory generated without noise (i.e, Kθ = Kr = 0) and the other lines
are five trajectories obtained by setting Kθ = Kr = 0.01m.
12
4.2
Arbitrary time-dependent ratio between linear and angular speed
We consider now a generic function µ(s). Specifically, we consider several kinds of dependence
on s obtaining very similar results. In the following, we illustrate the results obtained by setting
µ(s) = 10 s. As in the previous case, we consider the two cases of Brownian noises, i.e., charac-
terized by Kθ = Kr = 0.01m and by Kθ = Kr = 1m. Figure 3 displays the trajectory generated
without noise (blue line) together with five trajectories obtained by setting Kθ = Kr = 0.01m.
For the trajectory without noise, the final D(1)2 is equal to 0.1021 m2. In this case we directly
use the expressions in (15) in (20) and in appendix C by numerically computing all the inte-
grals. We obtain the following values:

D(1)2
= 0.1124 m2 and σ2
D(1)2 = 0.0026 m4 when
Kθ = Kr = 0.01m and

D(1)2
= 1.1443 m2 and σ2
D(1)2 = 1.4681 m4 when Kθ = Kr = 1m.
We ran 105 trials and we computed from them the mean value and the variance. Figure 4
displays the values of D(1)2 obtained in the first 104 trials together with the value of D(1)2 for
the trajectory without noise (red line), the value of

D(1)2
obtained by averaging D(1)2 on
the 104 trials (green line) and the value of

D(1)2
obtained by our formulas (black line). The
image on the left refers to the case Kθ = Kr = 0.01m and the image on the right to the case
Kθ = Kr = 1m.
Figure 4: Varying ratio between linear and angular speed. Values of D(1)2 for 104 trials together
with the value of D(1)2 for the trajectory without noise (red line), the value of

D(1)2
obtained
by averaging D(1)2 on the trials (green line) and the value of

D(1)2
obtained by our formulas
(black line). Kθ = Kr = 0.01m and Kθ = Kr = 1m in the left and right image, respectively.
Table 2 reports the values of

D(1)2
and σ2
D(1)2 obtained by running 103, 104 and 105 trials. As
in the case of constant ratio between linear and angular speed, we remark that, as the number
of trials increases, the results converge to the values obtained by using our formulas.
13
Kθ = Kr (m)
trials

D(1)2
(m2)
σ2
D(1)2 (m4)
0.01
103
0.1139
0.0022
0.01
104
0.1107
0.0026
0.01
105
0.1123
0.0026
1
103
1.1094
1.2534
1
104
1.1605
1.5289
1
105
1.1494
1.4708
Table 2: Varying ratio between linear and angular speed. Values of

D(1)2
and σ2
D(1)2 obtained
by running 103, 104 and 105 trials.
5
Conclusion
In this paper we derived the statistics, up to any order, for 2D−Brownian unicycle dynamics.
The chosen kinematic constraint is modelled by the unicycle differential equation which is a very
general constraint for a 2D motion. According to this model, the mobile robot can freely rotates
but only one direction for the shift is allowed. This model is suitable to characterize the dynamics
of many wheeled robots: namely, the ones that satisfy the unicycle dynamics. Additionally, the
considered deterministic motion is very general. The expressions here provided for the statistics
hold for any deterministic trajectory that satisfies the mentioned kinematic constraint.
The
expressions contain multiple integrals over the deterministic trajectory.
We showed the power of the derived analytical results by investigating the diffusivity of the
considered Brownian motion for constant and arbitrary time-dependent linear and angular speed.
We want to remark that this paper provides the analytical expressions for the statistics, up
to any order, of a non-trivial Brownian motion, i.e. propelled by arbitrary force and torque. To
the best of our knowledge, analytical solutions for any order moment in the case of arbitrary
linear and angular speed have never been derived in the past.
14
A
Computation of ⟨u(s)p w(s)q⟩
We have:
⟨u(s)p w(s)q⟩= lim
N→∞
X
j1···jpk1···kq
D
ei(θj1+···+θjp−θk1−···−θkq )E
(24)

(δs + ϵj1) · · · (δs + ϵjp)(δs + ϵk1) · · · (δs + ϵkq)

where each index goes from 1 to N. Let us focus our attention on the quantity

(δs + ϵj1) · · · (δs + ϵjp)(δs + ϵk1) · · · (δs + ϵkq)

We can write this quantity as the sum of p+q+1 terms, i.e., δsp+qC0+δsp+q−1C1+δsp+q−2C2+
· · · + δsCp+q−1 + Cp+q. We trivially have C0 = 1. Concerning the remaining coefficients, we
remark, first of all, that only the ones with even index are different from 0. In other words,
C2n+1 = 0, n = 0, 1, · · · , ⌊p+q
2 ⌋(where ⌊x⌋is the largest integer not greater than x). Let us
compute C2n. This coefficient is the sum of
p+q
2n


elements, each of them being the average
of a product of 2n terms "ϵ". On the other hand, since the exponential in (24) is symmetric
with respect to the change ja ↔ja′, a, a′ = 1, · · · , p and with respect to the change kb ↔kb′,
b, b′ = 1, · · · , q but it is not symmetric with respect to the change ja ↔kb, we need to separate
the
p+q
2n


elements in several groups. Specifically, let us denote by l the number of ϵ with index
of type k. The number of elements belonging to this group is
p
2n−l

q
l


. Note that we set
a
b


= 0
when b > a.
Note also that P2n
l=0
p
2n−l

q
l


=
p+q
2n


.
We now remark that, because of the
statistical properties of ϵ, the average of the product of 2n terms ϵ is different from zero only
when their indexes are equal two-by-two. Hence, we have to consider all the combinations of
products of terms ϵ, whose indexes are equal two-by-two and which differ for at least one pair
of indexes. On the other hand, for a given element in the group characterized by a given n and
l, the effect in (24) depends on the number of pairs which are hetero (i.e., with an index of type
j and one of type k). Let us denote by m the number of pairs which are hetero. Note that m
has the same parity of l. Indeed, by definition we have l indexes of type k. Additionally, we are
using m indexes of type j and m indexes of type k to make m hetero pairs. Hence, l −m indexes
of type k remain and, with them, we have to make l−m
2
homo pairs of type k. Similarly, we have
to make 2n−l−m
2
homo pairs of type j. Hence, l−m
2
must be integer.
We compute the number of combinations of the products of 2n terms ϵ, with l indexes of
type k, where the indexes are equal two-by-two, with m pairs which are hetero and which differ
for at least one pair. We obtain:
2n−l
m

 l
m


m!(2n −l −m −1)!!(l −m −1)!!. Following this, we
can write (24) as follows:
⟨u(s)p w(s)q⟩= lim
N→∞
X
{j}p{k}q
⌊p+q
2
⌋
X
n=0
δsp+q−2n
(25)
2n
X
l=0

p
2n −l
q
l

min(l,2n−l)
X
m=0 odd/even
2n −l
m
 l
m

m!
(2n −l −m −1)!!(l −m −1)!!

ϵ2
j1

δj1j2

ϵ2
j3

δj3j4 · · ·
D
ϵ2
j2n−l−m−1
E
δj2n−l−m−1j2n−l−m

ϵ2
k1

δk1k2

ϵ2
k3

δk3k4 · · ·
D
ϵ2
kl−m−1
E
δkl−m−1kl−m
D
ϵ2
kl−m+1
E
δkl−m+1j2n−l−m+1 · · ·
15

ϵ2
kl

δklj2n−l
D
ei(θj1+···+θjp−θk1−···−θkq )E
where, for the brevity sake, we denoted by P
{j}p{k}q the sum PN
j1···jpk1···kq=1. Each average

ϵ2
provides Krδs. Hence, all together, they provide Kn
r δsn. Additionally, the number of Kronecker
deltas is n. Hence, in the limit of N →∞for each value of n we get a multiple integral of
dimension p + q −n. Note that, when the indexes are equal h-by-h (h > 2), the result in the
limit N →∞vanishes since, the power of δs, is larger than the number of sums. By a direct
computation in (25) we obtain:
⟨u(s)p w(s)q⟩=
⌊p+q
2
⌋
X
n=0
Kn
r
2n
X
l=0

p
2n −l
q
l

min(l,2n−l)
X
m=0 odd/even
(26)
2n −l
m
 l
m

m!(2n −l −m −1)!!(l −m −1)!! sm lim
N→∞
X
{js}σ{jd}ρ{ks}χ{kd}η
δsβ D
exp
n
i[θjs
1 + · · · + θjsσ + 2(θjd
1 +
+ · · · + θjd
ρ) −(θks
1 + · · · + θksχ) −2(θkd
1 + · · · + θkdη)]
oE
where ρ ≡2n−l−m
2
, σ ≡p −(2n −l), η ≡l−m
2 , χ ≡q −l and β ≡ρ + σ + η + χ = p + q −n −m.
We must compute the average of the exponential in (26) and then we compute the limit N →∞.
We remark that the various θ in the exponential contain random quantities (i.e., the δθ at
different time steps). In order to proceed we have to separate all the random quantities which are
independent. We start this separation by redefining the indexes in the sum P
{js}σ{jd}ρ{ks}χ{kd}η.
Specifically, we consider the new indexes j
sj
dk
sk
d which differ from jsjdkskd since they are
ordered (in increasing order). For instance, j
s
1 < j
s
2 < · · · < j
s
σ. Hence, the last sum in (26)
can be replaced with P
{js}σ{jd}ρ{ks}χ{kd}η →ρ!σ!η!χ! P
{js}σ{jd}ρ{k
s}χ{k
d}η. The four types of
indexes are not ordered among them. Hence, the sum includes all the possible combinations which
maintain the order only restricted to a single index type. For instance, a possible combination
is: k
d
1 < k
d
2 < j
s
1 < k
s
1 < j
s
2 < k
d
3 < j
d
1 < j
d
2 < · · · . We introduce the β ordered indexes
w1 < w2 < · · · < wβ. Additionally, let us denote with ∆b
a ≡Pb
c=a+1 δθc and ∆
b
a ≡Pb
c=a+1 δθc =
θ(bδs)−θ(aδs), where θ(s) is the deterministic trajectory. Finally, let us define α ≡2ρ+σ−χ−2η.
The sum in the exponential in (26) contains α(∆w1
0 +∆
w1
0 ). Then, depending on the combination
of the indexes j
sj
dk
sk
d, we have a different result for the (∆w2
w1+∆
w2
w1). Specifically, if w1 = j
d
1 the
sum in the exponential contains (α −2)(∆w2
w1 +∆
w2
w1). If w1 = j
s
1 it contains (α −1)(∆w2
w1 +∆
w2
w1).
If w1 = k
s
1 it contains (α + 1)(∆w2
w1 + ∆
w2
w1). If w1 = k
d
1 it contains (α + 2)(∆w2
w1 + ∆
w2
w1). We
introduce the vector c ≡[c1, · · · , cβ]T whose entries are −2, −1, 1 and 2. Specifically, it contains
η entries equal to 2, χ equal to 1, σ equal to −1 and ρ equal to −2. It is easy to realize that we
have
β
χ

β−χ
σ

ρ+η
η


vectors c. Finally, we define Φb ≡Pb
a=1 ca. Note that Φβ = −α. According
to this, the last sum in (26) can be written as follows:
eiαθ0ρ!σ!η!χ!
X
c
X
{w}β
δsβ·
·
*
exp
(
i
β−1
X
b=0
h
(α + Φb)(∆wb+1
wb
+ ∆
wb+1
wb
)
i)+
16
Hence, we have:
⟨u(s)p w(s)q⟩=
⌊p+q
2
⌋
X
n=0
Kn
r
2n
X
l=0

p
2n −l
q
l

·
(27)
·
min(l,2n−l)
X
m=0 odd/even
2n −l
m
 l
m

m!(2n −l −m −1)!!(l −m −1)!!·
·smeiαθ0ρ!σ!η!χ!
X
c
lim
N→∞
X
{w}β
δsβ·
·
*
exp
(
i
β−1
X
b=0
h
(α + Φb)(∆wb+1
wb
+ ∆
wb+1
wb
)
i)+
Now we can compute the average since we were able to separate all the independent quantities.
We have:
*
exp
(
i
β−1
X
b=0
h
(α + Φb)(∆wb+1
wb
+ ∆
wb+1
wb
)
i)+
=
(28)
= exp
(
i
β−1
X
b=0
h
(α + Φb)∆
wb+1
wb
i)
·
·
β−1
Y
b=0
D
exp
n
i
h
(α + Φb)∆wb+1
wb
ioE
=
= exp
(
i
β−1
X
b=0
h
(α + Φb)∆
wb+1
wb
i)
·
·
β−1
Y
b=0
exp

−(α + Φb)2(wb+1 −wb)Kθδs
2

=
e
Pβ−1
b=0

i(α+Φb)∆
wb+1
wb
−
(α+Φb)2(wb+1−wb)Kθδs
2

We used the equality (η = N (0, 1)):

eAη
= e
A2
2
(29)
with A = i(α+Φb)σb+1 and σ2
b+1 = (wb+1 −wb)Kθδs. By substituting equation (28) in (27) and
by taking the limit (N →∞), we finally obtain:
⟨u(s)p w(s)q⟩=
⌊p+q
2
⌋
X
n=0
Kn
r
2n
X
l=0

p
2n −l
q
l

·
(30)
·
min(l,2n−l)
X
m=0 odd/even
2n −l
m
 l
m

m!(2n −l −m −1)!!(l −m −1)!!·
17
·smeiαθ0ρ!σ!η!χ!
X
c
Z s
0
ds1
Z s
s1
ds2 · · ·
Z s
sβ−1
dsβ
e
Pβ−1
b=0

i(α+Φb)[θ(sb+1)−θ(sb)]−
(α+Φb)2(sb+1−sb)Kθ
2

which coincides with (17).
B
Computation of
D
u(s)p w(s)q ˜θrE
This computation follows the same initial steps carried out in appendix A. We obtain an ex-
pression equal to the one given in (27) but the mean value at the end must be replaced with
D
exp
n
i Pβ−1
b=0
h
(α + Φb)(∆wb+1
wb
+ ∆
wb+1
wb
)
io
· ˜θrE
. The deterministic part can be factorized out
the mean value. We need to calculate:
*
exp
(
i
β−1
X
b=0
h
(α + Φb)∆wb+1
wb
i)
· ˜θr
+
=
=
*
˜θr
β−1
Y
b=0
n
exp
h
i(α + Φb)∆wb+1
wb
io+
By using the multinomial theorem we can write:
˜θr =
 β
X
b=0
∆wb+1
wb
!r
=
X
γ
r! ·
β
Y
b=0

∆wb+1
wb
γb+1
γb+1!
where wβ+1 ≡N and the sum over γ is the sum over all the vectors γ = [γ1, · · · , γβ+1], where γb
(b = 1, · · · , β + 1) are positive integers satisfying the constraint: Pβ+1
b=1 γb = r. Hence, we obtain
*
˜θr
β−1
Y
b=0

ei(α+Φb)∆
wb+1
wb
+
=
(31)
X
γ
r!
*
∆wβ+1
wβ
γβ+1
γβ+1!
β−1
Y
b=0
ei(α+Φb)∆
wb+1
wb

∆wb+1
wb
γb+1
γb+1!
+
=
=
X
γ
r!
Qβ
b=0 γb+1!
D
∆wβ+1
wβ
γβ+1E
·
·
β−1
Y
b=0

ei(α+Φb)∆
wb+1
wb

∆wb+1
wb
γb+1
We use the following two standard results for a normal distribution:
D
∆wβ+1
wβ
γβ+1E
=

0
γβ+1
odd
σγβ+1
β+1 (γβ+1 −1)!!
γβ+1
even
(32)
and
D
exp
h
i(α + Φb)∆wb+1
wb
i 
∆wb+1
wb
γb+1E
=
(33)
18
= σγb+1
b+1 exp

−(α + Φb)2σ2
b+1
2

·
·
⌊
γb+1
2
⌋
X
a=0
γb+1! (i(α + Φb)σb+1)γb+1−2a
a!(γb+1 −2a)!2a
with σ2
b+1 ≡Kθδs(wb+1−wb). Equation (33) is obtained starting from (29) and by differentiating
γb+1 times with respect to A. Hence, the expression in (31), becomes:
*
˜θr
β−1
Y
b=0
n
exp
h
i(α + Φb)∆wb+1
wb
io+
=
(34)
=
X
γ
r!σγβ+1
β+1 (γβ+1 −1)!!
Qβ
b=0 γb+1!
β−1
Y
b=0

σγb+1
b+1 e−
(α+Φb)2σ2
b+1
2
·
·
⌊
γb+1
2
⌋
X
a=0
γb+1! (i(α + Φb)σb+1)γb+1−2a
a!(γb+1 −2a)!2a



where now the sum over γ only includes the vectors γ whose last entry is even. By taking the
limit N →∞, this expression becomes:
X
γ
r!(γβ+1 −1)!!
Qβ
b=0 γb+1!
(s −sβ)
γβ+1
2
K
r
2
θ ·
·
β−1
Y
b=0

(sb+1 −sb)
γb+1
2 e−
(α+Φb)2Kθ(sb+1−sb)
2
·
·
⌊
γb+1
2
⌋
X
a=0
γb+1!

i(α + Φb)
p
Kθ(sb+1 −sb)
γb+1−2a
a!(γb+1 −2a)!2a





and we finally obtain:
D
u(s)p w(s)q ˜θ(s)rE
=
⌊p+q
2
⌋
X
n=0
Kn
r
2n
X
l=0

p
2n −l
q
l

·
(35)
·
min(l,2n−l)
X
m=0 odd/even
2n −l
m
 l
m

m!(2n −l −m −1)!!(l −m −1)!!
·smeiαθ0ρ!σ!η!χ!
X
c
X
γ
r!(γβ+1 −1)!!K
r
2
θ
Qβ
b=0 γb+1!
Z s
0
ds1
Z s
s1
ds2·
· · ·
Z s
sβ−1
dsβ(s −sβ)
γβ+1
2
ei Pβ−1
b=0 (α+Φb)[θ(sb+1)−θ(sb)]
β−1
Y
b=0

(sb+1 −sb)
γb+1
2 e−
(α+Φb)2Kθ(sb+1−sb)
2
·
19
⌊
γb+1
2
⌋
X
a=0
γb+1!

i(α + Φb)
p
Kθ(sb+1 −sb)
γb+1−2a
a!(γb+1 −2a)!2a





which coincides with (18).
C
Computation of

D(s)4
=

u(s)2 w(s)2
We use equation (17) with p = q = 2 to obtain the analytical expression of

D(s)4
=

u(s)2w(s)2
.
The first sum on n includes the three terms n = 0, 1, 2.
When n = 0, the
second sum on l and the third sum on m only include one term, i.e., l = m = 0. All the factors
before the sum on c are 1 with the exception of (p −2n + l)! = 2 and (q −l)! = 2. Additionally,
β = 4 and the sum on c consists of six vectors, i.e., [−1, −1, 1, 1], [−1, 1, −1, 1], [−1, 1, 1, −1],
[1, −1, −1, 1], [1, −1, 1, −1], [1, 1, −1, −1]. Hence we obtain the following term:
C000 = 8
Z s
0
ds1
Z s
s1
ds2
Z s
s2
ds3
Z s
s3
ds4
n
e−
Kθ(−s1−3s2+3s3+s4)
2
cos
−θ1 −θ2 + θ3 + θ4


+
+e−
Kθ(−s1+s2−s3+s4)
2

cos
−θ1 + θ2 −θ3 + θ4


+ cos
−θ1 + θ2 + θ3 −θ4

o
where θj ≡θ(sj), j = 1, 2, 3, 4.
Let us consider now the terms with n = 1. In this case, l = 0, 1, 2. For each value of l we
have in general several values of m, which must have the same parity of l. Hence, when l = 1,
we only have m = 1. When l = 2, we can have m = 0, 2. On the other hand, when l = m = 2,
the factor
2n−l
m


=
0
2


vanishes. By an explicit computation, it is possible to see that C120 is
the conjugate of C100. Hence, their sum is real and it is:
C100 + C120 = 4Kr
Z s
0
ds1
Z s
s1
ds2
Z s
s2
ds3
n
e−
Kθ(−4s1+3s2+s3)
2
cos
2θ1 −θ2 −θ3


+
+e−
Kθ(−s1+s3)
2
cos
−θ1 + 2θ2 −θ3


+ e−
Kθ(−s1−3s2+4s3)
2
cos
−θ1 −θ2 + 2θ3

o
The last term with n = 1 is the one with l = m = 1. By a direct computation we obtain:
C111 = 8Krs
Z s
0
ds1
Z s
s1
ds2e−
Kθ(−s1+s2)
2
cos
−θ1 + θ2


Finally, when n = 2 the sum on l includes the five values: l = 0, 1, 2, 3, 4. On the other hand,
only the term l = 2 does not vanish since the product
p
2n−l

q
l


is zero in all the other cases.
When l = 2 we have two possible values for m, i.e., m = 0, 2. We obtain the following two
contributions:
C220 = 2K2
r
Z s
0
ds1
Z s
s1
ds2e−2Kθ(−s1+s2) cos
−2θ1 + 2θ2


C222 = 2K2
rs2
20
References
[1] Sepideh Bazazi, Pawel Romanczuk, Sian Thomas, Lutz Schimansky-Geier, Joseph J. Hale,
Gabriel A. Miller, Gregory A. Sword, Stephen J. Simpson and Iain D. Couzin, Nutritional
state and collective motion: from individuals to mass migration, Proc. R. Soc. B 278, 356
(2010).
[2] Callan, David (2009), A combinatorial survey of identities for the double factorial,
arXiv:0906.1317
[3] S. Challa and Y. Bar-Shalom, "Nonlinear filter design using Fokker-Planck-Kolmogorov prob-
ability density evolutions", IEEE Transactions on Aerospace and Electronic Systems, vol. 36,
no. 1, 2000
[4] I. D. Couzin et al., Collective memory and spatial sorting in animal groups, J. Theor. Biol.
218, 1 (2002).
[5] F. E. Daum, "New nonlinear filters and exact solutions of the Fokker-Planck equation",
Proceedings of IEEE American Control Conference, pp.884 -888 1986
[6] F. E. Daum, "Solution of the Zakai equation by separation of variables", IEEE Transactions
on Automatic Control, vol. AC-32, no. 10, pp.941 -943 1987
[7] F. E. Daum, "Nonlinear filters: beyond the Kalman filter", Aerospace and Electronic Systems
Magazine, IEEE, Volume:20, Issue: 8, 2005
[8] R. Dreyfus, J. Baudry, M. L. Roper, M. Fermigier, H. A. Stone, and J. Bibette, Microscopic
artificial swimmers, Nature London, 437, 862, 2005
[9] Freidlin Mark and Hu Wenqing, On Diffusion in Narrow Random Channels, J Stat Phys,
152:136–158, 2013
[10] Hanggi Peter, Marchesoni Fabio, Artificial Brownian motors: Controlling transport on the
nanoscale, Reviews of Modern Physics 81 (1): 387, 2009
[11] J. C. Latombe, Robot motion planning, Kluwer Academic Publishers, 1991
[12] E. Lauga and T. R. Powers, The hydrodynamics of swimming microorganisms, Rep. Prog.
Phys. 72, 096601, 2009
[13] M. Leoni, J. Kotar, B. Bassetti, P. Cicuta, and M. C. Lagomarsino, A basic swimmer at low
Reynolds number, Soft Matter 5, 472, 2009
[14] Long, A.W., Wolfe, K.C., Mashner, M.J., Chirikjian, G.S., “The Banana Distribution is
Gaussian: A Localization Study with Exponential Coordinates,” Robotics: Science and Sys-
tems, Sydney, NSW, Australia, July 09 - July 13, 2012.
[15] A. Martinelli, The odometry error of a mobile robot with a synchronous drive system, IEEE
Transactions on Robotics and Automation, Volume 18, Issue 3, 399–405, Jun 2002.
[16] A. Martinelli, The accuracy on the parameter estimation of an odometry system of a mobile
robot, ICRA ’02. IEEE International Conference on Robotics and Automation, 1378–1383
vol.2, 2002.
[17] H. Niwa, Self-organizing dynamic model offish schooling, J. Theor. Biol. 171, 123 (1994).
21
[18] Bernt Oksendal, Stochastic Differential Equations:
An Introduction with Applications,
Springer, Berlin, 2003
[19] Papoulis, A. Probability, Random Variables, and Stochastic Processes, 2nd ed. New York:
McGraw-Hill, 1984.
[20] Park, W., Kim, J.S., Zhou, Y., Cowan, N.J., Okamura, A.M., Chirikjian, G.S., “Diffusion-
based motion planning for a nonholonomic flexible needle model,” ICRA, Barcelona, Spain,
April 2005
[21] Park, W., Liu, Y., Zhou, Y., Moses, M., Chirikjian, G.S., “Kinematic State Estimation
and Motion Planning for Stochastic Nonholonomic Systems Using the Exponential Map,”
Robotica 26(4), 419-434. July-August 2008
[22] Park, W., Reed, K. B., Okamura, A. M., Chirikjian, G. S., “Estimation of model parameters
for steerable needles,” ICRA, Anchorage, Alaska, 2010.
[23] Reimann Peter, Brownian motors: noisy transport far from equilibrium, Physics Reports,
361, 57–265, 2002
[24] Hannes Risken, The Fokker-Planck Equation Methods of Solution and Applications,
Springer Series in Synergetics Volume 18, 1989
[25] A. Sokolov, I. S. Aranson, J. O. Kessler, and R. E. Goldstein, Concentration Dependence of
the Collective Dynamics of Swimming Bacteria, Phys. Rev. Lett. 98, 158102
[26] P. Tierno, R. Golestanian, I. Pagonabarraga, and F. Sagues, Magnetically actuated colloidal
microswimmers, J. Phys. Chem. B. 112, 16525 2008
[27] P Tierno, R Golestanian, I Pagonabarraga, F Sagués, Controlled swimming in confined
fluids of magnetically actuated colloidal rotors, Phys Rev Lett 101 218304, 2008
[28] R. Trouilloud, T. S. Yu, A. E. Hosoi and E. Lauga, Soft swimming: Exploiting deformable
interfaces for low-Reynolds number locomotion, Phys. Rev. Lett. 101, 048102, 2008
[29] Turner Jonathan, Brownian Motion Applied to Human Intersections, thesis at Texas State
University-San Marcos, Dept. of Mathematics, 2012
[30] Zhou, Y., Chirikjian, G.S., “Probabilistic Models of Dead-Reckoning Error in Nonholonomic
Mobile Robots,” ICRA’03, Taipei, Taiwan, September, 2003.
[31] Zhou Y., Chirikjian G.S., Planning for Noise-Induced Trajectory Bias in Nonholonomic
Robots with Uncertainty, 2004 IEEE International Conference on Robotics and Automation
(ICRA 2004), pp.4596-4601, New Orleans, Louisiana, 2004.
22
Complete analytic solution to Brownian unicycle dynamics
Agostino Martinelli
October 8, 2018
Abstract
This paper derives a complete analytical solution for the probability distribution of the
configuration of a non-holonomic vehicle that moves in two spatial dimensions by satisfying
the unicycle kinematic constraints and in presence of Brownian noises. In contrast to pre-
vious solutions, the one here derived holds even in the case of arbitrary linear and angular
speed. This solution is obtained by deriving the analytical expression of any-order moment
of the probability distribution. To the best of our knowledge, an analytical expression for
any-order moment that holds even in the case of arbitrary linear and angular speed, has
never been derived before. To compute these moments, a direct integration of the Langevin
equation is carried out and each moment is expressed as a multiple integral of the determin-
istic motion (i.e., the known motion that would result in absence of noise). For the special
case when the ratio between the linear and angular speed is constant, the multiple integrals
can be easily solved and expressed as the real or the imaginary part of suitable analytic
functions. As an application of the derived analytical results, the paper investigates the
diffusivity of the considered Brownian motion for constant and for arbitrary time-dependent
linear and angular speed.
Keywords: Brownian motion, Stochastic dynamics, Fokker-Plank equation without detailed
balance, Stochastic processes, Localization.
1
arXiv:1501.03300v1  [cs.RO]  14 Jan 2015
Contents
1
Introduction
3
2
Stochastic motion model
4
3
Computation of the probability distribution
5
3.1
First and second-order statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
3.2
Computation of any-order moment . . . . . . . . . . . . . . . . . . . . . . . . . .
7
4
Diffusivity for constant and for arbitrary time-dependent linear and angular
speed
9
4.1
Constant ratio between linear and angular speed . . . . . . . . . . . . . . . . . .
9
4.2
Arbitrary time-dependent ratio between linear and angular speed . . . . . . . . .
13
5
Conclusion
14
A Computation of ⟨u(s)p w(s)q⟩
15
B Computation of
D
u(s)p w(s)q ˜θrE
18
C Computation of

D(s)4
=

u(s)2 w(s)2
20
2
1
Introduction
In recent years a great interest has been devoted to investigate the statistical properties of the
motion of active (or self-propelled) particles. These particles differ from passive particles since
they move under the action of an internal force. Typical examples of these particles are Brownian
motors [9, 10, 23], biological and artificial microswimmers [8, 12, 13, 25, 26, 27, 28], macroscopic
animals [1, 4, 17] and even pedestrians [29].
In the last decade, this investigation has also been addressed in the framework of mobile
robotics, in particular by analyzing the 2D cases of the unicyle, cart, and car [30, 31], and the
3D case of flexible needle steering [20, 22]. In these works, analytical solutions are obtained by
directly solving the corresponding Fokker-Planck equations by using the Fourier transform. An-
other method proposed for obtaining a closed-form solution is the use of exponential coordinates
[14, 21]. In all these works, the benefit of having closed-form solutions is clearly illustrated by
introducing motion planning methods based on them.
In this paper we want to deal with wheeled robots, which represent a special case of active
particles moving in a two dimensional space. Their configuration is characterized by a 3D-vector
(2 components characterize their position and one their orientation). On the other hand, in most
of cases, a wheeled robot cannot move along any direction since it must satisfy the so-called
non-holonomic constraints [11]. Our stochastic motion model (presented in section 2) differs
from the one considered in [14, 21, 30, 31] since we assume independent errors acting on the
linear and angular components. The model adopted in [14, 21, 30, 31] refers to the case of a
differential drive system where the errors acting on each wheel are independent.
Starting from our stochastic motion model, we compute the analytical expression of any-
order moment of the probability distribution (section 3). Specifically, each moment is expressed
as a multiple integral of the deterministic motion performed by the mobile robot (i.e., the known
motion that would result in absence of noise). In other words, the expression is a multiple integral
on two time-dependent functions, which describe the time behaviour of the deterministic linear
and angular speed. This allows obtaining any order moment for any deterministic motion, i.e.,
when the mobile robot is propelled by arbitrary time-dependent linear and angular speed. For
the special case when the ratio between these two speeds is constant, the multiple integrals can be
expressed as the real or the imaginary part of suitable analytic functions. In section 4 we show the
power of the derived analytical results by investigating the diffusivity of the considered Brownian
motion for constant and for arbitrary time-dependent linear and angular speed. Conclusions are
provided in section 5.
The interest in deriving the statistical properties of the motion of a wheeled robot comes
from the possibility of improving its localization (both in precision and speed). The localization
is a fundamental problem in mobile robotics which must be solved to autonomously and safely
navigate. It is a non linear estimation problem. In most of cases, the localization is carried out
by using nonlinear filters (e.g., Extended and Unscented Kalman Filters, Particle filters, etc). On
the other hand, non linear filters are strictly connected with the solution of partial differential
equations (see [3, 5, 6, 7] and the next section 2). Very easily speaking, for a dynamic stochastic
process described by stochastic differential equations, it is possible to introduce a probability
distribution which provides the probability that the process takes a given set of values. This
probability distribution must satisfy a partial differential equation. In this sense, a non linear
filter could be considered as a numerical solution of this partial differential equation. In the
specific case of localization, the dynamic process is the motion of the robot together with the
noisy data delivered by its sensors.
We already computed the statistics up to the second order for a similar motion model [15,
16] (in the expressions in [15] there are some typos, corrected in [16]). Here, we extend the
3
computation in [15] by including any-order statistics. We also show that the expressions hold for
a much more general model (the one described in section 2) and not only for a simple specific
odometry error model (as shown in [15]). As a result, this paper provides a complete analytical
solution to a Fokker-Planck equation for which the detailed balance condition is not verified1.
Note that the goal of this paper is to analytically compute the moments of the probability
distribution up to any order for the configuration of a non-holonomic vehicle. As it will be seen,
this is a very hard task from a computational point of view. These analytical results could be
very useful to improve the localization which is, as previously mentioned, a numerical solution
of a partial differential equation.
2
Stochastic motion model
We consider a mobile robot which moves in a 2D−environment. The configuration of the mobile
robot is characterized by its position and orientation. In Cartesian coordinates, we denote them
by r ≡[x, y]T and θ, respectively. We assume that the motion of the mobile robot satisfies the
unicycle model [11], which is the most general nonholonomic model for 2D−motions, where the
shift of the mobile robot can only occur along one direction ([cos θ, sin θ]T ):


˙x =
v cos θ
˙y =
v sin θ
˙θ =
ω
(1)
The quantities v = v(t) and ω = ω(t) are the linear and angular speed, respectively. We assume
that these functions are known and we call the motion that results from them the deterministic
motion. Now we want to introduce a stochastic model that generalizes equation (1) by accounting
a Gaussian white noise. We assume that also the noise satisfies the same nonholonomic constraint.
The stochastic differential equations are:


dx(t) =
cos θ(t) [v(t)dt + α(t)dwr(t)]
dy(t) =
sin θ(t) [v(t)dt + α(t)dwr(t)]
dθ(t) =
ω(t)dt + β(t)dwθ(t)
(2)
where [dwr(t), dwθ(t)]T is a standard Wiener process of dimension two [18].
The functions
α(t) and β(t) are modeled according to the following physical requirement (diffusion in the
overdamped regime). We require that, when the mobile robot moves during the infinitesimal
time interval dt, the shift is a random Gaussian variable, whose variance increases linearly with
the traveled distance. This is obtained by setting α(t) =
p
Kr|v(t)|, with Kr a positive parameter
which characterizes our system (interaction robot-environment). Similarly, we require that the
rotation accomplished during the same time interval dt is a random Gaussian variable, whose
variance increases linearly with the traveled distance. Hence, we set β(t) =
p
Kθ|v(t)|, with Kθ
another positive parameter which characterizes our system (interaction robot-environment). The
equations in (2) are the Langevin equations in the overdamped limit. From now on, we assume
that the mobile robot can only move ahead2, i.e., v(t) ≥0. We remind the reader that v(t) is a
deterministic and known function of time. The curve length ds = v(t)dt is the curve length of the
deterministic motion, which is independent of the Wiener process [dwr(t), dwθ(t)]T . Hence, we
1Note that, deriving any-order moment, corresponds to analytically derive the probability density distribution
[19].
2Note that this constraint does not limit the robot motion: for instance, a go and back motion is easily obtained
by accomplishing 180deg rotations
4
are allowed to use the deterministic curve length ds = v(t)dt instead of the time t. By denoting
the ratio between the angular and the linear speed by µ, equations in (2) read:


dx(s) =
cos θ(s) [ds + K
1
2r dwr(s)]
dy(s) =
sin θ(s) [ds + K
1
2r dwr(s)]
dθ(s) =
µ(s)ds + K
1
2
θ dwθ(s)
(3)
The associated Smoluchowski equation is:
∂p
∂s = −∇· (D1p) + ∇· (D2∇p) = 0
(4)
where:
• p = p(x, y, θ; s(t)) is the probability density for the mobile robot at the configuration
(x, y, θ) and at time t;
• D1 = [cos θ, sin θ, µ]T is the drift vector;
• D2 = 1
2


Kr cos2 θ
Kr cos θ sin θ
0
Kr cos θ sin θ
Kr sin2 θ
0
0
0
Kθ

is the diffusion tensor.
Our goal is to obtain the probability density p(x, y, θ; s). This will be done in the next two
sections.
3
Computation of the probability distribution
The probability density p(x, y, θ; s) satisfies the Smoluchowski partial differential equation in
(4), which is a special case of the Fokker-Planck equation. Since the detailed balance condition
is not satisfied [24], we follow a different procedure to compute p(x, y, θ; s). Specifically, we
use the Langevin equation in (3) to compute the moments, up to any order, of the probability
distribution. First of all, we remark that the third equation in (3) is independent of the first two,
is independent of θ and is linear in dwθ. As a result, the probability density only in terms of
θ, i.e., the probability Pθ(θ; s) ≡
R
dx
R
dy p(x, y, θ; s), is a Gaussian distribution with mean
value θ0 +
R s
0 ds′µ(s′) and variance Kθs, where θ0 is the initial orientation. This same result
could also be obtained by integrating (4) in x and y and by using the divergence theorem. In
other words, we have the following distribution for the orientation at a given value of s:
θ(s) = N(θ(s), Kθs)
(5)
where N(·, ·) denotes the normal distribution with mean value and variance the first and the
second argument; θ(s) ≡θ0 +
R s
0 ds′ µ(s′).
Note that, according to the unicycle model, a
trajectory is completely characterized by its starting point and by the orientation vs the curve
length, i.e., by the function θ(s). In the following, we will call the deterministic function θ(s),
the deterministic trajectory.
5
3.1
First and second-order statistics
Let us consider the first equation in (3). We compute the expression of x(s) by a direct integra-
tion. We divide the interval (0, s) in N equal segments, δs ≡
s
N . We have:
x(s) = lim
N→∞
N
X
j=1
(δs + ϵj) cos θj
(6)
where ϵj is a random Gaussian variable satisfying ⟨ϵj⟩= 0, ⟨ϵj ϵk⟩= δjkKrδs for j, k = 1, · · · , N
(δjk is the Kronecker delta) and θj ≡θ(jδs). On the other hand, we have:
θj = θ(jδs) +
j
X
m=1
δθm ≡θj + ∆θj
(7)
where, according to our stochastic model in (3), δθj is a random Gaussian variable satisfy-
ing ⟨δθj⟩= 0, ⟨δθj ϵk⟩= 0 and ⟨δθj δθk⟩= δjkKθδs for j, k = 1, · · · , N.
As a result,
∆θj ≡Pj
m=1 δθm is also a random Gaussian variable satisfying ⟨∆θj⟩= 0, ⟨∆θj ∆θk⟩=
 Kθjδs
if j ≤k
Kθkδs
if j > k .
Starting from (6) and (7) and by remarking that ⟨cos ∆θj⟩= e−
kθjδs
2
and ⟨sin ∆θj⟩= 0, it is easy to obtain the mean value of x(s). We have:
⟨x(s)⟩= lim
N→∞
N
X
j=1
(δs + ⟨ϵj⟩) ⟨cos θj⟩=
(8)
= lim
N→∞
N
X
j=1
δs cos θj e−
kθjδs
2
=
Z s
0
ds′ cos θ(s′) e−
Kθs′
2
Similarly, it is possible to obtain the mean value of y(s):
⟨y(s)⟩=
Z s
0
ds′ sin θ(s′) e−
Kθs′
2
(9)
Finally, in a similar manner, but with some more computation, we obtain all the second-order
moments:

x(s)2
=
Z s
0
ds′
Z s−s′
0
ds′′e−
Kθs′′
2
{(1 + χc(s′))
(10)
cos[θ(s′ + s′′) −θ(s′)] −χs(s′) sin[θ(s′ + s′′) −θ(s′)]
	
+
+Kr
2

s +
Z s
0
ds′ χc(s′)


y(s)2
=
Z s
0
ds′
Z s−s′
0
ds′′e−
Kθs′′
2
{(1 −χc(s′))
(11)
cos[θ(s′ + s′′) −θ(s′)] + χs(s′) sin[θ(s′ + s′′) −θ(s′)]
	
+
+Kr
2

s −
Z s
0
ds′ χc(s′)

⟨x(s) y(s)⟩=
Z s
0
ds′
Z s−s′
0
ds′′e−
Kθs′′
2
{χs(s′)
(12)
6
cos[θ(s′ + s′′) −θ(s′)] + χc(s′) sin[θ(s′ + s′′) −θ(s′)]
	
+
+Kr
2
Z s
0
ds′ χs(s′)
σxθ(s) ≡⟨x(s) θ(s)⟩−⟨x(s)⟩θ(s) = 2Kθ
∂⟨y(s)⟩
∂Kθ
(13)
σyθ(s) ≡⟨y(s) θ(s)⟩−⟨y(s)⟩θ(s) = −2Kθ
∂⟨x(s)⟩
∂Kθ
(14)
where χc(s′) ≡cos[2θ(s′)] e−2Kθs′ and χs(s′) ≡sin[2θ(s′)] e−2Kθs′.
Obtaining the expression of higher-order moments will demand more tricky computation and
will be dealt separately, in the next subsection. Here we conclude by considering the quantity
D(s)2 ≡x(s)2 + y(s)2, which provides the time-evolution of the square of the distance of the
mobile robot from its initial position. From (10) and (11) we obtain a simple expression for its
mean value:

D(s)2
=
(15)
Krs + 2
Z s
0
ds′
Z s−s′
0
ds′′e−
Kθs′′
2
cos[θ(s′ + s′′) −θ(s′)]
3.2
Computation of any-order moment
We introduce the following two complex random quantities:
u(s) ≡lim
N→∞
N
X
j=1
(δs + ϵj)eiθj; w(s) ≡lim
N→∞
N
X
j=1
(δs + ϵj)e−iθj
(16)
From (6) and (16) it is immediate to realize that x(s) = u(s)+w(s)
2
. Similarly we have y(s) =
u(s)−w(s)
2i
. Hence, in order to compute any-order moment which involves x(s) and y(s) it suffices
to compute ⟨u(s)p w(s)q⟩for any p, q ∈N. The computation of this quantity requires several
tricky steps, which are provided in appendix A. The key is to separate all the independent random
quantities in order to compute their mean values. This is obtained by arranging all the sums in
a suitable manner (see appendix A for all the details). We have:
⟨u(s)p w(s)q⟩=
⌊p+q
2
⌋
X
n=0
Kn
r
2n
X
l=0

p
2n −l
q
l

(17)
X
m
2n −l
m
 l
m

m!(2n −l −m −1)!!(l −m −1)!! sm
ei(p−q)θ0
2n −l −m
2

!
l −m
2

!(p −2n + l)!(q −l)!
X
c
Z s
0
ds1
Z s
s1
ds2 · · ·
Z s
sβ−1
dsβexp
(β−1
X
b=0
[i(p −q + Φb)·
·[θ(sb+1) −θ(sb)] −(p −q + Φb)2(sb+1 −sb)Kθ
2

where:
7
• ⌊p+q
2 ⌋is the largest integer not greater than p+q
2 ;
• the second sum on l (i.e., P2n
l=0) is restricted to the values of l for which 2n −l ≤p and
l ≤q (or, equivalently, we are using the convention that
x
y


= 0 when y > x);
• the sum on m (i.e., P
m) goes from 0 to the minimum between l and 2n −l and it is
restricted to the integers m with the same parity of l (hence, both 2n−l−m
2
and l−m
2
are
integers);
• the symbols "!" and "!!" denote the factorial and the double factorial, respectively [2] (note
that 0! = 0!! = (−1)!! = 1);
• β = p + q −n −m is the dimension of the remaining multiple integral (note that a multiple
integral of dimension m has already been computed and provided the term sm);
• the sum on c, i.e., P
c, is the sum over all the vectors c of dimension β, whose entries
are −2, −1, 1 and 2: specifically, each vector c has l−m
2
entries equal to 2, q −l equal
to 1, p −2n + l equal to −1 and 2n−l−m
2
equal to −2 (note that the sum P
c consists of
 β
q−l

 β−q+l
p−2n+l

n−m
l−m
2


addends);
• Φb ≡Pb
a=1 ca (note that Φβ = q −p);
• s0 ≡0 and θ(s0) = θ0 (note that s0 is not a variable of integration).
In order to complete the derivation of the statistics for our problem, we need to compute any-
order moment which also involves the orientation θ. We provide the formula for the quantity
D
u(s)p w(s)q ˜θ(s)rE
, where ˜θ(s) ≡θ(s) −θ(s). The details of this computation are provided in
appendix B. We have, for any p, q, r ∈N:
D
u(s)p w(s)q ˜θ(s)rE
=
⌊p+q
2
⌋
X
n=0
Kn
r
2n
X
l=0

p
2n −l
q
l

(18)
X
m
2n −l
m
 l
m

m!(2n −l −m −1)!!(l −m −1)!! sm ei(p−q)θ0
2n −l −m
2

!
l −m
2

!(p −2n + l)!(q −l)!
X
c
X
γ
r!(γβ+1 −1)!!K
r
2
θ
Qβ
b=0 γb+1!
Z s
0
ds1
Z s
s1
ds2 · · ·
Z s
sβ−1
dsβ(s −sβ)
γβ+1
2
exp
(
i
β−1
X
b=0
(p −q + Φb)[θ(sb+1) −θ(sb)]
)
·
·
β−1
Y
b=0

(sb+1 −sb)
γb+1
2 exp

−(p −q + Φb)2Kθ(sb+1 −sb)
2

⌊
γb+1
2
⌋
X
a=0
γb+1!

i(p −q + Φb)
p
Kθ(sb+1 −sb)
γb+1−2a
a!(γb+1 −2a)!2a





8
where the sum over γ, i.e. P
γ, is the sum over all the vectors γ = [γ1, γ2, · · · , γβ, γβ+1], where γb
(b = 1, · · · , β) are positive integers and γβ+1 is a positive integer with even parity. Additionally,
they satisfy the constraint: Pβ+1
b=1 γb = r.
4
Diffusivity for constant and for arbitrary time-dependent
linear and angular speed
In this section we want to illustrate the power of the formulas derived in the previous section,
which hold for any time-dependent linear and angular speed, to compute the statistics in several
specific cases. In particular, we compare the results obtained by using our formulas with the
results that it is possible to obtain by running Monte Carlo simulations. We generate many trials
of random motion via Monte Carlo simulations according to the motion model in (3). The initial
configuration is the same for all the trials and it is x(0) = y(0) = θ(0) = 0. The final curve
length is also the same and it is s = 1. We consider several values for the two parameters Kθ and
Kr and we set µ(s) constant in 4.1 and variable in 4.2. We will see that, as the number of trials
increases, the statistical properties of the motion directly obtained from the trials approach the
ones obtained by using our formulas. In particular, we will refer to the motion diffusivity and we
compute the moments up to the forth order. We start by computing for each trial, the square of
the final distance from the origin, i.e. the distance at s = 1:
D(1)2 = x(1)2 + y(1)2
(19)
Once we have the previous quantity for many trials, we compute its mean value and its variance.
On the other hand, our formulas allow us to directly obtain both the mean value and the
variance. Specifically, equation (15) provides the mean value. Regarding the variance we have:
σD(s)2 ≡

D(s)4
−

D(s)22. Additionally, D(s)4 = x(s)4 + 2x(s)2y(s)2 + y(s)4 = u(s)2w(s)2,
where we used x(s) = u(s)+w(s)
2
and y(s) = u(s)−w(s)
2i
. The quantity

u(s)2w(s)2
can be easily
computed by using equation (17) with p = q = 2 (see appendix C for the details). Its analytical
expression includes the following six terms: n = 0, l = 0, m = 0 (C000); n = 1, l = 0, m = 0
(C100); n = 1, l = 1, m = 1 (C111); n = 1, l = 2, m = 0 (C120); n = 2, l = 2, m = 0 (C220) and
n = 2, l = 2, m = 2 (C222). We have:

D(s)4
= C000 + C100 + C111 + C120 + C220 + C222
(20)
See appendix C for the analytical expression of the previous terms. In the next two subsections,
we compare the results for the mean value and the variance of D(1)2 obtained by running many
trials and by using our analytic expressions. In 4.1 we refer to a constant ratio between the linear
and the angular speed while in 4.2 we refer to a generic varying ratio. Note that this analysis is
a test for the 4th order statistics of the considered Brownian motion.
4.1
Constant ratio between linear and angular speed
When µ(s) = µ0 we have that θ(s) is linear in s:
θ(s) = θ0 + µ0s
(21)
and, the corresponding active trajectory, is a circumference with radius
1
µ0 (note that when the
angular speed vanishes, the radius becomes infinite and the active trajectory becomes a straight
line). The computation of the integrals in (17) and (18) is immediate. As a result, the expression
9
of

D(s)2
can be easily provided in terms of µ0, s, Kr and Kθ. Let us introduce the following
complex quantity:
z ≡−Kθ
2 + iµ0
(22)
By a direct analytical computation of the double integral in (15) we easily obtain:

D(s)2
= Krs + 2ℜ
ezs −1 −zs
z2

(23)
where the symbol ℜ{·} denotes the real part of a given complex quantity.
We do not provide here the analytical expression for

D(s)4
. We only remark that all the
integrals appearing in the expressions in appendix C can be computed in a similar manner and
are expressed as the real part of suitable analytic functions of the following complex quantities:
z, −3Kθ
2
+ iµ0 and −2Kθ + 2iµ0.
Figure 1: Trajectories generated in the case of constant ratio between linear and angular speed.
The blue line is the trajectory generated without noise (i.e, Kθ = Kr = 0) and the other lines
are five trajectories obtained by setting Kθ = Kr = 0.01m.
In the following, we illustrate some numerical results obtained by setting µ(s) = µ0 = 5.
We consider two cases of Brownian noises, the former is characterized by Kθ = Kr = 0.01m
and the latter by Kθ = Kr = 1m. Figure 1 displays the trajectory generated without noise
(blue line) together with five trajectories obtained by setting Kθ = Kr = 0.01m.
For the
trajectory without noise, the final D(1)2 is equal to 0.0573 m2. We use the expression in (23)
to compute the mean value and a similar expression, which depends on the three mentioned
complex quantities z, −3Kθ
2
+ iµ0 and −2Kθ + 2iµ0, to compute the variance. We obtain the
10
Kθ = Kr (m)
trials

D(1)2
(m2)
σ2
D(1)2 (m4)
0.01
103
0.0688
0.0011
0.01
104
0.0664
0.0012
0.01
105
0.0679
0.0012
1
103
1.1189
1.3736
1
104
1.1124
1.2767
1
105
1.1126
1.3035
Table 1: Constant ratio between linear and angular speed. Values of

D(1)2
and σ2
D(1)2 obtained
by running 103, 104 and 105 trials.
following values:

D(1)2
= 0.0680 m2 and σ2
D(1)2 = 0.0012 m4 when Kθ = Kr = 0.01m and

D(1)2
= 1.1130 m2 and σ2
D(1)2 = 1.3052 m4 when Kθ = Kr = 1m.
We ran 105 trials and we computed from them the mean value and the variance. Figure 2
displays the values of D(1)2 obtained in the first 104 trials together with the value of D(1)2 for
the trajectory without noise (red line), the value of

D(1)2
obtained by averaging D(1)2 on
the 104 trials (green line) and the value of

D(1)2
obtained by our formulas (black line). The
image on the left refers to the case Kθ = Kr = 0.01m and the image on the right to the case
Kθ = Kr = 1m.
Figure 2: Constant ratio between linear and angular speed. Values of D(1)2 for 104 trials together
with the value of D(1)2 for the trajectory without noise (red line), the value of

D(1)2
obtained
by averaging D(1)2 on the trials (green line) and the value of

D(1)2
obtained by our formulas
(black line). Kθ = Kr = 0.01m and Kθ = Kr = 1m in the left and right image, respectively.
Table 1 reports the values of

D(1)2
and σ2
D(1)2 obtained by running 103, 104 and 105 trials.
It is possible to see that, as the number of trials increases, the results converge to the values
obtained by using our formulas.
11
Figure 3: Trajectories generated in the case of varying ratio between linear and angular speed.
The blue line is the trajectory generated without noise (i.e, Kθ = Kr = 0) and the other lines
are five trajectories obtained by setting Kθ = Kr = 0.01m.
12
4.2
Arbitrary time-dependent ratio between linear and angular speed
We consider now a generic function µ(s). Specifically, we consider several kinds of dependence
on s obtaining very similar results. In the following, we illustrate the results obtained by setting
µ(s) = 10 s. As in the previous case, we consider the two cases of Brownian noises, i.e., charac-
terized by Kθ = Kr = 0.01m and by Kθ = Kr = 1m. Figure 3 displays the trajectory generated
without noise (blue line) together with five trajectories obtained by setting Kθ = Kr = 0.01m.
For the trajectory without noise, the final D(1)2 is equal to 0.1021 m2. In this case we directly
use the expressions in (15) in (20) and in appendix C by numerically computing all the inte-
grals. We obtain the following values:

D(1)2
= 0.1124 m2 and σ2
D(1)2 = 0.0026 m4 when
Kθ = Kr = 0.01m and

D(1)2
= 1.1443 m2 and σ2
D(1)2 = 1.4681 m4 when Kθ = Kr = 1m.
We ran 105 trials and we computed from them the mean value and the variance. Figure 4
displays the values of D(1)2 obtained in the first 104 trials together with the value of D(1)2 for
the trajectory without noise (red line), the value of

D(1)2
obtained by averaging D(1)2 on
the 104 trials (green line) and the value of

D(1)2
obtained by our formulas (black line). The
image on the left refers to the case Kθ = Kr = 0.01m and the image on the right to the case
Kθ = Kr = 1m.
Figure 4: Varying ratio between linear and angular speed. Values of D(1)2 for 104 trials together
with the value of D(1)2 for the trajectory without noise (red line), the value of

D(1)2
obtained
by averaging D(1)2 on the trials (green line) and the value of

D(1)2
obtained by our formulas
(black line). Kθ = Kr = 0.01m and Kθ = Kr = 1m in the left and right image, respectively.
Table 2 reports the values of

D(1)2
and σ2
D(1)2 obtained by running 103, 104 and 105 trials. As
in the case of constant ratio between linear and angular speed, we remark that, as the number
of trials increases, the results converge to the values obtained by using our formulas.
13
Kθ = Kr (m)
trials

D(1)2
(m2)
σ2
D(1)2 (m4)
0.01
103
0.1139
0.0022
0.01
104
0.1107
0.0026
0.01
105
0.1123
0.0026
1
103
1.1094
1.2534
1
104
1.1605
1.5289
1
105
1.1494
1.4708
Table 2: Varying ratio between linear and angular speed. Values of

D(1)2
and σ2
D(1)2 obtained
by running 103, 104 and 105 trials.
5
Conclusion
In this paper we derived the statistics, up to any order, for 2D−Brownian unicycle dynamics.
The chosen kinematic constraint is modelled by the unicycle differential equation which is a very
general constraint for a 2D motion. According to this model, the mobile robot can freely rotates
but only one direction for the shift is allowed. This model is suitable to characterize the dynamics
of many wheeled robots: namely, the ones that satisfy the unicycle dynamics. Additionally, the
considered deterministic motion is very general. The expressions here provided for the statistics
hold for any deterministic trajectory that satisfies the mentioned kinematic constraint.
The
expressions contain multiple integrals over the deterministic trajectory.
We showed the power of the derived analytical results by investigating the diffusivity of the
considered Brownian motion for constant and arbitrary time-dependent linear and angular speed.
We want to remark that this paper provides the analytical expressions for the statistics, up
to any order, of a non-trivial Brownian motion, i.e. propelled by arbitrary force and torque. To
the best of our knowledge, analytical solutions for any order moment in the case of arbitrary
linear and angular speed have never been derived in the past.
14
A
Computation of ⟨u(s)p w(s)q⟩
We have:
⟨u(s)p w(s)q⟩= lim
N→∞
X
j1···jpk1···kq
D
ei(θj1+···+θjp−θk1−···−θkq )E
(24)

(δs + ϵj1) · · · (δs + ϵjp)(δs + ϵk1) · · · (δs + ϵkq)

where each index goes from 1 to N. Let us focus our attention on the quantity

(δs + ϵj1) · · · (δs + ϵjp)(δs + ϵk1) · · · (δs + ϵkq)

We can write this quantity as the sum of p+q+1 terms, i.e., δsp+qC0+δsp+q−1C1+δsp+q−2C2+
· · · + δsCp+q−1 + Cp+q. We trivially have C0 = 1. Concerning the remaining coefficients, we
remark, first of all, that only the ones with even index are different from 0. In other words,
C2n+1 = 0, n = 0, 1, · · · , ⌊p+q
2 ⌋(where ⌊x⌋is the largest integer not greater than x). Let us
compute C2n. This coefficient is the sum of
p+q
2n


elements, each of them being the average
of a product of 2n terms "ϵ". On the other hand, since the exponential in (24) is symmetric
with respect to the change ja ↔ja′, a, a′ = 1, · · · , p and with respect to the change kb ↔kb′,
b, b′ = 1, · · · , q but it is not symmetric with respect to the change ja ↔kb, we need to separate
the
p+q
2n


elements in several groups. Specifically, let us denote by l the number of ϵ with index
of type k. The number of elements belonging to this group is
p
2n−l

q
l


. Note that we set
a
b


= 0
when b > a.
Note also that P2n
l=0
p
2n−l

q
l


=
p+q
2n


.
We now remark that, because of the
statistical properties of ϵ, the average of the product of 2n terms ϵ is different from zero only
when their indexes are equal two-by-two. Hence, we have to consider all the combinations of
products of terms ϵ, whose indexes are equal two-by-two and which differ for at least one pair
of indexes. On the other hand, for a given element in the group characterized by a given n and
l, the effect in (24) depends on the number of pairs which are hetero (i.e., with an index of type
j and one of type k). Let us denote by m the number of pairs which are hetero. Note that m
has the same parity of l. Indeed, by definition we have l indexes of type k. Additionally, we are
using m indexes of type j and m indexes of type k to make m hetero pairs. Hence, l −m indexes
of type k remain and, with them, we have to make l−m
2
homo pairs of type k. Similarly, we have
to make 2n−l−m
2
homo pairs of type j. Hence, l−m
2
must be integer.
We compute the number of combinations of the products of 2n terms ϵ, with l indexes of
type k, where the indexes are equal two-by-two, with m pairs which are hetero and which differ
for at least one pair. We obtain:
2n−l
m

 l
m


m!(2n −l −m −1)!!(l −m −1)!!. Following this, we
can write (24) as follows:
⟨u(s)p w(s)q⟩= lim
N→∞
X
{j}p{k}q
⌊p+q
2
⌋
X
n=0
δsp+q−2n
(25)
2n
X
l=0

p
2n −l
q
l

min(l,2n−l)
X
m=0 odd/even
2n −l
m
 l
m

m!
(2n −l −m −1)!!(l −m −1)!!

ϵ2
j1

δj1j2

ϵ2
j3

δj3j4 · · ·
D
ϵ2
j2n−l−m−1
E
δj2n−l−m−1j2n−l−m

ϵ2
k1

δk1k2

ϵ2
k3

δk3k4 · · ·
D
ϵ2
kl−m−1
E
δkl−m−1kl−m
D
ϵ2
kl−m+1
E
δkl−m+1j2n−l−m+1 · · ·
15

ϵ2
kl

δklj2n−l
D
ei(θj1+···+θjp−θk1−···−θkq )E
where, for the brevity sake, we denoted by P
{j}p{k}q the sum PN
j1···jpk1···kq=1. Each average

ϵ2
provides Krδs. Hence, all together, they provide Kn
r δsn. Additionally, the number of Kronecker
deltas is n. Hence, in the limit of N →∞for each value of n we get a multiple integral of
dimension p + q −n. Note that, when the indexes are equal h-by-h (h > 2), the result in the
limit N →∞vanishes since, the power of δs, is larger than the number of sums. By a direct
computation in (25) we obtain:
⟨u(s)p w(s)q⟩=
⌊p+q
2
⌋
X
n=0
Kn
r
2n
X
l=0

p
2n −l
q
l

min(l,2n−l)
X
m=0 odd/even
(26)
2n −l
m
 l
m

m!(2n −l −m −1)!!(l −m −1)!! sm lim
N→∞
X
{js}σ{jd}ρ{ks}χ{kd}η
δsβ D
exp
n
i[θjs
1 + · · · + θjsσ + 2(θjd
1 +
+ · · · + θjd
ρ) −(θks
1 + · · · + θksχ) −2(θkd
1 + · · · + θkdη)]
oE
where ρ ≡2n−l−m
2
, σ ≡p −(2n −l), η ≡l−m
2 , χ ≡q −l and β ≡ρ + σ + η + χ = p + q −n −m.
We must compute the average of the exponential in (26) and then we compute the limit N →∞.
We remark that the various θ in the exponential contain random quantities (i.e., the δθ at
different time steps). In order to proceed we have to separate all the random quantities which are
independent. We start this separation by redefining the indexes in the sum P
{js}σ{jd}ρ{ks}χ{kd}η.
Specifically, we consider the new indexes j
sj
dk
sk
d which differ from jsjdkskd since they are
ordered (in increasing order). For instance, j
s
1 < j
s
2 < · · · < j
s
σ. Hence, the last sum in (26)
can be replaced with P
{js}σ{jd}ρ{ks}χ{kd}η →ρ!σ!η!χ! P
{js}σ{jd}ρ{k
s}χ{k
d}η. The four types of
indexes are not ordered among them. Hence, the sum includes all the possible combinations which
maintain the order only restricted to a single index type. For instance, a possible combination
is: k
d
1 < k
d
2 < j
s
1 < k
s
1 < j
s
2 < k
d
3 < j
d
1 < j
d
2 < · · · . We introduce the β ordered indexes
w1 < w2 < · · · < wβ. Additionally, let us denote with ∆b
a ≡Pb
c=a+1 δθc and ∆
b
a ≡Pb
c=a+1 δθc =
θ(bδs)−θ(aδs), where θ(s) is the deterministic trajectory. Finally, let us define α ≡2ρ+σ−χ−2η.
The sum in the exponential in (26) contains α(∆w1
0 +∆
w1
0 ). Then, depending on the combination
of the indexes j
sj
dk
sk
d, we have a different result for the (∆w2
w1+∆
w2
w1). Specifically, if w1 = j
d
1 the
sum in the exponential contains (α −2)(∆w2
w1 +∆
w2
w1). If w1 = j
s
1 it contains (α −1)(∆w2
w1 +∆
w2
w1).
If w1 = k
s
1 it contains (α + 1)(∆w2
w1 + ∆
w2
w1). If w1 = k
d
1 it contains (α + 2)(∆w2
w1 + ∆
w2
w1). We
introduce the vector c ≡[c1, · · · , cβ]T whose entries are −2, −1, 1 and 2. Specifically, it contains
η entries equal to 2, χ equal to 1, σ equal to −1 and ρ equal to −2. It is easy to realize that we
have
β
χ

β−χ
σ

ρ+η
η


vectors c. Finally, we define Φb ≡Pb
a=1 ca. Note that Φβ = −α. According
to this, the last sum in (26) can be written as follows:
eiαθ0ρ!σ!η!χ!
X
c
X
{w}β
δsβ·
·
*
exp
(
i
β−1
X
b=0
h
(α + Φb)(∆wb+1
wb
+ ∆
wb+1
wb
)
i)+
16
Hence, we have:
⟨u(s)p w(s)q⟩=
⌊p+q
2
⌋
X
n=0
Kn
r
2n
X
l=0

p
2n −l
q
l

·
(27)
·
min(l,2n−l)
X
m=0 odd/even
2n −l
m
 l
m

m!(2n −l −m −1)!!(l −m −1)!!·
·smeiαθ0ρ!σ!η!χ!
X
c
lim
N→∞
X
{w}β
δsβ·
·
*
exp
(
i
β−1
X
b=0
h
(α + Φb)(∆wb+1
wb
+ ∆
wb+1
wb
)
i)+
Now we can compute the average since we were able to separate all the independent quantities.
We have:
*
exp
(
i
β−1
X
b=0
h
(α + Φb)(∆wb+1
wb
+ ∆
wb+1
wb
)
i)+
=
(28)
= exp
(
i
β−1
X
b=0
h
(α + Φb)∆
wb+1
wb
i)
·
·
β−1
Y
b=0
D
exp
n
i
h
(α + Φb)∆wb+1
wb
ioE
=
= exp
(
i
β−1
X
b=0
h
(α + Φb)∆
wb+1
wb
i)
·
·
β−1
Y
b=0
exp

−(α + Φb)2(wb+1 −wb)Kθδs
2

=
e
Pβ−1
b=0

i(α+Φb)∆
wb+1
wb
−
(α+Φb)2(wb+1−wb)Kθδs
2

We used the equality (η = N (0, 1)):

eAη
= e
A2
2
(29)
with A = i(α+Φb)σb+1 and σ2
b+1 = (wb+1 −wb)Kθδs. By substituting equation (28) in (27) and
by taking the limit (N →∞), we finally obtain:
⟨u(s)p w(s)q⟩=
⌊p+q
2
⌋
X
n=0
Kn
r
2n
X
l=0

p
2n −l
q
l

·
(30)
·
min(l,2n−l)
X
m=0 odd/even
2n −l
m
 l
m

m!(2n −l −m −1)!!(l −m −1)!!·
17
·smeiαθ0ρ!σ!η!χ!
X
c
Z s
0
ds1
Z s
s1
ds2 · · ·
Z s
sβ−1
dsβ
e
Pβ−1
b=0

i(α+Φb)[θ(sb+1)−θ(sb)]−
(α+Φb)2(sb+1−sb)Kθ
2

which coincides with (17).
B
Computation of
D
u(s)p w(s)q ˜θrE
This computation follows the same initial steps carried out in appendix A. We obtain an ex-
pression equal to the one given in (27) but the mean value at the end must be replaced with
D
exp
n
i Pβ−1
b=0
h
(α + Φb)(∆wb+1
wb
+ ∆
wb+1
wb
)
io
· ˜θrE
. The deterministic part can be factorized out
the mean value. We need to calculate:
*
exp
(
i
β−1
X
b=0
h
(α + Φb)∆wb+1
wb
i)
· ˜θr
+
=
=
*
˜θr
β−1
Y
b=0
n
exp
h
i(α + Φb)∆wb+1
wb
io+
By using the multinomial theorem we can write:
˜θr =
 β
X
b=0
∆wb+1
wb
!r
=
X
γ
r! ·
β
Y
b=0

∆wb+1
wb
γb+1
γb+1!
where wβ+1 ≡N and the sum over γ is the sum over all the vectors γ = [γ1, · · · , γβ+1], where γb
(b = 1, · · · , β + 1) are positive integers satisfying the constraint: Pβ+1
b=1 γb = r. Hence, we obtain
*
˜θr
β−1
Y
b=0

ei(α+Φb)∆
wb+1
wb
+
=
(31)
X
γ
r!
*
∆wβ+1
wβ
γβ+1
γβ+1!
β−1
Y
b=0
ei(α+Φb)∆
wb+1
wb

∆wb+1
wb
γb+1
γb+1!
+
=
=
X
γ
r!
Qβ
b=0 γb+1!
D
∆wβ+1
wβ
γβ+1E
·
·
β−1
Y
b=0

ei(α+Φb)∆
wb+1
wb

∆wb+1
wb
γb+1
We use the following two standard results for a normal distribution:
D
∆wβ+1
wβ
γβ+1E
=

0
γβ+1
odd
σγβ+1
β+1 (γβ+1 −1)!!
γβ+1
even
(32)
and
D
exp
h
i(α + Φb)∆wb+1
wb
i 
∆wb+1
wb
γb+1E
=
(33)
18
= σγb+1
b+1 exp

−(α + Φb)2σ2
b+1
2

·
·
⌊
γb+1
2
⌋
X
a=0
γb+1! (i(α + Φb)σb+1)γb+1−2a
a!(γb+1 −2a)!2a
with σ2
b+1 ≡Kθδs(wb+1−wb). Equation (33) is obtained starting from (29) and by differentiating
γb+1 times with respect to A. Hence, the expression in (31), becomes:
*
˜θr
β−1
Y
b=0
n
exp
h
i(α + Φb)∆wb+1
wb
io+
=
(34)
=
X
γ
r!σγβ+1
β+1 (γβ+1 −1)!!
Qβ
b=0 γb+1!
β−1
Y
b=0

σγb+1
b+1 e−
(α+Φb)2σ2
b+1
2
·
·
⌊
γb+1
2
⌋
X
a=0
γb+1! (i(α + Φb)σb+1)γb+1−2a
a!(γb+1 −2a)!2a



where now the sum over γ only includes the vectors γ whose last entry is even. By taking the
limit N →∞, this expression becomes:
X
γ
r!(γβ+1 −1)!!
Qβ
b=0 γb+1!
(s −sβ)
γβ+1
2
K
r
2
θ ·
·
β−1
Y
b=0

(sb+1 −sb)
γb+1
2 e−
(α+Φb)2Kθ(sb+1−sb)
2
·
·
⌊
γb+1
2
⌋
X
a=0
γb+1!

i(α + Φb)
p
Kθ(sb+1 −sb)
γb+1−2a
a!(γb+1 −2a)!2a





and we finally obtain:
D
u(s)p w(s)q ˜θ(s)rE
=
⌊p+q
2
⌋
X
n=0
Kn
r
2n
X
l=0

p
2n −l
q
l

·
(35)
·
min(l,2n−l)
X
m=0 odd/even
2n −l
m
 l
m

m!(2n −l −m −1)!!(l −m −1)!!
·smeiαθ0ρ!σ!η!χ!
X
c
X
γ
r!(γβ+1 −1)!!K
r
2
θ
Qβ
b=0 γb+1!
Z s
0
ds1
Z s
s1
ds2·
· · ·
Z s
sβ−1
dsβ(s −sβ)
γβ+1
2
ei Pβ−1
b=0 (α+Φb)[θ(sb+1)−θ(sb)]
β−1
Y
b=0

(sb+1 −sb)
γb+1
2 e−
(α+Φb)2Kθ(sb+1−sb)
2
·
19
⌊
γb+1
2
⌋
X
a=0
γb+1!

i(α + Φb)
p
Kθ(sb+1 −sb)
γb+1−2a
a!(γb+1 −2a)!2a





which coincides with (18).
C
Computation of

D(s)4
=

u(s)2 w(s)2
We use equation (17) with p = q = 2 to obtain the analytical expression of

D(s)4
=

u(s)2w(s)2
.
The first sum on n includes the three terms n = 0, 1, 2.
When n = 0, the
second sum on l and the third sum on m only include one term, i.e., l = m = 0. All the factors
before the sum on c are 1 with the exception of (p −2n + l)! = 2 and (q −l)! = 2. Additionally,
β = 4 and the sum on c consists of six vectors, i.e., [−1, −1, 1, 1], [−1, 1, −1, 1], [−1, 1, 1, −1],
[1, −1, −1, 1], [1, −1, 1, −1], [1, 1, −1, −1]. Hence we obtain the following term:
C000 = 8
Z s
0
ds1
Z s
s1
ds2
Z s
s2
ds3
Z s
s3
ds4
n
e−
Kθ(−s1−3s2+3s3+s4)
2
cos
−θ1 −θ2 + θ3 + θ4


+
+e−
Kθ(−s1+s2−s3+s4)
2

cos
−θ1 + θ2 −θ3 + θ4


+ cos
−θ1 + θ2 + θ3 −θ4

o
where θj ≡θ(sj), j = 1, 2, 3, 4.
Let us consider now the terms with n = 1. In this case, l = 0, 1, 2. For each value of l we
have in general several values of m, which must have the same parity of l. Hence, when l = 1,
we only have m = 1. When l = 2, we can have m = 0, 2. On the other hand, when l = m = 2,
the factor
2n−l
m


=
0
2


vanishes. By an explicit computation, it is possible to see that C120 is
the conjugate of C100. Hence, their sum is real and it is:
C100 + C120 = 4Kr
Z s
0
ds1
Z s
s1
ds2
Z s
s2
ds3
n
e−
Kθ(−4s1+3s2+s3)
2
cos
2θ1 −θ2 −θ3


+
+e−
Kθ(−s1+s3)
2
cos
−θ1 + 2θ2 −θ3


+ e−
Kθ(−s1−3s2+4s3)
2
cos
−θ1 −θ2 + 2θ3

o
The last term with n = 1 is the one with l = m = 1. By a direct computation we obtain:
C111 = 8Krs
Z s
0
ds1
Z s
s1
ds2e−
Kθ(−s1+s2)
2
cos
−θ1 + θ2


Finally, when n = 2 the sum on l includes the five values: l = 0, 1, 2, 3, 4. On the other hand,
only the term l = 2 does not vanish since the product
p
2n−l

q
l


is zero in all the other cases.
When l = 2 we have two possible values for m, i.e., m = 0, 2. We obtain the following two
contributions:
C220 = 2K2
r
Z s
0
ds1
Z s
s1
ds2e−2Kθ(−s1+s2) cos
−2θ1 + 2θ2


C222 = 2K2
rs2
20
References
[1] Sepideh Bazazi, Pawel Romanczuk, Sian Thomas, Lutz Schimansky-Geier, Joseph J. Hale,
Gabriel A. Miller, Gregory A. Sword, Stephen J. Simpson and Iain D. Couzin, Nutritional
state and collective motion: from individuals to mass migration, Proc. R. Soc. B 278, 356
(2010).
[2] Callan, David (2009), A combinatorial survey of identities for the double factorial,
arXiv:0906.1317
[3] S. Challa and Y. Bar-Shalom, "Nonlinear filter design using Fokker-Planck-Kolmogorov prob-
ability density evolutions", IEEE Transactions on Aerospace and Electronic Systems, vol. 36,
no. 1, 2000
[4] I. D. Couzin et al., Collective memory and spatial sorting in animal groups, J. Theor. Biol.
218, 1 (2002).
[5] F. E. Daum, "New nonlinear filters and exact solutions of the Fokker-Planck equation",
Proceedings of IEEE American Control Conference, pp.884 -888 1986
[6] F. E. Daum, "Solution of the Zakai equation by separation of variables", IEEE Transactions
on Automatic Control, vol. AC-32, no. 10, pp.941 -943 1987
[7] F. E. Daum, "Nonlinear filters: beyond the Kalman filter", Aerospace and Electronic Systems
Magazine, IEEE, Volume:20, Issue: 8, 2005
[8] R. Dreyfus, J. Baudry, M. L. Roper, M. Fermigier, H. A. Stone, and J. Bibette, Microscopic
artificial swimmers, Nature London, 437, 862, 2005
[9] Freidlin Mark and Hu Wenqing, On Diffusion in Narrow Random Channels, J Stat Phys,
152:136–158, 2013
[10] Hanggi Peter, Marchesoni Fabio, Artificial Brownian motors: Controlling transport on the
nanoscale, Reviews of Modern Physics 81 (1): 387, 2009
[11] J. C. Latombe, Robot motion planning, Kluwer Academic Publishers, 1991
[12] E. Lauga and T. R. Powers, The hydrodynamics of swimming microorganisms, Rep. Prog.
Phys. 72, 096601, 2009
[13] M. Leoni, J. Kotar, B. Bassetti, P. Cicuta, and M. C. Lagomarsino, A basic swimmer at low
Reynolds number, Soft Matter 5, 472, 2009
[14] Long, A.W., Wolfe, K.C., Mashner, M.J., Chirikjian, G.S., “The Banana Distribution is
Gaussian: A Localization Study with Exponential Coordinates,” Robotics: Science and Sys-
tems, Sydney, NSW, Australia, July 09 - July 13, 2012.
[15] A. Martinelli, The odometry error of a mobile robot with a synchronous drive system, IEEE
Transactions on Robotics and Automation, Volume 18, Issue 3, 399–405, Jun 2002.
[16] A. Martinelli, The accuracy on the parameter estimation of an odometry system of a mobile
robot, ICRA ’02. IEEE International Conference on Robotics and Automation, 1378–1383
vol.2, 2002.
[17] H. Niwa, Self-organizing dynamic model offish schooling, J. Theor. Biol. 171, 123 (1994).
21
[18] Bernt Oksendal, Stochastic Differential Equations:
An Introduction with Applications,
Springer, Berlin, 2003
[19] Papoulis, A. Probability, Random Variables, and Stochastic Processes, 2nd ed. New York:
McGraw-Hill, 1984.
[20] Park, W., Kim, J.S., Zhou, Y., Cowan, N.J., Okamura, A.M., Chirikjian, G.S., “Diffusion-
based motion planning for a nonholonomic flexible needle model,” ICRA, Barcelona, Spain,
April 2005
[21] Park, W., Liu, Y., Zhou, Y., Moses, M., Chirikjian, G.S., “Kinematic State Estimation
and Motion Planning for Stochastic Nonholonomic Systems Using the Exponential Map,”
Robotica 26(4), 419-434. July-August 2008
[22] Park, W., Reed, K. B., Okamura, A. M., Chirikjian, G. S., “Estimation of model parameters
for steerable needles,” ICRA, Anchorage, Alaska, 2010.
[23] Reimann Peter, Brownian motors: noisy transport far from equilibrium, Physics Reports,
361, 57–265, 2002
[24] Hannes Risken, The Fokker-Planck Equation Methods of Solution and Applications,
Springer Series in Synergetics Volume 18, 1989
[25] A. Sokolov, I. S. Aranson, J. O. Kessler, and R. E. Goldstein, Concentration Dependence of
the Collective Dynamics of Swimming Bacteria, Phys. Rev. Lett. 98, 158102
[26] P. Tierno, R. Golestanian, I. Pagonabarraga, and F. Sagues, Magnetically actuated colloidal
microswimmers, J. Phys. Chem. B. 112, 16525 2008
[27] P Tierno, R Golestanian, I Pagonabarraga, F Sagués, Controlled swimming in confined
fluids of magnetically actuated colloidal rotors, Phys Rev Lett 101 218304, 2008
[28] R. Trouilloud, T. S. Yu, A. E. Hosoi and E. Lauga, Soft swimming: Exploiting deformable
interfaces for low-Reynolds number locomotion, Phys. Rev. Lett. 101, 048102, 2008
[29] Turner Jonathan, Brownian Motion Applied to Human Intersections, thesis at Texas State
University-San Marcos, Dept. of Mathematics, 2012
[30] Zhou, Y., Chirikjian, G.S., “Probabilistic Models of Dead-Reckoning Error in Nonholonomic
Mobile Robots,” ICRA’03, Taipei, Taiwan, September, 2003.
[31] Zhou Y., Chirikjian G.S., Planning for Noise-Induced Trajectory Bias in Nonholonomic
Robots with Uncertainty, 2004 IEEE International Conference on Robotics and Automation
(ICRA 2004), pp.4596-4601, New Orleans, Louisiana, 2004.
22