Iván De Jesús

Salgado Ramos

Andrey Polyakov

Signal differentiation by means of algebraic techniques and sliding mode techniques

The numerical di¤erentiation has been considered an important problem to deal in control, signal processing, applied mathematics, among others. Several techniques based on Luenberger observers, homogeneity, algebraic techniques have been developed considering the problems that can be associated to a real signal. In this manuscript, a qualitative study was made to show the behavior of di¤erent di¤erentiation techniques. The di¤erentiation algorithms were based mainly in the algebraic techniques and polynomial approximation. A second approach was selected as the well known homogeneous di¤erentiator. Four di¤erent di¤erentiator techniques were studied: the algebraic di¤erentiator, a di¤erentiator based on Legendre polynomials, a di¤erentiator based on Jacobi polynomials and the Homogeneous di¤erentiator. The problems inherit to signal di¤erentiation included in the analysis are the algorithm complexity and algorithm accuracy under noise and quantization.

Introduction

Numerical di¤erentiation (ND) is the computation of values of the derivative of a function f only from some given values of this function. It has been well known that numerical di¤erentiation should be avoided whenever possible because it tends to make matters rough and generally give values of d dt f much less accurate than those of f . It has been reported that unstable characteristics of numerical di¤erentiation process may lead to large errors in the computed results [Choi, 2009]. Unfortunately, in many applications it is necessary to estimate the derivative of a function given the noisy values of it. Actually, in real-time applications most of errors come from round-o¤ and truncation.

ND is ill-posed in the sense that a small error in measurement data can induce a large error in the approximate derivatives. Therefore, various numerical methods have been developed to obtain stable algorithms more or less sensitive to additive noise. Several results has been obtained from different theories like observer design in the control literature, digital …ltering in signal processing and system identi…cation [Dayan, 2011]. According to [Dayan, 2011], most of signal di¤erentiators fall into 8 main categories, that is, the …nite di¤erence methods, the Savitzky Golay methods, the wavelet di¤erentiation methods, the Fourier transform methods [START_REF] Jauberteau | [END_REF], the molli…cation methods, the Tikhonov regularization, the algebraic methods and the di¤erentiation by integration methods. By this kind of problems, an analysis of di¤erent di¤erentiation techniques under the presences of noises and with the problem of quantization is made in this manuscript [Levant, 1998], [Liu et al., 2014], [START_REF] Mboup | [END_REF]. The rest of the manuscript is organized as follow, in the next section, the di¤erentiation techniques are explained, including, the Euler approximation, the homogeneous di¤erentiator, algebraic di¤erentiators and the di¤erentiators based on the weak derivative. Then, in section III, simulation results are presented, the di¤erent di¤erentiators are compared in terms of computational aspects (like complexity, memory size among others) and accuracy and robustness in the presence of noisy measurements and quantization. Finally in section IV, some conclusions are given.

Design of di¤erentiators

The problem of numerical di¤erentiation can be summarized as:

Problem 1 Given a signal y (t), to estimate _ y (t) with t 2 R; using the information on the signal y (s) for s 2 [t T; t], where T is the size of the window.

The information on y (t) is obtained on-line in a discrete time, i.e. y (t i ) is known for t i = ih, i = 1; 2; 3::: and h > 0 is the sampled period.

The signal y (t) is corrupted by measurements noises whose statistics are unknown

The model of the system that generates the output signal y (t) is unavailable.

Euler Di¤erentiation

The simplest method for obtaining the …rst order time derivative of a signal y (t), is the so called …nite di¤erence method. It consist in the following approximation to the derivative _ y (t) in the time instant t i _ y e (t i ) = y (t) y (t i ) t t i This method is quite sensitive to the presence of noise perturbations in the signal to be processed.

Homogeneous di¤erentiator

For obtaining the …rst derivative of a signal, the Homogeneous di¤erentiator [START_REF] Perruquetti | [END_REF], [Polyakov et al., 2014] is described by the following equation

dx 1 dt = k 1 jx 1 yj 1 1+ sign (x 1 y) + x2 dx 2 dt = k 2 jx 1 yj 1 2 1+ sign (x 1 y)
The well known linear high-gain di¤erentiator and the exact high order sliding mode di¤erentiator [Levant, 1998] can be obtained as a partial cases of the homogeneous di¤erentiator, when the homogeneity degree tends to 0 and -1 respectively.

Algebraic Signal di¤erentiation (Jacobi approximation)

The algebraic di¤erentiator [START_REF] Mboup | [END_REF], [Liu et al., 2014] involving Jacobi Polynomials and algebraic methods is describred as follows. Consider a noise-corrupted measurement ( ) is the i th order shifted Jacobi orthogonal polynomial (See the appendix).

x $ : I ! / R, x $ (t) = x (t) + $ (t)
In order to approximate x on [t o T; t o ], the classical polynomial approximation by taking the N + 1 …rst terms in the Jacobi series expansion given in (2.1) was considered, the obtained N th order polynomial is denoted by D (0) ; ;T;N x (t o T ). Thus,

8 2 [0; 1] , D (0) ; ;T;N x (t 0 T ) = N X i=0 D P ( ; ) i ( );x(t 0 T ) E ; P ( ; ) i 2 ; P ( ; ) i ( ) :
(2.2)

Hence, the n th order derivative of x can be approximated by the one of D

; ;T;N x (t 0 T ). Then, the (causal) Jacobi di¤erentiator is de…ned as follows

D (n) ; ;T;q x (t 0 T ) := 1 ( T ) n d n d n n D (0) ; ;T;N x (t 0 T ) o ;
where q = N n 2 N. This di¤erentiator can also be obtained by taking the q + 1 …rst terms in the Jacobi series expansion of x (n) by a q th order polynomial on [t 0 T; t 0 ] :

8 2 [0; 1] , D (n)
; ;T;q x (t 0 T ) = Moreover, it can be given by the following integral formula:

8 2 [0; 1] , D (n)
; ;T;q x (t 0 T ) = 1 ( T ) n 1 Z 0 Q ; ;n;q; ( ) x (t 0 T ) d ; (2.3)
where

Q ; ;n;q; ( ) = ! ; ( ) q X i=0 C ; ;n;i P ( +n; +n) i ( ) P ( ; ) n+i ( ) with C ; ;n;i = ( + + 2n + i) ( + + 2n + i) (n + i + 1) ( + n + i + 1) ( + n + i + 1) :
Finally, we substitute x in (2.3) by x $ so as to obtain the Jacobi di¤erentiator D ($)

; ;T;q x (t 0 T ) in the noisy case.

Weak derivative based approach

In the following de…nition the concept of weak derivative is described,

De…nition 2 (Clasiccal) D y 2 L 2 [a;b] is the weak derivative of y (t) 2 L 2 [a;b] if and only if hD y ; i = D y; _ E ; 8 2 C 1 [a;b] (a) = (b) = 0 (2.4)
Remark 3 From integration by parts

hD y ; i = y ( ) ( )j b a D y; _ E Note. h ; i is the scalar product in L 2 .
Proposition 4 The weak derivative of a signal y (t) 2 L 2 [a;b] is given by

Dy = 1 X i=0 c i '(1) i (t) Proof. De…ne a function S N = N X i=0 c i '(1) i (t), where, if N ! 1, S N ! Dy and hS N ; i ! hD y ; i , 8 2 C 1 [ 1;1]
Then, by the de…nition of the internal product [Hildebrand, 1987], one have

* N X i=0 c i '(1) i ; + = N X i=0 c i D '(1) i ; E Integrating by parts the next result is obtained N X i=0 c i D '(1) i ; E = N X i=0 c i '(1) i b a N X i=0 c i D 'i ; (1)

E

By the de…nition of ( ) in (2.4) the result is obtained.

Di¤erentiation of a signal from its discrete measurements by Legendre polynomials

Consider the function y (t) 2 L 2 [ 1;1] , the main problem to deal is to approximate the …rst order derivative of y (t) taking at discrete instans of time by means of a series of orthonormal functions de…ned as

'i (t) 2 L 2 [ 1;1] ; i = 1; 2; ::: (2.5)
The function y (t) could be approximated by

y (t) 1 X i=0 c i 'i (t) (2.6)
The derivative of the function y (t) by de…nition is given by

d dt y (t) := 1 X i=0 c i d dt 'i (t) = 1 X i=0 hy; 'i i d dt 'i (t) (2.7)
where, the dot product is given as

hx; yi = 1 Z 1 x ( ) y ( ) d (2.8)
with this de…nition, equation (2.7) can be rewritten as

d dt y (t) = 1 X i=0 1 Z 1 y ( ) 'i ( ) d dt 'i (t) d d dt y (t) = 1 Z 1 y ( ) 1 X i=0 'i ( ) d dt 'i (t) d
if the function is approximated with N orthonormal functions, one has

d dt y (t) 1 Z 1 y ( ) N X i=0 'i ( ) d dt 'i (t) d where 'i = 1 h' i ; ' i i 1=2 ' i (2.9)
where ' i can be any orthogonal basis. The derivative approximation in the point y (1) (by the windowing to be described below) is given by

d dt y (1) 1 = Z 1 y ( ) N X i=0 'i ( ) d dt 'i (1) d (2.10) De…ning i = d dt 'i (1) ; equation (2.10) takes the form d dt y (1) = 1 Z 1 y ( ) N X i=1 i 'i ( ) d (2.11)
The polynomials used in this approximation were selected as Legendre polynomials that are orthogonal functions over the interval [ 1; 1]. In this way, for signals de…ned in a di¤erent interval [a; b], the following transformation in time must be implemented

t = s a b a s b b a (2.12) Then, s = (b a) t + a + b 2 (2.13)
Applying the last equation, the following approximation is obtained

d dt y (1) 1 Z 1 y (s) N X i=1 i 'i ( ) d (2.14)
For the discrete-time domain, taking a windowing of Ñ samples the last equation is approximated by

d dt y (1) Ñ 1 X j=1 j+1 Z j ỹ (s) N X i=1 i 'i ( ) d = Ñ 1 X j=1 y (s) j+1 Z j N X i=1 i 'i ( ) d = Ñ 1 X j=1 y j N X i=1 i j+1 Z j 'i ( ) d (2.15) De…ning ij = j+1 Z j 'i ( ) d (2.16) valid for i i+1 with i = 1 + 2h b a i = 1 + 2 Ñ i (2.17)
h is the sampled period and q j is de…ned as

q j = N X i=1 i ij (2.18)
Finally, the derivative approximation takes de form

d dt y (1) Ñ 1 X j=0 y j q j (2.19) 2.4.2 Example for N = 3
From equation (2.9) the parameters 'i were obtaining for N = 3, for ' 0 one gets

h' 0 ; ' 0 i = 1 Z 1 d = 2, '0 = 1 p 2 (2.20) For ' 1 one gets h' 1 ; ' 1 i = 1 Z 1 2 d = 3 3 1 1 = 2 3 '1 = r 3 2 (2.21)
For ' 2 one gets h' 2 ; '

2 i = 1 Z 1 1 4 (3 2 1) 2 d = 1 4 1 Z 1 (9 4 6 2 + 1) d = = 1 4 9 5 5 2 3 + 1 1 = 1 4 9 5 (1) 5 2 (1) 3 + 1 1 4 9 5 ( 1) 5 2 ( 1) 3 1 = 2 5 '2 = r 5 2 1 2 (3 2 1) (2.22)
Finally, for ' 3 one gets

h' 3 ; ' 3 i = 1 Z 1 1 4 (5 3 3 ) 2 d = 1 4 1 Z 1 (25 6 30 4 + 9 2 ) d = 1 4 25 7 7 6 5 + 3 3 1 1 ! = 1 4 25 7 6 + 3 1 4 25 7 + 6 3 = 1 4 50 7 12 + 6 = 2 7 '3 = r 7 2 ' 2 ( ) = r 7 2 1 2 (5 3 3 ) (2.23)
The above calculations let us to the values of i

0 = 0 1 = r 3 2 2 = 3 r 5 2 3 = 6 r 7 2
(2.24)

Weak-derivatives of periodic functions

Optimal estimation in L 2

The objective is to minimize

Dy Dy L 2 ! min Dy2M (2.25)
where M is a linear subspace of L 2 .

Remark 7 The optimal estimate of Dy is the unique solution of (2.25) if and only if D Dy Dy; z

E = 0, 8z 2 M
Hence, if M = span f' 0 ; : : : ; ' N g where ' i are orthonormal, that is

' i ; ' j = 1; i = j 0; i 6 = j then Dy = N X i=0 D Dy; ' i E ' i 2.
5.2 Estimation of weak derivatives of periodic functions

1. y 2 W k;2 [a;b] ; then 9D (k) y in the sense of de…nition 2. 2. If f' i g 1 i=0 is orthonormal basis in L 2 [a;b] , such that ' i 2 C 1 [a;b] and ' (j) i (a) = ' (j) i (b) for i = 1; : : : ; k 1, then D(k) y = N X i=0 D D(k) y; ' i E ' i (2.26)
is optimal estimation of D (k) y in the subspace span f' 0 ; : : : ; ' N g : 3. Since '

(k 1) i (a) = ' (k 1) i (b)
, then (2.26) can be rewritten as

D(k) y = ( 1) k N X i=0 D y; ' (k) i E ' i (2.27)

Estimation of the …rst derivative

Consider the function

y 2 W 1;2 L 2 [0;1]
Assume that y can be extended into a periodic function ỹ 2 W 1;2 [ 1;1] de…ned as (see …gure )

ỹ ( ) = y ( ) 2 [ 1; 0] y ( ) 2 [0; 1] (2.28) Let ' i be orthonormal basis of periodic functions in L 2 [ 1;1] such that ' i 2 C 1 [ 1;1] : ' ( 1) = ' (1)
Then by (2.27) we have

Dy ( ) = N X i=0 hy; _ ' i i ' i ( ) , 2 [ 1; 1]
is the optimal estimate of Dỹ in the subspace span f' 0 ; : : : ; ' N g, N 0.

Comparison between signal di¤erentiators

The function used in simulations (…gure 3.1) was chosen as 15t) the available output was considered. i. e.

y (t) = e sin(

Di¤erent sampled times:

Di¤erent levels of quantization Two kinds of noise in the output

-Uniform distributed noise from [ 1; 1]
-Deterministic noise ( (t) = 0:005 sin 900t + 0:01 cos 600t)

The parameters used in simulation for each di¤erentiator were chosen as:

Homogeneous Di¤erentiator = 1 2 k 1 = 180 k 2 = 350 Algebraic Di¤erentiator (Jacobi) = 11 = 0:1 q = 4
where q is the number of polinomios used. The parameter selection was done using the notation presented in table 1 in the appendix.

Algebraic Di¤erentiator (Legendre) q = 4 and Weak di¤erentiator ' i ( ) = sin (i ) ; i = 1 : q q = 8 = 3 in the same way, q is the number of elements in the cosines basis, is the delay for computing the derivative. For the di¤erentiators based on algebraic The parameters for the Jacobi and weak di¤erentiators were chosen in order to obtain a good approximation of the signal derivative without having a big delay. By this reason the number of elements in the basis is not equal in all the cases.

For comparing the behavior of each di¤erentiation technique, several scripts were run in Matlab in a computer with the following characteristics: Intel Core i7-4700MQ @ 2.40 MHz. 16 Gb RAM, Operating system of 64 bits.

The …rst analysis was done in terms of computational aspects like complexity, memory size and complexity tuning of each di¤erentiators The second part of the analysis was done in terms of accuracy, where several sampled times were tested in order to know what di¤erentiator has a better performance in presence of bounded perturbations and quantization. The objective of all the simulations was to approximate the derivative of an arbitrary unknown signal.

Computational aspects

The computational aspects involve complexity of the algorithm, the memory required to approximated one point of the signal in time and the complexity of tuning the algorithm.

Complexity

In this study it was de…ned as the number of basic operations that have to be implemented in order to obtain the derivative in an instant of time. In the following table the number of operations required for each di¤erentiator to calculate one sampled time of simulation is summarized. A second characteristic to know the complexity of the algorithm was selected as the time that the computer needed to obtain one second of simulation.

In the following table, this analysis is described with being the sampled period. The results are given for whole interval Homog. k 1 , k 2 , and Algebraic (Jacobi) q, , Algebraic (Legendre) q Weak q and 3.2 Accuracy 3.2.1 Noisy free.

7:35x10

1. Sampled time h = 0:001; m = 100:

In …gure (3.2), the simulations results obtained for each di¤erentiator are shown. When the signal is not a¤ected by any noise in its output, the four di¤erentiators can estimate the derivative of the nonlinear signal when it is sampled with small sampled time. The homogeneous approximation and the Euler approximation converged faster than the other algorithms. The best accuracy was obtained when the Euler methods is applied for small sampled period. In …gure (3.3) the estimation error de…ned as the absolute value of the error is shown. Clearly, a less error is appreciated with the Euler approximation.

2. Sampled time h = 0:002, T = 50:

When the sampled time is selected as 0.002 and the window is reduced according to this new sampled time …gure (3.4), the di¤erent di¤erentiators still reproduce the derivative error. However, error in the Euler approximation increased, The estimation errors are shown in …gure (3.5). The Jacobi di¤erentiator without any modi…cation in its parameters becomes to be inaccurate. The three di¤erentiators based on weak derivative still having good approximation properties.

3. Sampled time h = 0:01; T = 10 A third example was chosen taking a smaller sampled time. See …gure (3). Also the window was reduced to a value of 10. In this way, the Jacobi di¤erentiator without any modi…cation of its parameters lose its accuracy. The best estimation under these conditiones was obtained with the weak de…rrentiator with periodic functions with a …rst order hold in the available output. The di¤erent estimation errors are described in …gure (3.6).:

Quantization.

A second test was made including the problem of quantization. In this analysis three di¤erent levels of quantization were tested. In this section the signal was considered without any noise. Two levels of quantization were tested: 0.005 and 0.02 were etsd 1. Level of quantization of 0.005, h = 0:001 and T = 100. In …gure (3.7) the derivative estimation of the same nonlineal signal is depicted. Even when the level of quantization was chosen small enough in order to obtain a good approximation the Euler approximation and the gomogeneous di¤erentiator becomes to have some inaccurate in the estimation 2. Quantization of 0.02, h = 0:001 and T = 100. When the signal is quantized only with 50 values between two consecutive integers, the Euler approximation and the homogeneous approximation did not estimate the trajectories of the nonlinear function (see …gure 3.9). The accuracy of the algebraic approaches decreased with this level of quantization, however the approximation remains acceptable. They are more robust to this problem. This fact can be appreciated in …gure (3.10). Two classes of noise where tested to see the robustness of each di¤erentiator. Distributed uniform noise and deterministic noises were applied. The simulations were obtained for di¤erent cases, …rst, each di¤erentiator was tested only with noise in the measurements without considering the problem of quantization. Then, a second group of simulations was made adding the problem of quantization to the noisy signal.

Signal with an uniform distributed noise 1. Signal without quantization and h = 0:001. The results for each di¤erentiator are shown in …gure (3.11). In …gure (3.12) the estimation errors are described. The Euler approximation lost accuracy in the presence of noises. The polynomial approximations, as they were calculated using the …rst 100 values to calculate one instant of time produce a less oscillated signal.

2. Signal with quantization of 0.02 and h = 0:001. In …gure (3.13), the results for the noisy signal with quantization of 0.001 can be appreciated, the results are quite similar to the case without quantization.

The less a¤ected algorithms are the ones based in algebraic techniques.

In …gure (3.14) the estimation errors are shown.

Deterministic noise

1. Signal without quantization and h = 0:001: In …gure (3.15) the simulations results for the di¤erentiator algorithm a¤ected by a deterministic noise described by ( (t) = 0:005 sin 900t + 0:01 cos 600t are presented. The four algorithms reproduced the derivative. However, the Euler and homogeneous algorithms were less robust than the other approximators. 2. Signal with quantization of 0.02 and sampling period 0.001. In …gure (3.17) the simulation results for the di¤erentiators are shown when the signal is quantized. This small quantization value a¤ect so much the behavior of the four di¤erentiators in comparison with the previous simulation when the signal was not quantized. The estimation error is described in …gure (3.18).

Conclusions

In this manuscript, several di¤erentiation techniques were tested in order to see the advantages and disadvantages of each technique. Four di¤erent techniques were implemented considering several problems that always are presented in real applications. The techniques based on algebraic techniques have better performance in the presence of noise in the measurements. Moreover, if the problem of quantization is included in the signal di¤erentiation, the di¤erentiators based on algebraic technique also bring some better results. The fourth di¤erentiation technique was selected as the homogeneous di¤erentiation, selecting its parameters in order to have a sliding mode differentiator.

Between the algebraic techniques, the one based on Jacobi polynomials showed better behavior, however, it is important to note, that this technique is a generalization of the Legendre di¤erentiator and Tchevyshev di¤erentiator. However, the parameter induced a delay in the di¤erentiation process. If a Fourier basis is used, the estimation of the …rst order derivative presents a good approximation when the signal is not corrupted by noise. In the followin tables, the …nal conclusions are given.

  b] be the Sobolev space of absolute continuous periodic functions: y (a) = y (b) if there exist 2 C 1 [a;b] : (a) = (b), by the de…nition of (2.8) and integrating by parts we have h _ y; i = y ( ) ( )j 1 [a;b] : (a) = (b) then, Dy is called the weak derivative of periodic function y 2 W 1;2 [a;b] : y (a) = y (b). Remark 6 If y 2 C k [a;b] ; then Dy = _ y .

Figure 2

 2 Figure 2.1: Transformation of signal y into the ỹ ( ) in the interval [ 1; 1]

Figure 3

 3 Figure 3.1: Function chosen in simulation

Figure 3 . 2 :

 32 Figure 3.2: Derivative estimation sampled at 0.001 seconds.

Figure 3 . 3 :

 33 Figure 3.3: Derivative estimation error with sampled time of 0.001

Figure 3 . 6 :

 36 Figure 3.6: Derivative estimation error

Figure 3 . 8 :

 38 Figure 3.8: Derivative estimation error for quantization level of 0.002 and sampled period of 0.001

Figure ( 3

 3 .16) shows the euclidean norm of the estimation errors. The Algebraic di¤erentiators almost have the same performance.

Figure 3 .

 3 Figure 3.11: Derivative estimation with the signal corrupted by distributed uniform noise from -0.1 to 0.01 and sampled period of 0.001

Figure 3 .

 3 Figure 3.16: Derivative estimation error with the signal corrupted by deterministic noise, sampled period of 0.001

  , where I is a …nite time open interval of / R + , x 2 C n (I) with n 2 N, and $ is an additive corrupting noise. The objective is to estimate the n th derivative of x using x $ . In the following, a class of algebraic di¤erentiator involving Jacobi polynomials is introduced.For any t o 2 I, D t 0 := t 2 / R + ; t 0 t 2 I , x can be locally given on [t o T; t o ] with T 2 D t 0 by the following Jacobi orthogonal series expansion

	D ; ;T;q x (t 0 T ) =	q X i=0	D	P i ( +1; +1) i P ( +1; +1) ( ); _ x(t 0 T ) +1; +1 E +n; +n 2	P	( +1; +1) i	( ) :	(2.1)
	where P	( ; ) i

Remark 8 If ' i ( 1) = ' i (1) = 0, Dy (1) = 0, then Dy ( 1) is the estimate of the derivative of = 1: Remark If ' i ( 1) = ' i (1), Dy(1) = 0, then Dy( 1) is the estimate of the derivative of = 1.

Hence

and by the property given in (2.28), the next result is obtained

Discrete estimation of the …rst derivative

Consider the equation

for the discrete time domain, taking a windowing of Ñ samples, the last equation is approximated by

h is the sampled period. With q j de…ned as

(2.33)

The optimal derivative approximation takes the form

First Order Approximation Consider the approximation given by

where the function y ( ) is de…ned as

Applying this output as the information of the di¤erentiator in the discrete time domain, the equation (2.34) becomes

because ' i ( ) do not depend of one has

y ( j ) and y ( j+1 ) are constants during the interval [ j ; j+1 ] and they do not depend of the variable , then

and they can be obtained integrating by parts, that is,

Appendix

Legendre polynomials

In mathematics, an orthogonal polynomial sequence is a family of polynomials such that any two di¤erent polynomials in the sequence are orthogonal to each other under some inner product. The most widely used orthogonal polynomials are the classical orthogonal polynomials, consisting of the Hermite polynomials, the Laguerre polynomials, the Jacobi polynomials together with their special cases the Gegenbauer polynomials, the Chebyshev polynomials, and the Legendre polynomials [Hildebrand, 1987].

The Legendre polynomials are given by the relation

that is often called the Rodrigues formula for P r (t). From the preceding equation, it follows that 1 Z 1 P r (t) P s (t) dt = 0, r 6 = s where r and s are nonnegative integers.

Jacobi polynomials

The Jacobi orthogonal polynomial de…ned on [0; 1] as follows


Giovanni Neglia
email: giovanni.neglia@inria.fr

Matteo Sereno
email: matteo.sereno@di.unito.it

Giuseppe Bianchi
email: giuseppe.bianchi@uniroma2.it

Geographical Load Balancing across Green Datacenters: a Mean Field Analysis *

Geographic Load Balancing" is a strategy for reducing the energy cost of data centers spreading across different terrestrial locations. In this paper, we focus on load balancing among micro-datacenters powered by renewable energy sources. We model via a Markov Chain the problem of scheduling jobs by prioritizing datacenters where renewable energy is currently available. Not finding a convenient closed form solution for the resulting chain, we use mean field techniques to derive an asymptotic approximate model which instead is shown to have an extremely simple and intuitive steady state solution. After proving, using both theoretical and discrete event simulation results, that the system performance converges to the asymptotic model for an increasing number of datacenters, we exploit the simple closed form model's solution to investigate relationships and trade-offs among the various system parameters.

Introduction

Providers such as Amazon, Google, Facebook, etc., are making a considerable effort to offer efficient, scalable, and reliable services. To achieve these goals such services need to be supported by massive datacenters and relevant infrastructures to distribute power and provide cooling. Power management is becoming a crucial issue. Not only power consumption is ever increasing with an increasing user base and service expansion, but, as pointed out by several studies, the power consumption of datacenters is largely wasted.

In this paper we consider a set of micro-datacenters which are additionally powered by renewable energy sources, e.g., photovoltaic (PV) panels. Due to the current high costs for storing energy, the best use of renewable energy is to consume it when it is produced. Hence, we would ideally wish to adapt each micro-datacenter's load to the instantaneous energy production. One way to address such goal is to federate several micro-datacenters each other, and use a central controller to dispatch jobs where renewable energy is available, so as to minimize the (non-renewable) energy cost. The possibility to manage more jobs obviously offers a higher flexibility. The law of large numbers guarantees indeed that the aggregated load will be more regular and then easier to exploit for smart load scheduling as it is the case in a big datacenter. But when local renewable sources are available, a microdatacenters' federation offers an additional advantage in comparison to a large datacenter: renewable energy production at different locations can be loosely correlated and then the aggregated energy production exhibits less variability.

Consider the following ideal case: a set of N identical datacenters, each with independent job arrival processes with rate λ and a single server with computing rate µ, and PV panels able to feed the datacenter a fraction s of the time. Compare it with a single datacenter which aggregates locally all the jobs as well as the computing and energy production infrastructure. The total normalized load for the federation of datacenter is ρ = (N λ)/(N µ) = λ/µ with a normalized variability (standard deviation of the number of working servers divided by the average number of working server) equal to (1 -ρ)/(ρN ). Similarly the federation can power through renewables a fraction s of its computing resources with a normalized variability equal to (1 -s)/(sN ), if the amounts of renewable energy produced at different datacenters can be considered independent. The single datacenter manages the same aggregate load with the same normalized variability, but the situation is different energywise. The single datacenter can be powered by renewables a fraction s of the time, but now the normalized variability is (1 -s)/s, if, as it is reasonable to assume in first approximation, all the PV panels at a given location produce (/do not produce) at the same time.

The example above is clearly over-simplified, it ignores the costs of job dispatchment among the micro-datacenters, the effect of fixed energy costs that are easier to optimize at a single datacenter, the possibility that renewable energy dynamics are too fast to be exploited by smart scheduling strategies, how revenues should be split among the datacenters, etc. Nevertheless, this example highlights a potential benefit from federating micro-datacenters, that is interesting to quantify. As we are going to show below, even simple models for job traffic and energy production lead soon to scenarios for which it is difficult to provide closed-form expressions for the energy cost of a federation of micro-datacenters. One may then need to rely on expensive simulations that hide the role played by the different parameters. For this reason, in this paper we propose a mean field (fluid) model that is asymptotically correct and allows us to derive simple formulas for the main performance metrics, like the expected energy costs of the system.

The paper is organized as follows. After a brief discussion of related work, we introduce the system model in section III, and we provide and justify with both theoretical and simulation results a mean field approximation in section IV. In section V, we exploit the resulting simple model to quantify performance and trade-offs emerging in scenarios characterized by variable renewable energy production across micro datacenters.

Related Work

In Geographic Load Balancing (GLB) systems user requests are initially accepted by front-end elements and then redirected by a scheduler to geographically distributed datacenters for processing. The scheduler's decisions may depend on several mutually interacting (and in some case conflicting) objectives such as minimizing the electricity cost, the carbon-footprint and the response time. The paper [START_REF] Qureshi | Cutting the Electric Bill for Internet-scale Systems[END_REF] is one of the first studies about GLB. In particular, it focusses the attention on the key issues fostering the use of GLB such as different energy markets (e.g., day-ahead and real-time markets), and temporal or geographical energy price variations. The GLB represents the combination of these basic ingredients with the use of energy related metrics in the scheduler decisions. In this manner it is possible to account for different workload conditions, and time and geographical variability of the electricity costs.

In the last years other studies addressed the same problem by adding different scheduling constraints and/or by optimizing different metrics (see [START_REF] Rao | Minimizing Electricity Cost: Optimization of Distributed Internet Data Centers in a Multi-Electricity-Market Environment[END_REF], [START_REF] Rao | Coordinated Energy Cost Management of Distributed Internet Data Centers in Smart Grid[END_REF], [START_REF] Yao | Data centers power reduction: A two time scale approach for delay tolerant workloads[END_REF], [START_REF] Guo | Cutting Down the Energy Cost of Geographically Distributed Cloud Data Centers by Using Energy Storage[END_REF], [START_REF] Lin | Online algorithms for geographical load balancing[END_REF], [START_REF] Liu | Greening Geographical Load Balancing[END_REF]). For instance, the papers [START_REF] Rao | Minimizing Electricity Cost: Optimization of Distributed Internet Data Centers in a Multi-Electricity-Market Environment[END_REF] and [START_REF] Liu | Greening Geographical Load Balancing[END_REF] introduce additional constraints for accounting QoS guarantees; while the interaction between GLB and smart grids, and then the exploitation of the workload demand-response capability have been addressed in [START_REF] Rao | Minimizing Electricity Cost: Optimization of Distributed Internet Data Centers in a Multi-Electricity-Market Environment[END_REF]. Furthermore, the interaction of energy storage systems and GLB has been addressed in [START_REF] Guo | Cutting Down the Energy Cost of Geographically Distributed Cloud Data Centers by Using Energy Storage[END_REF]. Indeed, storage systems can be used to smooth the variability of power supply and this is very important when the datacenters are powered by renewable sources.

Several studies pointed out that large datacenters are extremely expensive to maintain and this has encouraged the development of architectures that interconnects multiple micro-datacenters [START_REF] Bianchini | Leveraging Renewable Energy in Data Centers: Present and Future[END_REF]. This trend influenced our work because workload scheduling among a large number of interconnected datacenters gives rise to computational problems (e.g., see the summary of the techniques used in geographical load balancing in [START_REF] Rahman | A Survey on Geographic Load Balancing Based Data Center Power Management in the Smart Grid Environment[END_REF]).

The works closest to ours are [START_REF] Liu | Greening Geographical Load Balancing[END_REF] and [START_REF] Lin | Online algorithms for geographical load balancing[END_REF], where geographical load balancing is driven by time-varying energy prices, that can be due to a significant local production from renewable sources. While in these papers energy prices are considered to be known in advance over some future timehorizon, in our case renewable energy production is a stochastic process and scheduling is decided on the basis of the current state of the system.

Problem

We consider a federation of N identical micro datacenters. The aggregated job arrival process at the federation is modeled as a Poisson process with rate N λ. The service time of each job is assumed to be exponential with expected value 1/µ. 1 Each datacenter is connected to the grid but it can be powered also by some renewable source. We consider here that the renewable source can be in two states: in state S (sunny) the energy produced by the source is able to power the whole datacenter, in state C (cloudy) the energy produced is negligible. Renewable states evolve according to a continuous time Markov Chain. Let ν C and ν S denote respectively the transition rates from S to C and from C to S. The model for the renewable source can be made arbitrarily more realistic by adding multiple states. For the moment we assume that the Markov chains associated to renewable sources at different datacenters evolve independently.

When a new job arrives the scheduler dispatches it i) to a datacenter that is available to process it and in state S (i.e., currently powered by renewables) if any, otherwise ii) to an available datacenter if any, and as last option iii) to a central waiting queue from which the job will be moved to the first available datacenter. The system is then operating as an M/M/N queue with the characteristic that available servers in state S get jobs with strict higher priority than other servers. Among the work conserving disciplines this intuitively minimizes the total expected energy cost.

The system can be described as a continuous time Markov Chain with state (J N (t), S N (t), B N S (t)), where J N (t) is the number of jobs in the system, S N (t) is the number of servers in state S, and B N S (t) is the number of servers busy (i.e., serving a job) and in state S, all at time t. The Markov chain has a very particular structure: for example J N (t) itself, representing the number of jobs in a M/M/N queue, evolves as a Markov chain. S N (t) is described by a simple Markov chain too. In particular the stationary dis-tributions of J N (t) and S N (t) can be derived easily in closed-form. Despite these properties, it is not easy to characterize B N S (t) and in particular we have not been able to derive in closed-form its stationary distribution. This is less surprising if we think about a similar problem for parallel queues where the simple join-the-shortest-queue policy couples the status of the different queues so that their stationary distribution can be expressed only as an infinite mixture of geometric distributions [START_REF] Adan | Analysis of the Asymmetric Shortest Queue Problem[END_REF] (there are many works on priority queues and/or shortest queue policies, see for instance [START_REF] Foley | Join the Shortest Queue: Stability and Exact Asymptotics[END_REF] [START_REF] Harchol-Balter | Multi-server queueing systems with multiple priority classes[END_REF]). Similarly, here our dispatching policy couples the two different states of a server (being busy and being powered by renewables) in a non-trivial way so that it is difficult to characterize the process B N S (t), as we need to quantify the energetic savings coming from the federation.

In order to study the system (J N (t), S N (t), B N S (t)) we could resort to simulations or to a numerical solution of the Markov Chain. In both cases the computational cost increases with the number of datacenters N . These difficulties are aggravated if more realistic and then more complex models for traffic arrival process or renewable energy evolution are considered with a potential explosion of the state space. Moreover, the effect of the different parameters can be more difficult to unveil using numerical methods. For these reasons, as it has been successfully done in other fields, we derive the fluid limit of the Markov chain of interest, that allow us to obtain simple closed-form expressions for the main performance metrics independently from the system size N .

Fluid Model

In this section we show that the stochastic dynamics of the Markov chains (J N (t), S N (t), B N S (t)) converge in probability to a deterministic process as N diverges. 2 More precisely, we will show that if 1/N (J N (0), S N (0), B N S (0)) converges to the constant values (j 0 , s 0 , b s,0 ) when N diverges, then there exists a vector of deterministic functions (j(t), s(t), b s (t)) such that (j(0), s(0), b s (0)) = (j 0 , s 0 , b s,0 ) and for any T > 0:

sup 0≤t≤T 1 N (J N (t), S N (t), B N S (t)) -(j(t), s(t), b s (t)) → P 0,
i.e., the rescaled process converges to (j(t), s(t), b s (t)). This kind of convergence results has become popular since the seminal work of Kurtz (see for example [START_REF] Kurtz | Limit Theorems and Diffusion Approximations for Density Dependent Markov Chains[END_REF]), that shows that the limiting process can be described by a system of differential equations: dx/dt = f (x(t)), where f () is called the limiting drift function. Classic results require f () to be a Lipschitz function. By carrying out the usual derivation of the fluid limit for the process (J N (t), S N (t), B N S (t)), the corresponding function f () will appear to be discontinuous and then it has not the Lipschitz property. Nevertheless, we can apply more recent and general results from [START_REF] Gast | Markov chains with discontinuous drifts have differential inclusion limits[END_REF] to show that the dynamics converge to the solution of a system of differential inclusions, i.e. where the function f () is replaced by a set valued function.

As we observed in the previous section, the processes J N (t) and S N (t) are themselves Markov chains.

Rather than studying the joint system (J N (t), S N (t), B N S (t)) we first derive the fluid limits for J N (t) and S N (t) and then move to consider the fluid limit for B N S (t). While we could directly consider the limit of the triplet, this approach can result easier to follow for the reader unfamiliar with fluid limits. Moreover, the results for J N (t) and S N (t) do not require the more complex machinery of differential inclusions, so this approach allows us to better highlight where difficulties arise for B N S (t). The Markov chain describing J N (t) is such that the transition from state J to state J + 1 occurs with rate N λ, while the transition from state J to state J -1 occurs with rate µJ, if J ≤ N , and with rate µN , if J > N . We consider now the scaled process J N (t)/N , whose transition rates from x to x + l/n can be expressed as N β l (x) where the functions β l (x) do not depend on N . In particular

β 1 (x) = λ, β -1 (x) = µx for x ≤ 1, β -1 (x) = µ for x > 1 and β l (x) = 0 otherwise. The rate of changes of J N (t)/N is then f (x) = 1 N (β 1 (x) -β -1 (x)) = λ -µx if x ≤ 1 λ -µ if x > 1
that is a Lipschitz function. This property and the fact that l |l|β l (x) < ∞ guarantee [START_REF] Kurtz | Limit Theorems and Diffusion Approximations for Density Dependent Markov Chains[END_REF] that if J N (0)/N converges to j 0 , J N (t)/N converges to the unique solution 3 of the following equation

dj dt = f (j(t)), j(0) = j 0 . (1)
Observe that j(t) > 1 corresponds to J N (t) > N and then a situation where all the N data centers are working and there are J N (t)-N jobs in the queue. Given that f (x) < 0 for x > 1, Eq. ( 1) shows that j(t) < 1 for large enough t and then after some transient the job queue is asymptotically empty and the number of jobs in the system coincides with the number of busy servers. Moreover, j(t) converges when t diverges: j * j(∞) = λ/µ = ρ. This value is the only accumulation point for the possible trajectories of j(t) and then it is also the stationary probability that a server is busy in the original Markov chain [START_REF] Benaïm | A Class of Mean Field Interaction Models for Computer and Communication Systems[END_REF] (as it is known from the analysis of the M/M/N queue). In a similar way, it is possible to show that if S N (0)/N converges to s 0 , S N (t)/N converges to the solution of the following equation

ds dt = ν S -(ν S + ν C )s(t), s(0) = s 0 , (2)
and when t diverges s(t) converges to s * ν S /(ν S + ν C ), that is the stationary probability that a given datacenter is powered by renewables.

It is clear that we would not have needed fluid models to derive the asymptotic probability that a datacenter is busy or that it is powered by the renewables, but the fluid models allow us to evaluate simply the transient dynamics for the percentage of busy datacenters and of datacenters powered by renewables. Moreover, they are required to characterize the quantity B S (t)/N that is needed to quantify how many datacenters work using the cheap renewable energy.

In Fig. 1 we show the Markov chain transitions affecting B S (t), i.e. the number of datacenters working and powered by renewables, when the number of jobs in the system J(t) is smaller than N . As we observed above, for N large enough J(t) < N holds with probability arbitrarily close to one after some finite time depending on J N (0). For this reason, we can for simplicity assume that the system is in this situation. Observe that the transition indicated in the figure by the dashed line is possible only for specific values of S and B S . If a new job arrives and there are idle datacenters in state S (i.e. B S < S) then the job will be assigned to one of them and B S will increase by one unit. Otherwise B S will stay constant. If we calculate the drift for B S (t)/N when (J(t), S(t), B S (t)) = (j, s, b s ) as done above we obtain that it is equal to

g(j, s, b s ) = λ -(ν S + ν C + µ)b s + ν S j if b s < s, -(ν S + ν C + µ)b s + ν S j if b s = s.
Unfortunately the function g() is not continuous, and then neither Lipschitz.

Nevertheless, [START_REF] Gast | Markov chains with discontinuous drifts have differential inclusion limits[END_REF] shows that when B N S (0)/N converges in probability to b s,0 , B N S (t)/N is related to the solutions of the following differential inclusion

db s dt = G(j, s, b s ) = {g(j, s, b s )} if b s < s, [g(j, s, s), g(j, s, s) + λ] if b s = s. (3)
b s (0) = b s,0

The set-valued function G(j, s, b s ) coincides with g(j, s, b s ) for b s < s, while G(j, s, s) is the interval obtained by the convexification of the accumulation points of g(j, s, b s ) when b s = s. Equation ( 3) admits at least a solution because G() is upper-semicontinuous and Theorem 5 in [START_REF] Gast | Markov chains with discontinuous drifts have differential inclusion limits[END_REF] shows that in such case

inf bs∈D sup 0≤t≤T B S (t) N -b s (t) → P 0,
where D is the set of solutions of Eq. ( 3). This result has practical utility if we can prove that the differential inclusion (3) has a unique solution. A standard sufficient condition for the uniqueness of the solution is the one side Lipschitz condition [START_REF] Kunze | Non-Smooth Dynamical Systems[END_REF], that unfortunately does not hold for G(). We suspect that Eq. ( 3) has a unique solution, but we have not been able to prove it. Nevertheless, we can prove that any possible solution converges to the same value as t diverges. This is enough to draw conclusions about the stationary distribution of our stochastic system. We start observing that for large t j(t) and s(t) are arbitrarily close respectively to the values j * = ρ and s * . It holds

g(j * , s * , s * ) = ν S ρ - ν S + ν C + µ ν S + ν C < 0,
because ρ < 1. If λ + g(j * , s * , s * ) < 0, then all the values of G(j * , s * , b s ) are negative when b s belongs to an opportune interval (s * -ǫ, s * ] and then any possible trajectory of b s (t) will be constrained to the interval [0, s * -ǫ],

where the differential inclusion (3) reduces to a usual differential equation with Lipschitz drift and then it admits a unique solution. This solution converges to (λ + ρν S )/(ν

S + ν C + µ) when t diverges. If λ + g(j * , s * , s * ) > 0, then for b s < s * g(j * , s * , b s ) > λ -(ν S + ν C + µ)s * + ν S j * = λ + g(j * , s * , s * ) > 0,
and any trajectory of b S (t) converges to s * , that is a stable point because in this case 0 ∈ G(j * , s * , s * ). Summarizing, it holds

b * s b s (∞) = s * , if λ + g(j * , s * , s * ) > 0 (λ + ρν S )/(ν S + ν C + µ), otherwise.
By observing that λ + g(j * , s * , s * ) > 0 is equivalent to s * < (λ + ρν S )/(ν S + ν C + µ), and replacing λ = ρµ we can write in a more compact way: Figure 2 shows how the stationary distribution of B N S /N converges to b * s as N increases. The quality of the fluid approximation is different for different values of the load ρ. In particular as far as ρ is far from the critical value for which s * = ρ(ν S + µ)/(ν S + ν C + µ), corresponding to the non differentiability in Eq. ( 4), the approximation is very accurate even for N = 20 datacenters. For the critical load, the federation should include an order of magnitude more datacenters to achieve a good level of approximation. For a given value of N the quality of the approximation improves (/worsen) the larger (/smaller) is the acute angle between the two segments determined by the fluid model, as it happens if ν C increases (/decreases).

b * s = min s * , ρ ν S + µ ν S + ν C + µ . ( 4

Exploiting the model

In this section we show how our simple fluid model can help quantifying the potential advantages of a federation of datacenters and the effect of the different parameters.

We start by discussing Eq. ( 4). The percentage b * s of datacenters working and powered by renewables is obviously limited by the percentage s * of datacenters powered by renewables, and by the percentage ρ of datacenters working, then b * s ≤ min{s * , ρ}. These two regimes appear also in Eq. ( 4) and we refer to them as the renewables-limited regime and the load-limited regime. In particular, Eq. ( 4) shows how close the dispatching algorithm can approach the bound min{s * , ρ} when, as we assumed, the job will be completed by the datacenter that started working on it. The factor (ν S + µ)/(ν S + ν C + µ) multiplying ρ takes into account the fact that a datacenter may change status from S to C (or the other way around) after starting to process a job. These changes limit the utility of job scheduling.

Without the federation every datacenter receives a load ρ and can exploit renewables a fraction s * of the time. Then the percentage of time a datacenter works and is powered by renewables is ρs * , that is smaller than b * s from Eq. (4):

min s * , ρ ν S + µ ν S + ν C + µ > ρs * , because ρ < 1 and (ν S +µ)/(ν S +ν C +µ) > ν S /(ν S +ν C ) = s * .
The difference between the left hand side and the right hand side of the inequality times N quantifies how many additional datacenters work powered by renewables thanks to the federation in comparison to the situation when there is no federation. In what follows we compare the corresponding average energy costs, by normalizing the energy cost per time unit to 1 (/0) when the datacenter is (/is not) powered by renewables. The average energy cost per time unit and per datacenter is then: We focus on the relative cost reduction achieved by the federation in comparison to the uncoordinated case, i.e., on (c nf -c f )/c nf . Fig. 3 shows how the relative cost changes as renewables' dynamics become faster for two different values of the load ρ. We set ν S = ν C so that s * is constant and equal to 0.5. Eq. ( 4) shows that, when ρ and s * are constant, b * s changes only for the effect of the ratio (ν S + ν C )/µ. In other words, it is not important how fast the quantity of renewable energy produced changes, but how much faster it changes than the job completion time. Intuitively, if this ratio is very large, the scheduling is not effective, because a datacenter changes its status S/C many times before completing the job, so that the job takes advantage of renewables' energy on average a fraction s * of the time, independently from the status of the datacenter when the job execution started. Fig. 3 shows indeed that the advantage of the federation converges to 0 as the ratio (ν S + ν C )/µ diverges. This behaviour is common to both the load values considered. When ρ = 0.5, the system is always in the load-limited regime (b * s = ρ(ν S + µ)/(ν S + ν C + µ)) and the advantage of the federation always decreases as the ratio (ν S + ν C )/µ increases. When ρ = 0.65 the system is initially in the renewables-limited regime, so that the relative gain of the federation is limited by the average availability of renewables' energy and the gain is independent on the speed of their dynamics. This situation corresponds to the initial horizontal part of the corresponding curve. As the speed of renewables' dynamics further increases, the scheduling is no more able to effectively follow them and the system enters in the load-limited regime. The relative improvement from the federation in this regime is independent from ρ, so that both curves in Fig. 3 overlap.

c f ρ -b * s ,
Our analysis shows significant reduction of energy costs achievable by the federation of different datacenters, but, until now, we have assumed that the states of the renewables' sources at the different datacenters are independent. This is not true in general. For example, production from PV panels or wind turbines are clearly positively correlated at nearby locations. When energy quantities produced at the datacenters are positively correlated, the improvement from scheduling is reduced. In order to quantify the effect of positive correlation, we consider the following simple model. We assume that the Markov chain determining the state of a renewable source (S or C) is modulated by an underlying Markov chain that is common to all the different sources. In particular, as a toy-example, we consider a Markov chain with two states G and B. The transition rates ν S and ν C of each renewable source depend now on the particular state of the modulating Markov chain and we denote them ν S,G , ν C,G , ν S,B and ν C,B . We consider that

s * G = ν S,G ν S,G + ν C,G > ν S,B ν S,B + ν C,B = s * B
and then states G and B correspond respectively to good and bad weather (at least for the purpose of renewable energy production). It is possible to extend simply our previous analysis, if we assume that the dynamics of the modulating Markov chain are much slower than those of the modulated chain and of job execution (i.e., max{ν G , ν B } << min{ν S,G , ν C,G , ν S,B , ν C,B , µ}).

In such case, the average percentage of datacenters working and powered by datacenters can be obtained through a weighted sum of what would happen without modulation as follows

b * s ≈ π G b * s,G + π B b * s,B ,
where

π G = ν G /(ν G + ν B ), π B = ν B /(ν G + ν B ) and b * s,G (resp. b * s,B
) is calculated from Eq. (4) replacing the rates ν S and ν C by ν S,G and ν C,G (resp. ν S,B and ν C,B ). As we anticipated, the modulating Markov chain correlates the state of the renewable sources. We can quantify this effect by using the correlation coefficient η defined as

η = (s * G ) 2 π G + (s * B ) 2 π B -(s * G π G + s * B π B ) 2 (s * G π G + s * B π B )(1 -s * G π G -s * B π B )
.

As a sanity check, we observe that if ν S,G = ν S,B and ν C,G = ν C,B (i.e., the modulating Markov Chain has no effect on the renewables' state evolution), then η = 0. If instead we have that s * G = 1 and s * B = 0, then η = 1, because all the datacenters are in state S when the modulating Markov chain is in state G and in state C when the modulating Markov chain is in state B.

(ν G = ν B = 0.00001, ν S,G + ν C,G = ν S,B + ν C,B = 0.002, ν S,G = ν C,B , µ = 1).
Figure 4 shows the relative cost reduction due to the federation versus the correlation η. In the specific setting considered, the average percentage of time renewables can power datacenters is constant: s * = π G s * G +π B s * B = 0.5. Then, as the correlation increases s * G increases and s * B decreases of the same amount. As expected, the benefit from the federation is maximum when renewable sources evolve independently (η = 0) and null when at any time they are all in the same state (η = 1). The benefit is non-increasing in η, but, depending on the load ρ, there is a more or less wide range of correlation values for which the benefit does not depend on η. In order to justify this result, we write the specific expression of b * s neglecting for simplicity the rates ν S,G , ν C,G , ν S,B and ν C,B when summed to µ, that is much larger. Under this approximation, the setting ρ = 0.5 = s * corresponds to the case when the system is at the boundaries between the two regimes for η = 0. When the correlation increases, the system is i) in the load-limited regime in good weather (state G) with a value b * s,G almost constant and equal to ρ and ii) in the renewables-limited regime in bad weather (state B) with a value b * s,B decreasing in η. As a consequence the corresponding curve is decreasing. When ρ > s * = 0.5, the system is the renewables-limited regime in both states G and B when η = 0, and then b * s,G = s * G and b * s,B = s * B . As η increases, the increase of s * G is exactly compensated by the decrease of s * B so that the system exhibits the same relative improvement until η is so large that the system enters in the renewables-limited regime when in bad weather and then the improvement decreases again. Finally, when ρ < s * = 0.5, the system is initially in the load-limited regime in both states, and then b * s,G = b * s,B = ρ, independent on η. Again, the improvement does not depend on η until η becomes so large that the system enters in the renewables-limited regime when in bad weather.

As we have shown, our simple fluid model reveals the existence of two different regimes and helps to understand and quantify their non-trivial interaction as the parameters change.

Conclusions

The paper proposes a model of geographical load balancing strategies for a collection of federated (micro) datacenters powered by renewable energy sources. In our strategy the scheduler uses a selection criterion that prioritizes datacenters where renewable energy is currently produced. For this kind of system we use mean field techniques to derive a simple approximate model that allows us to derive several performance measures. First, asymptotic convergence is proven and the quality of the approximation for finite size systems is evaluated through an ad-hoc simulator. Then, we use the simple fluid model to quantify the effect of the different system parameters and to understand the different tradeoffs.

Figure 1 :

 1 Figure 1: Transitions that bring to a change in B S for J ≤ N .

Figure 2 :

 2 Figure 2: Percentage of datacenters powered by renewables : Fluid Model (b * s ) vs Simulation Averages (ν S = ν C = 0.01 µ = 1, N = 20, 100, 500).

  with the federation c nf ρ -ρs * , without the federation

Figure 3 :

 3 Figure 3: Cost reduction due to the federation vs speed of renewables' dynamics (s * = 0.5).

Figure 4 :

 4 Figure 4: Cost reduction due to the federation vs renewables' correlation(ν G = ν B = 0.00001, ν S,G + ν C,G = ν S,B + ν C,B = 0.002, ν S,G = ν C,B , µ = 1).

(

  National Research Agency, ANR), and in part by the Università Italo-Francese, call Galileo 2015-2016, ref. G15-133.

While we need an underlying Markovian process to correctly derive the asymptotic fluid model, empirical results show that the fluid model does not heavily depend on many of these assumptions.

In what follows, convergence of random variables is always "in probability." We omit to repeat it at each time.

Continuity of the right hand side guarantees the existence of the solution and Lipschitz property guarantees uniqueness.

Acknowledgements

This work was supported in part by the "Investments for the Future" Program reference #ANR-11-LABX-0031-01, funded by the French Government


Bernardo Amorim

Daniel Figueiredo
email: daniel@land.ufrj.br

Giulio Iacobelli
email: giulio]@land.ufrj.br

Giovanni Neglia
email: giovanni.neglia@inria.fr

Growing Networks Through Random Walks Without Restarts *

Network growth and evolution is a fundamental theme that has puzzled scientists for the past decades. A number of models have been proposed to capture important properties of real networks. In an attempt to better describe reality, more recent growth models embody local rules of attachment, however they still require a primitive to randomly select an existing network node and then some kind of global knowledge about the network (at least the set of nodes and how to reach them). We propose a purely local network growth model that makes no use of global sampling across the nodes. The model is based on a continuously moving random walk that after s steps connects a new node to its current location, but never restarts. Through extensive simulations and theoretical arguments, we analyze the behavior of the model finding a fundamental dependency on the parity of s, where networks with either exponential or a conditional power law degree distribution can emerge. As s increases parity dependency diminishes and the model recovers the degree distribution of Barabási-Albert preferential attachment model. The proposed purely local model indicates that networks can grow to exhibit interesting properties even in the absence of any global rule, such as global node sampling.

Introduction

The growth and evolution of networks is a fundamental problem in Network Science specially in the light that networks are constantly changing over time. Explaining how and why different real networks grow and evolve the way they do has kept researchers busy for the past decades. Not surprising, various mathematical models for network growth and evolution have been proposed in the literature, either ad-hoc models tailored to specific domains, or general models aiming to capture general principles. A celebrated general network growth model is the Barabási-Albert (BA) model [START_REF] Barabási | Emergence of scaling in random networks[END_REF] which embodies the principle of preferential attachment found in various real networks.

A recognized drawback of most proposed network growth and evolution models is the assumption of global information about the network [START_REF] Blanchard | Small world graphs by iterated local edge formation[END_REF][START_REF] Colman | Local rewiring rules for evolving complex networks[END_REF][START_REF] Gabel | Sublinear but never superlinear preferential attachment by local network growth[END_REF][START_REF] Li | Emergence of global preferential attachment from local interaction[END_REF]. For example, the BA model has a primitive to randomly select a node from the existing network according to the degree distribution. To relax such assumption, models that attach new nodes and edges to the existing network using local attachment rules, such as the Random Walk Model [START_REF] Saramäki | Scale-free networks generated by random walkers[END_REF][START_REF] Vázquez | Growing network with local rules: Preferential attachment, clustering hierarchy, and degree correlations[END_REF], have been proposed. Clearly, random walks require knowledge of the current node degree and its neighbors, a much more localized information. Moreover, it seems more reasonable that new nodes connect to nearby nodes (through some local process) rather than selecting new neighbors from the entire population (through some global process). However, the Random Walk Model studied in [START_REF] Saramäki | Scale-free networks generated by random walkers[END_REF][START_REF] Vázquez | Growing network with local rules: Preferential attachment, clustering hierarchy, and degree correlations[END_REF] and others [START_REF] Gabel | Sublinear but never superlinear preferential attachment by local network growth[END_REF][START_REF] Li | Emergence of global preferential attachment from local interaction[END_REF] still require a primitive to randomly select a node from the network (for the purpose of restarting the walker, for example) and are thus not purely local, because they need to know the number and the identity of all network nodes as well as a way to reach them. Such models have local attachment rules, but global "entry point" selection. More recently, models that have no global primitives have started to be explored [START_REF] Ikeda | Network formation determined by the diffusion process of random walkers[END_REF][START_REF] Ikeda | Network formed by movements of random walkers on a bethe lattice[END_REF]. A drawback of these other models is that they rely on an initial network already containing all nodes such as a lattice or a regular tree, that is then modified according to local rules, and thus are technically not growth models.

In this work we propose and explore a network growth model that is purely local, requiring no global selection over the nodes or any initial network. The model works as follows: 0. Start a network with a single node with a self-loop and place a random walk on this node.

1. Let the random walk take exactly s steps.

2. Connect a new node to the node where the walker resides.

3. Stop if the number of nodes in the network is n, otherwise go to Step 1.

Intuitively, the random walk moves around continuously and after every s steps a new node is added and connected to its current location. The new node immediately becomes part of the network and the walker sees no difference between it and any other node. Note that the model has two parameters s and n and grows an undirected tree (apart from the self-loop at the initial node) since every new node starts with degree one. Moreover, the random walk is uniform on the neighbors and is never restarted, thus its name NRRW (No Restart Random Walk) model. Figure 1 illustrates a sample path of network growth with s = 1 and s = 2. Can such purely local model give rise to interesting network structures such as networks that exhibit a power-law degree distribution? Interestingly, we uncover various non-trivial features of this model such as the fundamental dependency on the parity and magnitude of s and its relationship to the degree distribution. If s is odd and small we find that networks generated by the model tend to have very short-tailed degree distribution and very long distances. On the contrary, if s is even and small, networks exhibit a special kind of power law degree distribution (to be formalized later) and very short distances! As s increases the effect of parity decreases and networks exhibit a heavy-tailed degree distribution. Interestingly, with s large enough, the observed degree distribution follows a power law with exponent identical to the network generated by the BA model, recovering the effect of preferential attachment. We also rigorously prove that for s = 1 the random walk is transient and the degree of every node is bounded from above by a geometric distribution. Other interesting features will be highlighted in what follows.

The model here proposed is very related to the Random Walk Model [START_REF] Saramäki | Scale-free networks generated by random walkers[END_REF] which also allows a random walk to take s steps before connecting a new node. The key difference is that in [START_REF] Saramäki | Scale-free networks generated by random walkers[END_REF], after a new node is attached to the network, the random walk is restarted uniformly at random across all existing nodes in the network. Our random walk never restarts, and is therefore a purely local model. Interestingly, the authors of [START_REF] Saramäki | Scale-free networks generated by random walkers[END_REF] show (through simulations and approximations) that their model is closely related to the BA model and yields a power law degree distribution independently of s. However, recently this finding has been questioned and for s = 1 it was mathematically proven that this is not the case [START_REF] Cannings | Random walk attachment graphs[END_REF]. Our model and findings contributes to this debate and possibly sheds light on how both results could be reconciled (more on this on Section 6).

The remainder of this paper is organized as follows. Section 2 discusses the model and its intuitive behavior, as well as the connection with prior works. Section 3 presents the evolution of node degree induced by the model. Section 4 analyzes the depth of the tree generated by the model. Section 5 presents our theoretical findings for the case s = 1, showing the transient nature of the model in this case. Finally, we summarize our findings and present a brief discussion in Section 6.

Network Growth Model

As presented in Section 1, NRRW (No Restart Random Walk) model can be interpreted as a simple random walk that attaches a new node to its current location every s steps. Similar proposed random walk models for growing networks assume that the random walk restarts either after connecting a new node or adding some number of edges to the new node [START_REF] Saramäki | Scale-free networks generated by random walkers[END_REF][START_REF] Vázquez | Growing network with local rules: Preferential attachment, clustering hierarchy, and degree correlations[END_REF]. A restart consists of randomly selecting a node from the existing network (usually uniformly) and placing the random walk on that node. Despite the similarities, the lack of restarts makes NRRW fundamentally different from models with restart. In particular, the restart significantly reduces the correlation between consecutive node additions since it is very unlikely that the random walk will visit the previous new node when walking to add a new node. Intuitively, the random walk loses memory at every restart. Moreover, restarts have the drawback of assuming knowledge of all network nodes and random access to any such node, and is thus not a purely local growth model.

What is intuitively the behavior of NRRW? In a sense, when s is large the random walk will have little memory between node additions. However, this behavior is different from restarts since the random walk will not find itself on a node chosen uniformly at random but on a node chosen randomly proportional to its degree. 1 Thus, when s is large the NRRW seems similar to the BA model since new nodes connect to random nodes chosen proportionally to their degree. However, since s is fixed and the network grows, will NRRW indeed exhibit a behavior similar to BA model when s << n and then s will finally become small in comparison to the network size?

What about small values for s? Intuitively, the random walk will frequently stumble over the newly created nodes. Interestingly, this local behavior depends fundamentally on the parity of s. If s = 1 then the random walk can always walk to the newly created node and add a new node to it. Such behavior is just not possible if s = 2 and the walker is not on the root. This qualitative difference is not limited to s = 1 and s = 2. When s is odd the walker can always land to the most recently added node after s steps and then add a new node. For s even, this is impossible unless the walker does not traverse the loop at the initial node.

The above observation justifies why in the NRRW model we consider a single node with a self-loop as a starting point. If this was not the case, for any s even the random walk would only add nodes to the original node, trivially constructing a star since it can never step on a newly created node. The loop allows a change of "parity" with respect to the levels of the tree where new nodes can be added. In fact if the random walk is at level k of the tree and s is even, the random walk can only add new edges at the levels k + 2h for h = -k 2 , . . . , -1, 0, 1, 2, . . . until it does not traverse the loop. For s odd this is not necessary as the walker can step on a newly created node to be able to add to nodes in any level of the tree without returning to the root. Thus, yet another fundamental difference between s even and odd.

Will these differences between small and large s and even and odd s manifest themselves in structural properties of the trees generated by the model?

In particular, will the degree distribution fundamentally depend on s? In what follows we explore the degrees and other properties of the trees generated by the model showing in fact, that s plays a key role.

Simulations

In order to study the model we designed and implemented an efficient simulator (in C++) for the NRRW model which has as parameters s, n and r, with r denoting the number of independent runs. For each run, we start with a single node with a self loop, move the random walk s steps, connect a new node to its current location, and repeat. We collect statistics for the various properties merging the results across the r simulation runs, such as degree distribution (fraction of nodes with degree k across all runs). The worst case time complexity of a simulation run is O(sn log n) but the amortized time complexity is O(sn), as we use a growing vector to represent the neighbors of a node that doubles its capacity when needed. Thus, a walker step requires Θ(1) time and a node addition takes O(1) amortized time.

Degree Behavior

In this section we study the degree distribution of NRRW through extensive simulations illustrating its behavior and dependencies. Figure 2 shows the Complementary Cumulative Distribution Function (CCDF) of nodes' for various values for s. Surprisingly, when s is small (between 1 and 8) the respective degree distributions are fundamentally different, exhibiting a kind of power law for s even and an exponential tail for s odd. Note that when s = 1 we do not observe nodes with degree larger than 40 while for s = 2 a nonnegligible fraction of nodes have degree greater than 10 5 . We also observe opposite trends in the degree distribution as s increases. For s odd, increasing s yields a distribution with heavier tails, while for s even increasing s yields a distribution with a lighter tails. As s increases into a medium range (between 15 and 64) the trends continue and the two distributions approach each other. For even larger s (between 127 and 256) the degree distributions become very close, being almost indistinguishable. Interestingly, with large s the degree follows a power law distribution, suggesting that the effect of even s value dominates the dynamics. Moreover, for large s the CCDF exhibits a power law with exponent approximately -2 as it is also the case for the BA model which is based on linear preferential attachment. 2 This supports our initial intuition that when s is large, the random walk samples nodes (adding a new node and connecting to it) with probability proportional to their degree, behaving similarly to the BA model.

Figure 3 shows the degree distribution for s = 2 but over different values for n. Interestingly, note that independent of n the degree distribution exhibits the same power law exponent. However, as n increases the fraction of nodes greater than k becomes smaller for any fixed k > 0 (with the exception of the cut-off regime which occurs when k is near n). This implies that the fraction of nodes with k = 1, the minimum degree, is increasing with n. This is clear by observing d = 1 (leftmost point in x-axis) and noting that the fraction of nodes with degree greater than 1 is decreasing with n. Note that such behavior does not occur for s = 1 which maintains its degree distribution as n increases (the dots for different n values are barely distinguishable in plot).

If for s = 2 the fraction of nodes with degree 1 increases and converges to 1 as n goes to infinity, then we cannot claim that the degree distribution follows a power law. However, we can consider the degree distribution of the nodes that do not have degree 1. In particular, the conditional degree distribution, conditioned on D > 1, is shown in Figure 3. Note that the conditional degree distribution does not show dependance on n and moreover seems to follow a power law. This finding is quite interesting since the fraction of nodes with degree 1 can converge to 1 (as n → ∞) while the remainder of nodes can still follow a power law. This may shed new light on the contrasting results in [START_REF] Saramäki | Scale-free networks generated by random walkers[END_REF] and [START_REF] Cannings | Random walk attachment graphs[END_REF]. We return to this discussion in Section 6.

Figure 3 also shows the fraction of nodes with degree greater than 1, P(D > 1), as a function of n for different even s values (for s odd, it does not depend significantly on n as shown in the top right part of Figure 3 for s = 1). Note that for small s the fraction goes to zero reasonably fast (and thus, the fraction of nodes with degree 1 goes to one). As s increases the rate at which P(D > 1) decreases also decreases. Note that for large s (128 or 256) this decrease is barely noticeable, despite being present. Interestingly, as s odd increases, P(D > 1) decreases but without showing any dependency on n. When s = 257, P(D > 1) approaches the value shown in Figure 3 for s = 256 (result not shown due to space constraints). . Empirical degree CCDF conditioned on the degree being greater than 1; fraction of nodes with degree greater than 1 (bottom plots).

As shown, the NRRW model has a very particular behavior with respect to the degree distribution. In Section 6 we provide a further discussion with a few conjectures for its asymptotic behavior with n.

Level Behavior

We now investigate the level of the nodes on the trees generated by the NRRW model. 3 As we have shown above, the model dynamics has a fundamental dependence on the parity of s, specially when s is small. Indeed, this dependence also manifests itself on the level of the nodes. Figure 4 shows the level distribution (fraction of nodes at level larger than ) for a few small values of s separated into odd and even, respectively. The level distribution for s even decreases very fast. Note that although n = 10 6 , when s = 2, 90% of nodes are at level 4 or less and no node is at level 7 or higher. As s even increases the level distribution decreases relatively slower, with 90% of nodes found at level greater than 4 when s = 8. Still, no node is found at level greater than 10. The behavior is completely different for s odd, and the level distribution seems to be uniform (straight line on a linear-linear plot). For s = 1 the distribution has the heaviest tail with about 4 nodes per level, giving rise to 2.5 • 10 5 different levels. For s = 7 there are about 40 nodes per level, giving rise to 2.5 • 10 4 different levels. Interestingly, as s even increases the level distribution becomes heavier while as s odd increases the level distribution becomes lighter. Figure 4 also shows the level distribution for large s. Indeed, as s increases the level distributions for s even and s odd become more similar and the dependency on the parity diminishes. This behavior is similar to what observed for the degree distribution, illustrated in Section 3.

Note that from the level distribution we can infer the kind of trees that NRRW generates. When s is small and even, the trees generated are "fat and short", with most nodes near the root and a few with very large degrees. When s is small and odd, the trees are "thin and long" with few nodes spread across many levels and no node with large degree. As s increases, the two kind of trees move in each other's direction, becoming more and more similar. 5 Theoretical findings for s = 1

The numerical simulations with s odd and in particular with s = 1 suggest that trees grow in depth as the number of nodes increases. In particular, the growth in depth seems linear on the number of nodes. This is an indication that the random walk is continuously pushing the tree to lower depths just never to return to its origins. In a nutshell, the random walk is transient and visits each node in the tree only a relative small number of times, with high probability. The following Theorem rigorously formalizes this intuition.

Theorem 1. In the NRRW model with s = 1, the number of visits to a node is stochastically dominated by 1 plus a geometric random variable with support on Z >0 .

Proof. We consider here that the initial network consists of a single node with no self-loop. This simplifies the notation and does not compromise the main argument. Let r denote the initial node of the growing network hereafter referred to as the root of the tree (at any step the growing network is a tree) where the walker resides at time zero. Note that r is the only node at level zero. Let X n be the level (i.e. the distance from the root) of the node visited by the NRRW at step n. We call the process {X n , n ∈ Z ≥0 } the level process. Note that the random walk visits r the same number of times that the level process visits level zero. At step n > 0 the NRRW is in a node v n with at least two edges: the one the NRRW has arrived from and the new one added as a consequence of the NRRW's arrival. Let d n ≥ 2 denote the degree of node v n . If v n = r, the NRRW jumps from v n to a node with larger level with probability d n -1

d n ≥ 1 2
and with the complementary probability 1 d n ≤ 1 2 to a node with smaller level. If v n = r, then the level can obviously only increase. Note that, due to the fact that degrees keep changing because of the arrival of new edges, the level process is non-homogeneous (both in time and in space).

We now study the evolution of the level process every two steps, i.e. we consider the process Y n X 2n . Given that the network is a tree and X 0 = 0, the two-step level process can be seen as a non-homogeneous reflecting 'lazy' random walk on 2Z ≥0 = {0, 2, 4, . . .}. We denote by p k,h (n) the probability that the level at step n + 1 is h conditioned on the fact that it is k at step n.

Although the notation hides it, we observe that the probabilities p k,h (n) depend on the whole history of the NRRW until step n. The reason to consider the two-step level process is that we can get bounds on the transition probabilities p k,h (n) that allow a simple comparison with a (biased) homogeneous random walk. The bounds derived above for X n lead immediately to conclude that p k,k+2 (n) ≥ 1 for k ≥ 2, but we can get a tighter bound for p k,k-2 (n). If the NRRW is at level k, all the nodes on the path between its current position and the root r have degree at least 2. If it then moves to node v at level k -1, a new edge is attached to v, whose degree is now at least 3. The probability to move from v further closer to the root to a node with level k -2, is then at most

1 3 . It follows then that p k,k-2 (n) ≤ 1 2 1 3 = 1 6 for k ≥ 2.
We consider now a homogeneous biased lazy random walk (Y * n ) n≥0 on 2Z ≥0 starting from 0 with transition probabilities p *

k,k+2 = 1 4 for all k ∈ 2Z ≥0 and p * k,k-2 = 1 6 for k ∈ 2Z ≥0 and k = 0. We show that if (Y * n ) n≥0 also starts in 0 (Y * 0 = 0), it is stochastically dominated by (Y n ) n≥0 .
We prove it by coupling the two processes as follows. Let (ω n ) n≥0 be a sequence of independent uniform random variables over [0, 1]. We use them to generate sample paths for both processes (Y n ) n≥0 and (Y * n ) n≥0 as follows:

Y n+1 =      Y n -2, if ω n ∈ [0, p k,k-2 (n)) Y n + 2, if ω n ∈ [1 -p k,k+2 (n), 1] Y n otherwise Y * n+1 =      Y * n -2, if ω n ∈ [0, p * k,k-2 ) Y * n + 2, if ω n ∈ [1 -p * k,k+2 , 1] Y * n otherwise where p k,k-2 (n) and p * k,k-2 are 0 if k = 2.
= h + 1 ≤ k + 1 = Y n+1 . If h ≥ 2, then p * h,h-2 = 1 6 ≥ p k,k-2 (n) and if Y n decreases then Y * n must also decrease (Y * n+1 = h -1 ≤ k -1 = Y n+1
) n≥0 is transient since p * k,k+2 = 1/4 > p * k,k-2 = 1/6.
Thus, the probability of the first return time to level 0 is f 0 < 1. By the strong Markov property, the number of visits to level 0 is geometrically distributed on the set Z >0 with parameter equal to 1f 0 . Since a visit to level zero in (X n ) n≥0 (one level process) implies a visit to level zero in (Y n ) n≥0 (two level process), then it follows that the number of visits of (X n ) n≥0 to level zero is bounded by a geometric random variable and then even more so by 1 plus the same geometric random variable. Now let us consider any node v in the growing network. If the NRWW never visits v, then the degree of b is 1 and the thesis follows immediately. Otherwise, let consider the first time the NRWW visits v to be time t = 0 and let consider v to be the root of the current tree. We can retrace the same reasoning and conclude that the number of visits to v for t > 0 is bounded by a geometric random variable on Z >0 with parameter equal to 1f 0 . Then the total number of visits to v is bounded by 1 plus such random variable. This concludes the proof.

Corollary 1. In the NRRW model with s = 1, the degree distribution of any node is bounded by a geometric distribution. This follows since the degree of every node equals the number of visits of the random walk to the node plus 1 (the plus 1 accounts for the fact that any node joining the network, although not yet visited by the walker, has degree 1).

Discussion and Conclusion

As we have shown, the NRRW model exhibits interesting features that fundamentally impact the networks it generates. For s = 1 the random walk is transient and node degree is bounded by a geometric distribution (Theorem 1). For s = 2, the fraction of nodes with degree 1 seems to converge to 1 as n → ∞. However, the conditional degree distribution seems to follow a power law. Can such results be made mathematically rigorous? Other interesting questions also emerge from our analysis of the NRRW model. In particular, our numerical simulations seem to indicate that for any s even, the fraction of nodes with degree 1 will converge to 1 as n → ∞. On the other hand, our simulations also indicate that this is not the case for any s odd. So will there be a fundamental difference between a fixed but arbitrarily large even and odd s? It is hard to imagine that s = 2 10 and s = 2 10 -1 would have fundamentally different behavior, since in both cases the random walk moves quite a lot before adding a new node. Of course, any fixed s will be small as n → ∞. Thus, we make the following conjecture: Conjecture 1. For any fixed s even, the fraction of nodes with degree one converges to 1, as n → ∞. For any fixed s odd, the fraction of nodes with degree one converges to a number strictly less than 1, as n → ∞.

If true, such conjecture would imply that the degree distributions are also never identical, for any fixed s even or odd. However, our numerical results indicate that the conditional degree distribution (conditioned on degree being greater than one) for s even, seems to converge to a power law as n → ∞. On the other hand, for s = 1 we have proved that the random walk is transient and degree distribution is bounded by a geometric distribution (Theorem 1). Can fixed odd s values really generate power laws? If this is the case, then there would be a phase transition on s, from inducing a network with degree distribution with an exponential tail (s = 1) to a power law tail. Despite the numerical results indicating the heavy tail degrees for s = {127, 255}, we make the following conjecture: Conjecture 2. For any fixed s even, the conditional degree distribution is bounded from below by a power law, as n → ∞. For any fixed s odd, degree distribution is bounded from above by an exponential, as n → ∞.

Such conjectures consider that n diverges. In practice n must be finite when generating a network with NRRW model. Thus, for a fixed n, the differences induced by an even or odd s may diminish as s increases. In particular, the degree distribution generated by even and odd s values may become arbitrarily close as s increases, as we have observed in numerical simulations for a fixed n (Figure 2).

Last, we return to the recent dispute if the Random Walk model with restarts generates a power law degree distribution, independently of s [START_REF] Cannings | Random walk attachment graphs[END_REF][START_REF] Saramäki | Scale-free networks generated by random walkers[END_REF]. It has been mathematically proved that when s = 1 the fraction of nodes with degree one converges to 1, as n → ∞ [START_REF] Cannings | Random walk attachment graphs[END_REF]. At the same time, simulation results suggest that the degree distribution follows a kind of power law [START_REF] Saramäki | Scale-free networks generated by random walkers[END_REF]. We can attempt to reconcile such findings by leveraging our own findings on NRRW model. When s = 1 the new node is connected to a given existing node u if i) a neighbor of u is selected at the restart and then ii) the random walk moves to u. A node whose neighbors are all leaves would be selected with a probability proportional to its degree. Now it has been shown that when n → ∞ the fraction of nodes that are leaves converges to 1, then most of the neighbors of a non-leaf node are leaves and this node is essentially selected proportionally to its degree, similarly to the BA model embodying preferential attachment. Thus, the conditional degree distribution, leaving out degree 1 nodes, will follow a power law distribution with the same exponent as in the BA model. In some sense, this reconciles the findings of the two prior works [START_REF] Cannings | Random walk attachment graphs[END_REF][START_REF] Saramäki | Scale-free networks generated by random walkers[END_REF].

To conclude, as exemplified above, a fundamental understanding of NRRW model adds to our understanding of purely local network growth models. In particular, besides requiring a less strict assumption to operate, models that do not rely on any global primitive can also generate networks with rich and diverse structural properties.

Acknowledgement

Research conducted within the context of the THANES Associate Team, jointly supported by Inria (France) and FAPERJ (Brazil). This work has also been partially funded through research grants from CNPq and CAPES (Brazil).

Figure 1 :

 1 Figure 1: Examples of sample path for network growth for s = 1 (top) and s = 2 (bottom). The red square denotes the walker position. The snapshots represent the growing network just after the new node is connected.

Figure 2 :

 2 Figure 2: Empirical degree Complementary Cumulative Distribution Function (CCDF) for various values of s in loglog scale (n = 10 6 , r = 10 3 ).

  of non-leaves for each N; s even

Figure 3 :

 3 Figure 3: Empirical degree CCDF for various values of n in loglog scale (top plots). Empirical degree CCDF conditioned on the degree being greater than 1; fraction of nodes with degree greater than 1 (bottom plots).

Figure 4 :

 4 Figure 4: Empirical CCDF of the node level for different s values (n = 10 6 , r = 10 3 ).

  level k ≥ 0 and p k,k-2 (n)

  We start observing that if Y n and Y * n have the same value k, then every time Y n increases also Y * n increases because p * k,k+2 = 1 4 ≤ p k,k-2 (n). On the contrary if Y * n decreases (as it can happen only for k ≥ 2), then Y n may decrease or not because p k,k-2 (n) ≤ 1 6 = p * k,k-2 . It follows that if Y n and Y * n are at the same level, then Y * n+1 ≤ Y n+1 . We are going to prove by induction on n that Y * n+1 ≤ Y n+1 for every n. With a slight abuse of terminology we say that Y n increases (resp. decreases) if Y n+1 > Y n (resp. Y n+1 < Y n ). We start observing that indeed Y * 0 ≤ Y 0 , because both processes start in 0. Let us assume that Y * n = h ≤ k = Y n . For all values of h, every time Y * n increases also Y n increases because p * k,k+2 =

	1 4 ≤ p k,k+2 (n) and then Y * n+1

  ). It follows that for h ≥ 2 then Y * n+1 ≤ Y n+1 . The only case when Y n may decrease without Y * n decreasing is when h = 0 and k = 0, but in this case Y * n+1 = 0 and Y n+1 ≥ 0. This proves that Y * n+1 ≤ Y n+1 for every n. Given that Y * n ≤ Y n and both processes start at level zero, the number of visits of (Y n ) n≥0 to level zero is bounded by the number of visits of (Y * n ) n≥0 to level zero. The homogeneous biased lazy random walk (Y * n

Recall that the steady state distribution of a random walk on a fixed network is given by d i / ∑ j d j , where d i is the degree of node i.

Recall that if D follows a power law distribution, then P(D = k) ∼ k -α where α > 1 is the power law exponent, and it follows that P(D ≥ k) ∼ k -(α-1) . Thus, the CCDF has an exponent that is one unit less than the PDF.

Recall that the level of node on a tree is given by its distance to the root, and thus the root is at level zero.