Agostino Martinelli 
  
Complete analytic solution to Brownian unicycle dynamics

Keywords: Brownian motion, Stochastic dynamics, Stochastic processes, Localization

This paper derives a complete analytical solution for the probability distribution of the configuration of a non-holonomic mobile robot that moves in two spatial dimensions by satisfying the unicycle kinematic constraints. The proposed solution differs from previous solutions since it 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 deterministic 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.

Résumé : Ce document donne une solution analytique complète pour la distribution de probabilité de la configuration d'un robot mobile non-holonome, qui se déplace dans deux dimensions spatiales en satisfaisant les contraintes cinématiques du monocycle. La solution proposée est différente par rapport aux solutions précédentes car elle est obtenue par dérivation de l'expression analytique de tous les moments de la distribution de probabilité. Au meilleur de notre connaissance, une expression analytique pour les moments de tout les ordres qui est valide même dans le cas d'une vitesse linéaire arbitraire et d'une vitesse angulaire arbitraire, n'a jamais été dérivé auparavant. Pour calculer ces moments, une intégration directe de l'équation de Langevin est réalisée et chaque moment de probabilité est exprimé en intégrale multiple du mouvement déterministe (c'est á dire, du mouvement connu qui se traduirait par l'absence de bruit). Pour le cas particulier où le rapport entre la vitesse linéaire et angulaire est constante, les intégrales multiples peuvent être facilement résolus et exprimés comme la partie réelle et la partie imaginaire de fonctions analytiques appropriés. Comme application des résultats analytiques dérivés, on examine la diffusivité du mouvement brownien pour un rapport entre la vitesse linéaire et la vitesse angulaire constante et arbitraire.

Mots-clés : Mouvement brownien, Dynamique stochastique, Processus stochastiques, Localisation 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 [START_REF] Mark | On Diffusion in Narrow Random Channels[END_REF][START_REF] Peter | Artificial Brownian motors: Controlling transport on the nanoscale[END_REF][START_REF] Peter | Brownian motors: noisy transport far from equilibrium[END_REF], biological and artificial microswimmers [START_REF] Dreyfus | Microscopic artificial swimmers[END_REF][START_REF] Lauga | The hydrodynamics of swimming microorganisms[END_REF][START_REF] Leoni | A basic swimmer at low Reynolds number[END_REF][START_REF] Sokolov | Concentration Dependence of the Collective Dynamics of Swimming Bacteria[END_REF][START_REF] Tierno | Magnetically actuated colloidal microswimmers[END_REF][START_REF] Tierno | Controlled swimming in confined fluids of magnetically actuated colloidal rotors[END_REF][START_REF] Trouilloud | Soft swimming: Exploiting deformable interfaces for low-Reynolds number locomotion[END_REF], macroscopic animals [START_REF] Bazazi | Nutritional state and collective motion: from individuals to mass migration[END_REF][START_REF] Couzin | Collective memory and spatial sorting in animal groups[END_REF][START_REF] Niwa | Self-organizing dynamic model offish schooling[END_REF] and even pedestrians [START_REF] Turner | Brownian Motion Applied to Human Intersections[END_REF].

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 [START_REF] Zhou | Probabilistic Models of Dead-Reckoning Error in Nonholonomic Mobile Robots[END_REF][START_REF] Zhou | Planning for Noise-Induced Trajectory Bias in Nonholonomic Robots with Uncertainty[END_REF], and the 3D case of flexible needle steering [START_REF] Park | Diffusion-based motion planning for a nonholonomic flexible needle model[END_REF][START_REF] Park | Estimation of model parameters for steerable needles[END_REF]. In these works, analytical solutions are obtained by directly solving the corresponding Fokker-Planck equations by using the Fourier transform. Another method proposed for obtaining a closed-form solution is the use of exponential coordinates [START_REF] Long | The Banana Distribution is Gaussian: A Localization Study with Exponential Coordinates[END_REF][START_REF] Park | Kinematic State Estimation and Motion Planning for Stochastic Nonholonomic Systems Using the Exponential Map[END_REF]. 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 [START_REF] Latombe | Robot motion planning[END_REF]. Our stochastic motion model (presented in section 2) differs from the one considered in [START_REF] Long | The Banana Distribution is Gaussian: A Localization Study with Exponential Coordinates[END_REF][START_REF] Park | Kinematic State Estimation and Motion Planning for Stochastic Nonholonomic Systems Using the Exponential Map[END_REF][START_REF] Zhou | Probabilistic Models of Dead-Reckoning Error in Nonholonomic Mobile Robots[END_REF][START_REF] Zhou | Planning for Noise-Induced Trajectory Bias in Nonholonomic Robots with Uncertainty[END_REF] since we assume independent errors acting on the linear and angular components. The model adopted in [START_REF] Long | The Banana Distribution is Gaussian: A Localization Study with Exponential Coordinates[END_REF][START_REF] Park | Kinematic State Estimation and Motion Planning for Stochastic Nonholonomic Systems Using the Exponential Map[END_REF][START_REF] Zhou | Probabilistic Models of Dead-Reckoning Error in Nonholonomic Mobile Robots[END_REF][START_REF] Zhou | Planning for Noise-Induced Trajectory Bias in Nonholonomic Robots with Uncertainty[END_REF] 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 [START_REF] Challa | Nonlinear filter design using Fokker-Planck-Kolmogorov probability density evolutions[END_REF][START_REF] Daum | New nonlinear filters and exact solutions of the Fokker-Planck equation[END_REF][START_REF] Daum | Solution of the Zakai equation by separation of variables[END_REF][START_REF] Daum | Nonlinear filters: beyond the Kalman filter[END_REF] 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 [START_REF] Martinelli | The odometry error of a mobile robot with a synchronous drive system[END_REF][START_REF] Martinelli | The accuracy on the parameter estimation of an odometry system of a mobile robot[END_REF] (in the expressions in [START_REF] Martinelli | The odometry error of a mobile robot with a synchronous drive system[END_REF] there are some typos, corrected in [START_REF] Martinelli | The accuracy on the parameter estimation of an odometry system of a mobile robot[END_REF]). Here, we extend the computation in [START_REF] Martinelli | The odometry error of a mobile robot with a synchronous drive system[END_REF] 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 [START_REF] Martinelli | The odometry error of a mobile robot with a synchronous drive system[END_REF]). As a result, this paper provides a complete analytical solution to a Fokker-Planck equation for which the detailed balance condition is not verified 1 .

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.

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 [START_REF] Latombe | Robot motion planning[END_REF], 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 ):

  ẋ = v cos θ ẏ = 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)dw r (t)] dy(t) = sin θ(t) [v(t)dt + α(t)dw r (t)] dθ(t) = ω(t)dt + β(t)dw θ (t) (2) 
where [dw r (t), dw θ (t)] T is a standard Wiener process of dimension two [START_REF] Oksendal | Stochastic Differential Equations: An Introduction with Applications[END_REF]. 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) = K r |v(t)|, with K r 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) = 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 [dw r (t), dw θ (t)] T . Hence, we 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 2 r dw r (s)] dy(s) = sin θ(s) [ds + K 1 2 r dw r (s)] dθ(s) = µ(s)ds + K 1 2 θ dw θ (s) (3) 
The associated Smoluchowski equation is:

∂p ∂s = -∇ • (D 1 p) + ∇ • (D 2 ∇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;

• D 1 = [cos θ, sin θ, µ] T is the drift vector;

• D 2 = 1 2   K r cos 2 θ K r cos θ sin θ 0 K r cos θ sin θ K r sin 2 θ 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.

Computation of the probability distribution

The probability density p(x, y, θ; s) satisfies the Smoluchowski partial differential equation in [START_REF] Couzin | Collective memory and spatial sorting in animal groups[END_REF], which is a special case of the Fokker-Planck equation. Since the detailed balance condition is not satisfied [START_REF] Risken | The Fokker-Planck Equation Methods of Solution and Applications[END_REF], 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) ≡ dx dy p(x, y, θ; s), is a Gaussian distribution with mean value θ 0 + 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 + 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.

First and second-order statistics

Let us consider the first equation in [START_REF] Challa | Nonlinear filter design using Fokker-Planck-Kolmogorov probability density evolutions[END_REF]. We compute the expression of x(s) by a direct integration. We divide the interval (0, s) in N equal segments, δs ≡ s N . We have:

x(s) = lim N →∞ N j=1 (δs + j ) cos θ j ( 6 
)
where j is a random Gaussian variable satisfying j = 0, j k = δ jk K r δ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 m=1 δθ m ≡ θ j + ∆θ j (7) 
where, according to our stochastic model in (3), δθ j is a random Gaussian variable satisfying δθ j = 0, δθ j k = 0 and 6) and [START_REF] Daum | Nonlinear filters: beyond the Kalman filter[END_REF] 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:

δθ j δθ k = δ jk K θ δs for j, k = 1, • • • , N . As a result, ∆θ j ≡ j 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 (
x(s) = lim N →∞ N j=1 (δs + j ) cos θ j = (8) = lim N →∞ N j=1 δs cos θ j e -k θ jδs 2 = s 0 ds cos θ(s ) e -K θ s 2 
Similarly, it is possible to obtain the mean value of y(s):

y(s) = s 0 ds sin θ(s ) e -K θ s 2 (9) 
INRIA Finally, in a similar manner, but with some more computation, we obtain all the second-order moments:

x(s) 2 = s 0 ds s-s 0 ds e -K θ s 2 {(1 + χ c (s )) (10) cos[θ(s + s ) -θ(s )] -χ s (s ) sin[θ(s + s ) -θ(s )] + + K r 2 s + s 0 ds χ c (s ) y(s) 2 = s 0 ds s-s 0 ds e -K θ s 2 {(1 -χ c (s )) (11) cos[θ(s + s ) -θ(s )] + χ s (s ) sin[θ(s + s ) -θ(s )] + + K r 2 s - s 0 ds χ c (s ) x(s) y(s) = s 0 ds s-s 0 ds e -K θ s 2 {χ s (s ) (12) cos[θ(s + s ) -θ(s )] + χ c (s ) sin[θ(s + s ) -θ(s )] + + K r 2 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) K r s + 2 s 0 ds s-s 0 ds e -K θ s 2 cos[θ(s + s ) -θ(s )] RR n°8465

Computation of any-order moment

We introduce the following two complex random quantities:

u(s) ≡ lim N →∞ N j=1 (δs + j )e iθj ; v(s) ≡ lim N →∞ N j=1 (δs + j )e -iθj (16) 
From ( 6) and ( 16) it is immediate to realize that x(s) = u(s)+v(s)

2

. Similarly we have y(s) = u(s)-v(s)

2i

. Hence, in order to compute any-order moment which involves x(s) and y(s) it suffices to compute u(s) p v(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 v(s) q = p+q 2 n=0 K n r 2n l=0 p 2n -l q l ( 17 
) m 2n -l m l m m!(2n -l -m -1)!!(l -m -1)!! s m e i(p-q)θ0 2n -l -m 2 ! l -m 2 !(p -2n + l)!(q -l)! c s 0 ds 1 s s1 ds 2 • • • s s β-1 ds β exp β-1 b=0 [i(p -q + Φ b )• •[θ(s b+1 ) -θ(s b )] - (p -q + Φ b ) 2 (s b+1 -s b )K θ 2
where:

• p+q 2
is the largest integer not greater than p+q 2 ;

• the second sum on l (i.e., 2n 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., 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 [START_REF] Callan | A combinatorial survey of identities for the double factorial[END_REF] (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 s m );

• the sum on c, i.e., 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 c consists of

β q-l β-q+l p-2n+l n-m l-m 2 addends); INRIA • Φ b ≡ b a=1
c a (note that Φ β = q -p); • s 0 ≡ 0 and θ(s 0 ) = θ 0 (note that s 0 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 u(s) p v(s) q θ(s) r , where θ(s) ≡ θ(s) -θ(s). The details of this computation are provided in appendix B. We have, for any p, q, r ∈ N :

u(s) p v(s) q θ(s) r = p+q 2 n=0 K n r 2n l=0 p 2n -l q l ( 18 
) m 2n -l m l m m!(2n -l -m -1)!!(l -m -1)!! s m e i(p-q)θ0 2n -l -m 2 ! l -m 2 !(p -2n + l)!(q -l)! c γ r!(γ β+1 -1)!!K r 2 θ β b=0 γ b+1 ! s 0 ds 1 s s1 ds 2 • • • s s β-1 ds β (s -s β ) γ β+1 2 exp i β-1 b=0 (p -q + Φ b )[θ(s b+1 ) -θ(s b )] • • β-1 b=0 (s b+1 -s b ) γ b+1 2 exp - (p -q + Φ b ) 2 K θ (s b+1 -s b ) 2 γ b+1 2 a=0 γ b+1 ! i(p -q + Φ b ) K θ (s b+1 -s b ) γ b+1 -2a a!(γ b+1 -2a)!2 a     
where the sum over γ, i.e. γ , 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:

β+1 b=1 γ b = r.
4 Diffusivity for constant and for arbitrary timedependent 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 [START_REF] Challa | Nonlinear filter design using Fokker-Planck-Kolmogorov probability density evolutions[END_REF]. 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 K r 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 [START_REF] Martinelli | The odometry error of a mobile robot with a synchronous drive system[END_REF] provides the mean value. Regarding the variance we have:

σ D(s) 2 ≡ D(s) 4 -D(s) 2 2 . Additionally, D(s) 4 = x(s) 4 + 2x(s) 2 y(s) 2 + y(s) 4 = u(s) 2 w(s) 2
, where we used x(s) = u(s)+w(s) 2 and y(s) = u(s)-w(s)

2i

. The quantity u(s) 2 w(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 (C 000 ); n = 1, l = 0, m = 0 (C 100 ); n = 1, l = 1, m = 1 (C 111 ); n = 1, l = 2, m = 0 (C 120 ); n = 2, l = 2, m = 0 (C 220 ) and n = 2, l = 2, m = 2 (C 222 ).
We have:

D(s) 4 = C 000 + C 100 + C 111 + C 120 + C 220 + C 222 (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 4 th order statistics of the considered Brownian motion.

Constant ratio between linear and angular speed

When µ(s) = µ 0 we have that θ(s) is linear in s:

θ(s) = θ 0 + µ 0 s (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 [START_REF] Niwa | Self-organizing dynamic model offish schooling[END_REF] and ( 18) is immediate. As a result, the expression of D(s) 2 can be easily provided in terms of µ 0 , s, K r 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 [START_REF] Martinelli | The odometry error of a mobile robot with a synchronous drive system[END_REF] we easily obtain:

D(s) 2 = K r s + 2 e zs -1 -zs z 2 (23) 
where the symbol {•} denotes the real part of a given complex quantity.

INRIA

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 θ = K r = 0) and the other lines are five trajectories obtained by setting K θ = K r = 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 θ = K r = 0.01m and the latter by K θ = K r = 1m. Figure 1 displays the trajectory generated without noise (blue line) together with five trajectories obtained by setting K θ = K r = 0.01m. For the trajectory without noise, the final D(1) 2 is equal to 0.0573 m 2 . We use the expression in [START_REF] Peter | Brownian motors: noisy transport far from equilibrium[END_REF] to compute the mean value and a similar expression, which depends on the three mentioned complex quantities z, - 

K θ = K r (m) trials D(1) 2 (m 2 ) σ 2 D(

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., characterized by K θ = K r = 0.01m and by K θ = K r = 1m. Figure 3 displays the trajectory generated without noise (blue line) together with five trajectories obtained by setting K θ = K r = 0.01m. For the trajectory without noise, the final D(1) 2 is equal to 0.1021 m 2 . In this case we directly use the expressions in [START_REF] Martinelli | The odometry error of a mobile robot with a synchronous drive system[END_REF] in [START_REF] Park | Diffusion-based motion planning for a nonholonomic flexible needle model[END_REF] and in appendix C by numerically computing all the integrals. We obtain the following values: 2 reports the values of D(1) 2 and σ 2 D(1) 2 obtained by running 10 3 , 10 4 and 10 5 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.

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 

K θ = K r (m) trials D(1) 2 (m 2 ) σ 2 D(

INRIA

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 timedependent 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.

A Computation of u(s) p v(s) q

We have:

u(s) p v(s) q = lim N →∞ j1•••jpk1•••kq e i(θj 1 +•••+θj p -θ k 1 -•••-θ kq ) (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., δs p+q C 0 + δs p+q-1 C 1 + δs p+q-2 C 2 + • • • + δsC p+q-1 + C p+q . We trivially have C 0 = 1. Concerning the remaining coefficients, we remark, first of all, that only the ones with even index are different from 0. In other words,

C 2n+1 = 0, n = 0, 1, • • • , p+q 2 
(where x is the largest integer not greater than x). Let us compute C 2n . 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 j a ↔ j a , a, a = 1, • • • , p and with respect to the change

k b ↔ k b , b, b = 1, • • • , q but it
is not symmetric with respect to the change j a ↔ k b , 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 2n 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 [START_REF] Risken | The Fokker-Planck Equation Methods of Solution and Applications[END_REF] 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 [START_REF] Risken | The Fokker-Planck Equation Methods of Solution and Applications[END_REF] as follows:

u(s) p v(s) q = lim N →∞ {j}p{k}q p+q 2 n=0 δs p+q-2n (25) 
2n l=0 p 2n -l q l min(l,2n-l) m=0 odd/even 2n -l m l m m! (2n -l -m -1)!!(l -m -1)!! 2 j1 δ j1j2 2 j3 δ j3j4 • • • 2 j 2n-l-m-1 δ j 2n-l-m-1 j 2n-l-m 2 k1 δ k1k2 2 k3 δ k3k4 • • • INRIA 2 k l-m-1 δ k l-m-1 k l-m 2 k l-m+1 δ k l-m+1 j 2n-l-m+1 • • • 2 k l δ k l j 2n-l e i(θj 1 +•••+θj p -θ k 1 -•••-θ kq )
where, for the brevity sake, we denoted by {j}p{k}q the sum

N j1•••jpk1•••kq=1
. Each average 2 provides K r δs. Hence, all together, they provide K n r δs n . 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 [START_REF] Sokolov | Concentration Dependence of the Collective Dynamics of Swimming Bacteria[END_REF] we obtain:

u(s) p v(s) q = p+q 2 n=0 K n r 2n l=0 p 2n -l q l min(l,2n-l) m=0 odd/even (26) 2n -l m l m m!(2n -l -m -1)!!(l -m -1)!! s m lim N →∞ {j s }σ{j d }ρ{k s }χ{k d }η δs β exp i[θ j s 1 + • • • + θ j s σ + 2(θ j d 1 + + • • • + θ j d ρ ) -(θ k s 1 + • • • + θ k s χ ) -2(θ k d 1 + • • • + θ k d η )] 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 {j s }σ{j d }ρ{k s }χ{k d }η . Specifically, we consider the new indexes j s j d k s k d which differ from j s j d k s k d 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 {j s }σ{j d }ρ{k s }χ{k d }η → ρ!σ!η!χ! {j s }σ{j d }ρ{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: According to this, the last sum in [START_REF] Tierno | Magnetically actuated colloidal microswimmers[END_REF] can be written as follows:

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 w 1 < w 2 < • • • < w β . Additionally,
e iαθ0 ρ!σ!η!χ! c {w} β δs β • • exp i β-1 b=0 (α + Φ b )(∆ w b+1 w b + ∆ w b+1 w b )
Hence, we have:

u(s) p v(s) q = p+q 2 n=0 K n r 2n l=0 p 2n -l q l • (27) 
• min(l,2n-l)

m=0 odd/even 2n -l m l m m!(2n -l -m -1)!!(l -m -1)!!• •s m e iαθ0 ρ!σ!η!χ! c lim N →∞ {w} β δs β • • exp i β-1 b=0 (α + Φ b )(∆ w b+1 w b + ∆ w b+1 w b )
Now we can compute the average since we were able to separate all the independent quantities. We have:

exp i β-1 b=0 (α + Φ b )(∆ w b+1 w b + ∆ w b+1 w b ) = (28) = exp i β-1 b=0 (α + Φ b )∆ w b+1 w b • • β-1 b=0 exp i (α + Φ b )∆ w b+1 w b = = exp i β-1 b=0 (α + Φ b )∆ w b+1 w b • • β-1 b=0 exp - (α + Φ b ) 2 (w b+1 -w b )K θ δs 2 = e β-1 b=0 i(α+Φ b )∆ w b+1 w b - (α+Φ b ) 2 (w b+1 b )K θ δs 2 
We used the equality (η = N (0, 1)):

e Aη = e A 2 2 (29) 
INRIA with A = i(α + Φ b )σ b+1 and σ 2 b+1 = (w b+1 -w b )K θ δs. By substituting equation ( 28) in [START_REF] Tierno | Controlled swimming in confined fluids of magnetically actuated colloidal rotors[END_REF] and by taking the limit (N → ∞), we finally obtain:

u(s) p v(s) q = p+q 2 n=0 K n r 2n l=0 p 2n -l q l • (30) 
• min(l,2n-l)

m=0 odd/even 2n -l m l m m!(2n -l -m -1)!!(l -m -1)!!• •s m e iαθ0 ρ!σ!η!χ! c s 0 ds 1 s s1 ds 2 • • • s s β-1 ds β e β-1 b=0 i(α+Φ b )[θ(s b+1 )-θ(s b )]- (α+Φ b ) 2 (s b+1 -s b )K θ 2
which coincides with [START_REF] Niwa | Self-organizing dynamic model offish schooling[END_REF].

B Computation of u(s) p v(s) q θr

This computation follows the same initial steps carried out in appendix A. We obtain an expression equal to the one given in ( 27) but the mean value at the end must be replaced with exp i We use the following two standard results for a normal distribution:

∆ w β+1 w β γ β+1 = 0 γ β+1 odd σ γ β+1 β+1 (γ β+1 -1)!! γ β+1 even (32) and exp i(α + Φ b )∆ w b+1 w b ∆ w b+1 w b γ b+1 = (33) = σ γ b+1 b+1 exp - (α + Φ b ) 2 σ 2 b+1 2 • • γ b+1 2 a=0 γ b+1 ! (i(α + Φ b )σ b+1 ) γ b+1 -2a a!(γ b+1 -2a)!2 a
with σ 2 b+1 ≡ K θ δs(w b+1 -w b ). Equation ( 33) is obtained starting from (29) and by differentiating γ b+1 times with respect to A. Hence, the expression in [START_REF] Zhou | Planning for Noise-Induced Trajectory Bias in Nonholonomic Robots with Uncertainty[END_REF], becomes:

θr β-1 b=0 exp i(α + Φ b )∆ w b+1 w b = (34) = γ r!σ γ β+1 β+1 (γ β+1 -1)!! β b=0 γ b+1 ! β-1 b=0 σ γ b+1 b+1 e - (α+Φ b ) 2 σ 2 b+1 2 • • γ b+1 2 a=0 γ b+1 ! (i(α + Φ b )σ b+1 ) γ b+1 -2a a!(γ b+1 -2a)!2 a   
where now the sum over γ only includes the vectors γ whose last entry is even. By taking the limit N → ∞, this expression becomes:

γ r!(γ β+1 -1)!! β b=0 γ b+1 ! (s -s β ) γ β+1 2 K r 2 θ • • β-1 b=0 (s b+1 -s b ) γ b+1 2 e -(α+Φ b ) 2 K θ (s b+1 -s b ) 2 • • γ b+1 2 a=0 γ b+1 ! i(α + Φ b ) K θ (s b+1 -s b ) γ b+1 -2a a!(γ b+1 -2a)!2 a     
and we finally obtain: u(s) p v(s) q θ(s) r = where θ j ≡ θ(s j ), 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 C 120 is the conjugate of C 100 . Hence, their sum is real and it is: 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: 

C 220 = 2K 2

3K θ 2 +

 2 iµ 0 and -2K θ + 2iµ 0 , to compute the variance. We obtain the following values: D(1) 2 = 0.0680 m 2 and σ 2 D(1) 2 = 0.0012 m 4 when K θ = K r = 0.01m and D(1) 2 = 1.1130 m 2 and σ 2 D(1) 2 = 1.3052 m 4 when K θ = K r = 1m. We ran 10 5 trials and we computed from them the mean value and the variance.

Figure 2

 2 displays the values of D(1) 2 obtained in the first 10 4 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 10 4 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 θ = K r = 0.01m and the image on the right to the case K θ = K r = 1m.

Figure 2 :

 2 Figure 2: Constant ratio between linear and angular speed. Values of D(1) 2 for 10 4 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 θ = K r = 0.01m and K θ = K r = 1m in the left and right image, respectively.

Figure 3 :

 3 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 θ = K r = 0) and the other lines are five trajectories obtained by setting K θ = K r = 0.01m.

Figure 4 :

 4 Figure 4: Varying ratio between linear and angular speed. Values of D(1) 2 for 10 4 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 θ = K r = 0.01m and K θ = K r = 1m in the left and right image, respectively.

  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 We introduce the vector c ≡ [c 1 , • • • , 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, RR n°8465 we define Φ b ≡ b a=1 c a . Note that Φ β = -α.

  b ) • θr . The deterministic part can be factorized out the mean value. We need to calculate: α + Φ b )∆ w b+1 w b By using the multinomial theorem we can write: 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: β+1 b=1 γ b = r. Hence, we obtain θr

2 e i β- 1 b=0 1 b=0(s b+1 -s b ) γ b+1 2 e 2 • γ b+1 2 a=0γCds 4 e 2 cos -θ 1 - 2 cos -θ 1 +

 21122242121 2n -l -m -1)!!(l -m -1)!! •s m e iαθ0 ρ!σ!η!χ! c γ r!(γ β+1 -1)!!K (α+Φ b )[θ(s b+1 )-θ(s b )] β--(α+Φ b ) 2 K θ (s b+1 -s b ) b+1 ! i(α + Φ b ) K θ (s b+1 -s b ) γ b+1 -2a a!(γ b+1 -2a)!2 a Computation of D(s) 4 = u(s) 2 w(s) 2We use equation (17) with p = q = 2 to obtain the analytical expression of D(s) 4 = u(s) 2 w(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: -K θ (-s 1 -3s 2 +3s 3 +s 4 ) θ 2 + θ 3 + θ 4 + +e -K θ (-s 1 +s 2 -s 3 +s 4 ) θ 2 -θ 3 + θ 4 + cos -θ 1 + θ 2 + θ 3 -θ 4

ds 3 e 2 cos 2θ 1 -ds 2 e 2 cos -θ 1 + θ 2 Finally

 321222 -K θ (-4s 1 +3s 2 +s 3 ) θ 2 -θ 3 + +e -K θ (-s 1 +s 3 ) 2 cos -θ 1 + 2θ 2 -θ 3 + e -K θ (-s 1 -3s 2 +4s 3 ) 2 cos -θ 1 -θ 2 + 2θ 3The last term with n = 1 is the one with l = m = 1. By a direct computation we obtain:C 111 = 8K r s -K θ (-s 1 +s 2 ), 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

ds 2 e 2 INRIA 1 3 . 2 8 4 13 A 16 B 19 C

 221328131619 -2K θ (-s1+s2) cos -2θ 1 + 2θ 2 C 222 = 2K2 r s First and second-order statistics . . . . . . . . . . . . . . . . . . . 6 Computation of any-order moment . . . . . . . . . . . . . . . . . Diffusivity for constant and for arbitrary time-dependent linear and angular speed 9 4.1 Constant ratio between linear and angular speed . . . . . . . . . 10 4.2 Arbitrary time-dependent ratio between linear and angular speed 125 Conclusion Computation of u(s) p v(s) q Computation of u(s) p v(s) q θr Computation of D(s) 4 = u(s) 2 w(s)2 

Table 1 :

 1 Constant ratio between linear and angular speed. Values of D(1) 2 and σ 2 D(1) 2 obtained by running 10 3 , 10 4 and 10 5 trials.Table1reports the values of D(1) 2 and σ 2 D(1) 2 obtained by running 10 3 , 10 4 and 10 5 trials. It is possible to see that, as the number of trials increases, the results converge to the values obtained by using our formulas.

	0.01	10 3	0.0688	0.0011
	0.01	10 4	0.0664	0.0012
	0.01	10 5	0.0679	0.0012
	1	10 3	1.1189	1.3736
	1	10 4	1.1124	1.2767
	1	10 5	1.1126	1.3035

1) 2 (m 4 )

Table 2 :

 2 Varying ratio between linear and angular speed. Values of D(1) 2 and σ 2 D(1) 2 obtained by running 10 3 , 10 4 and 10 5 trials.

	1) 2 (m 4 )

RR n°8465

Note that, deriving any-order moment, corresponds to analytically derive the probability density distribution[START_REF] Papoulis | Probability, Random Variables, and Stochastic Processes[END_REF].INRIA

Note that this constraint does not limit the robot motion: for instance, a go and back motion is easily obtained by accomplishing 180deg rotations RR n°8465