Ido Cohen

Guy Gilboa

Shape Preserving Flows and the p-Laplacian Spectra

Introduction and main results

A main principle in image and signal processing is to represent the data in a meaningful manner, which allows effective processing for the type of signals and application at hand. Usually, a compact and sparse representation is desired. In filtering applications, we would often prefer a representation which well distinguishes between signal and noise. Linear transforms were extensively employed for this purpose. For example, Fourier transform represents harmonic functions in the time domain by delta functions in the frequency domain.

Fourier developed his theory while investigating the solution of the heat equation (or homogeneous linear diffusion), [START_REF] Fourier | Theorie analytique de la chaleur[END_REF] u t = ∆u, u(t = 0) = f. It is well known that the solution u(t) in an unbounded domain is a convolution of f with a Gaussian kernel of standard deviation σ = √ 2t, for a review see e.g. 1 [START_REF] Weickert | A review of nonlinear diffusion filtering[END_REF]. Moreover, if f is an eigenfunction, namely, admits the linear eigenvalue problem ∆f = λf, λ ∈ R, the variables x and t of the solution are separated, where u(x, t) = e λt • f (x).

(

) 1
We note that the heat equation is a gradient descent process with respect to the Dirichlet energy,

J D (u) = 1 2 |∇u| 2 dx.
The Dirichlet energy and its associated flow smooth sharp transitions. Therefore, non-quadratic energies and non-linear flows are more useful for image processing tasks. A widely used edge preserving functional in image processing is the Total Variation (TV) energy,

J T V (u) = |∇u| dx.
It can be perceived as the L 1 analog of the Dirichlet energy. The associated gradient descent, known as TV-flow [START_REF] Andreu | Minimizing total variation flow[END_REF], is a nonlinear diffusion process of the form,

u t = div ∇u |∇u| , u(t = 0) = f. (2)
It can be shown ( [START_REF] Steidl | On the equivalence of soft wavelet shrinkage, total variation diffusion, total variation regularization, and sides[END_REF][START_REF] Burger | Spectral decompositions using one-homogeneous functionals[END_REF]) that in the discrete 1D case the solution at time t coincides with the minimizer of the ROF model [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF],

arg min

u 1 2 u -f 2 + tJ T V (u).
In an analog manner to the linear case, we can define a nonlinear eigenfunction with respect to the TV subdifferential by, q = λf, -q ∈ ∂J T V (f ),

where ∂J(u) denotes the subdifferential of the functional J at u. Similarly to linear diffusion, if the initial condition admits (3) we obtain a simple analytic solution of the flow, where x and t are separated [START_REF] Andreu | Minimizing total variation flow[END_REF],

u(x, t) = (λt + 1) + • f (x). ( 4
)
The linear decay over time of ( 4) is the key to data representation by the TVtransform [START_REF] Gilboa | A total variation spectral framework for scale and texture analysis[END_REF], where a nonlinear eigenfunction (3) becomes a delta in the TV transform domain.

One can generalize the aforementioned energies to a the p-Dirichlet functional:

J p (u) = 1 p |∇u| p dx, ( 5
)
where p is classically in the range (1, ∞). The gradient descent with respect to [START_REF] Burger | Spectral representations of one-homogeneous functionals[END_REF] is the nonlinear diffusion flow (hereinafter p-flow):

u t = ∆ p u, u(0) = f, (6)
where ∆ p u is the p-Laplace operator defined by,

∆ p u = div |∇u| p-2 ∇u . (7)
We would like to investigate whether there are solutions to [START_REF] Burger | Spectral decompositions using one-homogeneous functionals[END_REF] in which x and t are separated, as in Eqs. ( 1) and (4). To address this issue we first consider a more general nonlinear flow

u t = P (u(t)), u(t = 0) = f, (P)
where P is a general nonlinear bounded operator. We seek conditions where the solution is in the form:

u(x, t) = a(t) • f (x). (ShP)
We denote these type of solutions as shape preserving flows, since there is no spatial change in u and the flow only modifies the contrast. The change of contrast is dictated by a(t). The function a(t) is induced by the operator P and can be arbitrary. In our study we are interested in flows that vanish with time.

We therefore term the function a(t) as the decay profile. We begin by formulating necessary and sufficient conditions for solutions of (P) to be in the form of (ShP). It turns out that the initial condition u(t = 0) = f has to be an eigenfunction of P , i.e. admits

P (f ) = λf, λ ∈ R.
Moreover, we found that at t = 0 the directional derivative of P (in the f direction) is proportional to f . When P is an homogeneous operator we can formulate analytic solutions for the flows. In the case of the p-Laplace operator, the decay profile is:

a(t) = [(2 -p)λt + 1] 1 2-p .
This decay profile holds for any value of p in the range [1, ∞), for any eigenfunction f and any eigenvalue λ. As an immediate consequence, the p-flow has a finite extinction time T for 1 ≤ p < 2,

T = 1 (p -2)λ .
Following these results, a new family of nonlinear transforms is introduced. These transforms connect between eigenfunctions of the homogeneous operator and delta functions in the transform domain. Definitions for spectrum, filtering and data representation are formulated. Preliminary experiments illustrate this framework.

The plan of the paper is as follows: In Section 2 necessary and sufficient conditions are given for (P) to be shape preserving. Then we focus our discussion on homogeneous operators and show they can induce shape preserving flows. For these flows we formulate their decay profile and extinction time. In Section 3 we generalized the spectral one homogeneous representation [START_REF] Burger | Spectral decompositions using one-homogeneous functionals[END_REF] to p-homogeneous spectral representation for 1 < p < 2. In Section 4 some numerical results are presented.

The theory of shape preserving flows

In this section we develop the theoretical basis for the p-spectra. We begin by formulating necessary and sufficient conditions for the flow (P) to be shape preserving. Then we show that operators of any homogeneity address these conditions.

Necessary and sufficient conditions for shape preserving flows

In the following theorem we define precise conditions for shape preservation under the assumption that in some initial time range the flow is infinitely differentiable, with respect to time.

Theorem 1 (Shape preserving flow -necessary and sufficient condition). Given a (non)linear scale-space flow, Eq. (P), which exists for all t ∈ R + . We assume u(t) is infinitely differentiable with respect to t in the range D = [0, T ) ⊆ R + . The solution u(t) is shape preserving (Eq. (ShP)) in the range D iff for any n = 0, 1, .. we obtain

d n dt n P (u(t)) t=0 = c n • f, (8)
where c n ∈ R is some constant.

Proof.

⇒ [We assume the solution is shape preserving (admits Eq. (ShP)) and prove Eq. ( 8)]. Let us note that P (u(t)) is infinitely differentiable in the range D using Eq. (P) and the fact that u(t) is infinitely differentiable. By substituting (P) and (ShP) in Eq. ( 8) we can calculate the nth derivative of P . For any n = 0, 1, .. we get:

d n dt n P (u(t)) t=0 = d n+1 dt n+1 u(t) t=0 = a (n+1) (t) t=0 • f = c n • f (x).
⇐[We assume Eq. ( 8) and prove the solution is in the form of (ShP)]. The solution u(t) can be represented as a Taylor series at any time t in D because it is infinitely differentiable. Without loss of generality, we express u(t) at t = 0:

u(t) = ∞ n=0 1 n! d n dt n u(t) t=0 t n =f + ∞ n=0 1 (n + 1)! d n dt n P (u(t)) t=0 t n+1 =f + ∞ n=0 c n f (n + 1)! t n+1 = 1 + ∞ n=1 α n t n f = a(t)f. (9)
Let us discuss some implications of the above theorem. If the flow is shape preserving then (8) holds for any natural n, in particular for n = 0,

P (f ) = λ • f, (10)
where λ = c 0 . Therefore, f being an eigenfunction of P is a necessary condition for a shape preserving flow. Let us further examine [START_REF] Cohen | Energy dissipating flows for solving nonlinear eigenpair problems[END_REF] for n = 1,

d dt P (u(t)) t=0 = lim dt→0 P (u(t + dt)) -P (u(t)) dt t=0 .
When dt → 0 we can substitute P (u(t + dt) by P (u(t) + dt • u t ),

d dt P (u(t)) t=0 = lim dt→0 P (u(t) + dt • u t ) -P (u(t)) dt t=0 .
Substituting u t by using (P) yields:

d dt P (u(t)) t=0 = lim dt→0 P (u(t) + dt • P (u(t))) -P (u(t)) dt t=0 .
Evaluating the above expression at t = 0 we get,

d dt P (u(t)) t=0 = lim dt→0 P (f + dt • P (f )) -P (f ) dt .
But according to Eq. ( 10) f is an eigenfunction of P , therefore

d dt P (u(t)) t=0 = lim dt→0 P (f + λdt • f ) -P (f ) dt . (11)
One can observe that Eq. ( 11) is the directional derivative of the operator P at point f in the direction f (multiplied by λ). From Theorem 1 it is linearly dependent of f . That conclusion is valid for derivative of any order. It turns out that eigenfunctions play a crucial rule in shape preserving flows.

As to the Taylor series at t = 0 (Eq. ( 9)). This series uses the time derivatives of P (u) of any order at t = 0. These derivatives must be right-hand-sided since t ∈ [0, T ). The Taylor series is valid as u(t) is infinitely differentiable in D. The rest of this section is dedicated to a family of operators for which Theorem 1 holds.

Homogeneous operators

Here we focus on homogeneous operator of any positive degree. We show that if the operator P is homogeneous the solution of (P) is a shape preserving flow if and only if the initial condition u(t = 0) = f is an eigenfunction in the form of [START_REF] Fourier | Theorie analytique de la chaleur[END_REF].

We begin by examining the general scale space flow (P), showing that if the operator P is homogeneous then its time derivatives of any order are homogeneous as well (Lemma 1). We conclude that when the initial condition is an eigenfunction of P the conditions of Theorem 1 hold. Moreover, we obtain an analytic expression for the decay profile of the flow (P) (Theorem 2). Consequently, a simple expression for the extinction time for any eigenfunction is formulated for homogeneous operators of degree in the range (0, 1) (Prop. 2).

The homogeneity of the operator can be written in a general manner by,

P (au) = g(a) • P (u), ∀a ∈ R, (12)
given some function g(a). However, since we should require restrictions such as g(ab) = g(a)g(b) (∀a, b) and that the sign of a is preserved, we assume g(a)

takes the simple polynomial form,

g(a) = a|a| p-2 , (13)
where p -1 is the degree of homogeneity. We use this definition as it coincides with the homogeneity of the p-Laplacian operator, which will be investigated more thoroughly later.

Lemma 1. Given the scale-space (P) and u(t) exists and infinitely differentiable. If P is a g(a) homogeneous operator then

(d n /dt n )P is a homogeneous operator with g(a) • g(a)
a n -homogeneity.

Proof. Let us first show the homogeneity of the operator ∂P (u)/∂u. On one hand we know that:

∂P (au) ∂u = g(a) ∂P (u) ∂u .
On the other hand:

∂P (au) ∂u = a ∂P (au) ∂au .
Thus we get:

∂P (au) ∂(au) = g(a) a ∂P (u) ∂u . (14)
Another identity, which holds for u(t) admitting (P), based on the chain rule, is,

d dt P (u) = ∂P (u) ∂u • u t = ∂P (u) ∂u • P (u). ( 15
)
We can now proceed and prove our claim by induction. According to Eqs. ( 14) and ( 15) we have,

d dt P (au) = ∂P (au) ∂(au) • P (au) =g(a) g(a) a ∂P (u) ∂u • P (u) =g(a) g(a) a • d dt P (u). (16)
Therefore, the homogeneity of the first order time derivative of P is g(a)

• g(a)
a . We note that a is in general a function of time, a(t). According to the assumption of the induction, the homogeneity of d n-1 /dt n-1 P is g(a) • (g(a)/a) n-1 . Let us denote the (n -1)th time derivative of P (u) as Q(u). Then, in a similar manner to Eq. ( 15), we have

d dt Q(u(t)) = ∂Q(u) ∂u • P (u).
Also, based on the assumption of the induction and similarly to Eq. ( 14) we get that the homogeneity of ∂Q/∂u is (g(a)/a) n . We therefore have,

d n dt n P (au) = d dt Q(au) = ∂Q(au) ∂(au) P (au) = g(a) a n • ∂Q(u) ∂(u) P (au) =g(a) g(a) a n • ∂Q(u) ∂(u) P (u) =g(a) g(a) a n • d dt Q(u),
which proves the assertion of the lemma.

As a conclusion of Theorem 1 and Lemma 1 we are able to state the following theorem.

Theorem 2. If P is a p -1 homogeneous operator, where g(a) as in [START_REF] Gilboa | Semi-inner-products for convex functionals and their use in image decomposition[END_REF], and u(t) is infinitely differentiable at t = 0, then:

1. The solution of the scale space Eq. (P) exists for all t ∈ [0, T ) and is shape preserving iff the initial condition u(t = 0) = f is an eigenfunction in the sense [START_REF] Fourier | Theorie analytique de la chaleur[END_REF].

2. The flow has the following analytic solution (for t ∈ [0, T )):

u(t) =a(t) • f, a(0) = 1, (17a)
a(t) = [(2 -p)λt + 1] 1 2-p (17b)
as long as a(t) = 0. Or equivalently, it is the solution of the following ODE,

d dt a(t) = g(a(t)) • λ, a(0) = 1.
Proof.

1. According to Theorem 1 a shape preserving flow and Eq. ( 8) are equivalent. ⇒ [We assume [START_REF] Cohen | Energy dissipating flows for solving nonlinear eigenpair problems[END_REF] and prove u(t = 0) = f is an eigenfunction].

Trivially, by assigning n = 0 in [START_REF] Cohen | Energy dissipating flows for solving nonlinear eigenpair problems[END_REF].

⇐ [We assume the initial condition is an eigenfunction and prove [START_REF] Cohen | Energy dissipating flows for solving nonlinear eigenpair problems[END_REF]].

We prove (8) by induction. The assumption of the induction is trivially valid for n = 0, as f is an eigenfunction.

For any n, we assume (d n /dt n )P (u(t)) t=0 = c n • f and compute the (n + 1)th time derivative of P (u(t)). Let us denote P (n) (u(t)) := (d n /dt n )P (u(t)). P (n) (u(t)) is differentiable since u(t) is infinitely differentiable and using Eq. (P), thus we can directly compute P (n+1) at t = 0 by

P (n+1) (u(t)) t=0 = lim dt→0 P (n) (u(t + dt)) -P (n) (u(t)) dt t=0 = lim dt→0 P (n) (u(t) + P (u(t))dt) -P (n) (u(t)) dt t=0 = lim dt→0 P (n) ((1 + λdt)f ) -P (n) (f ) dt .
Using Lemma 1 we have

P (n+1) (u(t)) t=0 = lim dt→0 g(1 + λdt) • g(1+λdt) 1+λdt n • P (n) (f ) -P (n) (f ) dt = lim dt→0 g(1 + λdt) • g(1+λdt) 1+λdt n -1 dt c n f.
Recall that g(1) = 1 and therefore g(1) g( 1)

1 n = 1
and we have = lim

dt→0 g(1 + λdt) • g(1+λdt) 1+λdt n -g(1) g(1)
1

n dt c n f.
The above can be expressed by the following derivative

= d dτ g(τ ) g(τ ) τ n τ =1 • λc n f = [(n + 1)g (1)
-n] λc n f, and using g (1) = p -1 yields

c n+1 = [(n + 1)(p -1) -n] λc n .
2. Using the expression for c n+1 we can directly compute the decay profile and reach (17b) (see Appendix A). However, a simpler approach is to view the solution as a first order ODE. Since the solution is in the form of (ShP), we can express the first time derivative as,

u t = d dt a(t) • f.
In addition, u(t) admits the scale space equation,

u t =P (a(t) • f ) =g(a(t))λf.
Thus, we get d dt a(t) = g(a)λ.

This is separable equation and using (13) yields,

d dt a(t) = a|a| p-2 λ.
Under the assumption a > 0 we have

a(t) 2-p = (2 -p)λt + C.
Based on the initial condition a(0) = 1 gives

a(t) = [(2 -p)λt + 1] 1 2-p .
Let us examine the conditions for a shape preserving flow to blow up in finite time. From the analytic solution (17b) we can state the following corollary:

Corollary 1 (Flows with finite blowup time). If the conditions of Theorem 2 hold, p > 2 and λ > 0, then the solution blows up in finite time T , where

T = 1 (p -2)λ .
This is an immediate consequence of (17b) for a negative power as the base tends to zero. As an example, the p-Laplacian (7) is a negative semidefinite operator. Therefore, its eigenvalues are negative, where the existence of eigenfunctions was proved in [START_REF] García Azorero | Existence and nonuniqueness for the p-laplacian[END_REF]. Consequently, the negative p-Laplacian induces positive eigenvalues and the inverse p-flow, u t = -∆ p (u), for p > 2, blows up in finite time.

Moreover, we can check when the flow vanishes. We refer to the first time T for which u(T ) = 0 as the extinction time. In a similar manner to Cor. 1 we can state a corollary for the extinction time.

Corollary 2 (Flows with finite extinction time). If the conditions of Theorem 2 hold, p < 2 and λ < 0, then the solution has a finite extinction time

T = 1 (p -2)λ . ( 18
)
In the conditions of Cor. 2, the power in Eq. ( 17b) is positive and the base is zero for t = T . As for any positively homogeneous operator P (0) = 0 we get that u is identically zero for t ≥ T . Following the above discussion, the p-flow vanishes in finite time for 1 < p < 2. In addition, if p → 2 Eq. (17a) becomes a(t) = e λt , and the extinction time is infinite as expected for the heat equation. If p > 2 the solution does not reach a steady state. This coincides with results in [START_REF] Xin | Extinction and positivity of the solutions for a -laplacian equation with absorption on graphs[END_REF].

Another interesting case for the p-flow is when p = 1, where we obtain the TV-flow. For the case p = 1 we get a linear decay profile (as proved in [START_REF] Bellettini | The total variation flow in rn[END_REF]). This is a subgfradient flow and the conditions of Theorems 1 and 2 do not hold. However, we see that [START_REF] Li | Numerical methods for fractional calculus[END_REF] still holds. Based on the decay profile of the T V -flow a new concept of nonlinear transform was introduced in [START_REF] Gilboa | A total variation spectral framework for scale and texture analysis[END_REF]. This transform was defined as a second time derivative of u(t),

φ T V (t) = t • u tt , (19)
for which an eigenfunction becomes a delta in the transform domain. This concept was extended to any absolutely one-homogeneous functionals (which have zero-homogeneous subgradients, p = 1) in [START_REF] Burger | Spectral decompositions using one-homogeneous functionals[END_REF] or to the weighted problem introduced in [START_REF] Bungert | Solution paths of variational regularization methods for inverse problems[END_REF]. This was used in several image-processing applications, e.g. for denoising [START_REF] Moeller | Learning nonlinear spectral filters for color image reconstruction[END_REF], segmentation [START_REF] Zeune | Multiscale segmentation via bregman distances and nonlinear spectral analysis[END_REF] and image fusion [START_REF] Hait | Spectral total-variation local scale signatures for image manipulation and fusion[END_REF].

3 The p-Laplace Spectra

In the last decade there is a rising interest in p-flows for 1 < p < 2 (see e.g. [START_REF] Kuijper | p-laplacian driven image processing[END_REF][START_REF] Wei | A p-laplace equation model for image denoising[END_REF][START_REF] Chen | Variable exponent, linear growth functionals in image restoration[END_REF]). Therefore, it would be beneficial to develop a representation method based on these flows.

Extinction time for the p-flow

In the previous section an expression for the time of extinction of the p-flow was given for eigenfunctions, (1 < p < 2). We address now the issue of finite extinction time for arbitrary initial conditions f . Upper bounds on the extinction time were proposed in several studies. For example, in [START_REF] Vázquez | Smoothing and decay estimates for nonlinear diffusion equations: equations of porous medium type[END_REF] the authors establish an upper bound (for 1 < p < 2) based on Barenblatt solution [START_REF] Barenblatt | On self-similar motions of a compressible fluid in a porous medium[END_REF] . In [START_REF] Dibenedetto | Degenerate parabolic equations[END_REF] (see pages 188,189) the authors provide a bound for the continuous setting in a bounded domain in R N . In Proposition 3.12 in [START_REF] Burger | Spectral decompositions using one-homogeneous functionals[END_REF] an upper bound is derived for p = 1 (T V -flow). In addition, it is worth mentioning that in [START_REF] Vázquez | The dirichlet problem for the fractional p-laplacian evolution equation[END_REF] the behavior of the p-flow near the time of extinction is presented as an open problem. In Appendix B we show a different approach to proving finite extinction time of the p-flow using arguments similar to the ones of [START_REF] Burger | Spectral decompositions using one-homogeneous functionals[END_REF] in a semi-discrete setting.

The p-Transform

In a similar manner to the T V -transform [START_REF] Gilboa | A total variation spectral framework for scale and texture analysis[END_REF] and to the one-homogeneous transform [START_REF] Burger | Spectral decompositions using one-homogeneous functionals[END_REF], we would like to express an eigenfunction as a delta in time by applying a linear operator to the flow. We propose to take the time derivative of order α for the p-flow, where

α = 1 2 -p + 1. ( 20
)
As α is not necessarily a natural number, fractional derivatives are required. There are many different definitions of the fractional calculus operators, see e.g. [START_REF] Li | Numerical methods for fractional calculus[END_REF][START_REF] Ortigueira | Fractional calculus for scientists and engineers[END_REF][START_REF] Luchko | Fractional fourier transform and some of its applications[END_REF][START_REF] Khalil | A new definition of fractional derivative[END_REF]. All definitions coincide with derivatives in the classical sense (of integer orders) but can differ for non-integer orders. We chose to use a well known definition of Sonin and Laurent (see [START_REF] Ortigueira | Fractional calculus for scientists and engineers[END_REF]), which is a generalization of Cauchy integral.

Definition 1 (Fractional derivative). The fractional derivative of order α of the function u at point t is given by,

D (α) {u(t)} := Γ(α + 1) 2πj c u(τ ) (τ -t) 1+α dτ ( 21
)
where c is a closed curve around t, Γ is the extension of the factorial function where Γ(n) = (n -1)! when n is natural.

We can now define the p-transform and the inverse p-transform as follows:

Definition 2 (The p-transform). The p-transform is defined by,

φ(t) = (-1) α-1 t α-1 Γ(α) D α {u(t)}, (22)
where u(t) is the solution of (6), p ∈ (1, 2), α is as in (20) and • is a ceil operator.

Definition 3 (Inverse p-transform).

f = ∞ 0 φ(t)dt. (23)
Let us show the reconstruction property for α ∈ N .

Proposition 1 (Reconstruction). The function f can be reconstructed by the inverse transform (Eq. ( 23)), for α ∈ N .

Proof. In this proof we limit ourselves to natural α = n, i.e. p = 2n-3 n-1 where n ∈ N . In this case the p-transform is:

φ(t) = (-1) n-1 t n-1 (n -1)! u (n) (t), n = 2, 3, ...
Then using integration by parts we have

∞ 0 φ(t)dt =(-1) n-1 ∞ 0 t n-1 (n -1)! u (n) (t)dt =(-1) n-1     t n-1 (n -1)! u (n-1) (t) ∞ 0 =0 - ∞ 0 t n-2 (n -2)! u (n-1) (t)dt     =(-1) n     t n-2 (n -2)! u (n-2) (t) ∞ 0 =0 - ∞ 0 t n-3 (n -3)! u (n-2) (t)dt     =... =(-1) 2n-3 ∞ 0 u (t)dt =(-1) 2n-2 f =f 3.3 The p-Spectra Definition 4 (The p-Spectrum). S(t) = f, φ(t) .
The Parseval-type theorem holds,

S(t)dt = f, φ(t) dt = f, φ(t)dt = f, f = f 2 .
Theorem 3 (Spectral response of an eigenfunction). The p-transform, [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], of an eigenfunction f with eigenvalue λ, when α ∈ N , is:

φ(t) = δ t + α -1 λ • f. (24)
Proof. The solution of ( 17) can be written as,

u(t) =    λ α-1 t + 1 α-1 f t ≤ T 0 t > T.
When α ∈ N we can formulate the α-th time derivative of u(t) as

D (α) {u(t)} = (α -1)! λ α -1 α-1 δ t + α -1 λ • f.
Assigning this in [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] gives,

φ(t) =(-1) α-1 t α-1 (α -1)! (α -1)! λ α -1 α-1 δ t + α -1 λ • f =(-1) α-1 t α-1 λ α -1 α-1 δ t + α -1 λ • f.
As the delta function is non zero for t = -(α -1)/λ, this is equivalent to

φ(t) =(-1) α-1 - α -1 λ α-1 λ α -1 α-1 δ t + α -1 λ • f =δ t + α -1 λ • f.

Filtering

The filtering is identical to the previous work on TV-transform (see e.g. [START_REF] Burger | Spectral representations of one-homogeneous functionals[END_REF]).

The filtering is done in the time domain, as follows,

φ(t) = φ(t) • h(t), (25)
where LP, HP, BP and BS are defined as in [START_REF] Gilboa | A total variation spectral framework for scale and texture analysis[END_REF]. The filtered signal is respectively defined:

f h = ∞ 0 φ(t)dt. (26)

Experiments

In this section we demonstrate numerically some of the theoretical findings presented earlier. The p-flow is implemented by simple explicit method using a standard 3 × 3 spatial stencil. The step size is different for different values of p. Unless stated otherwise, our experiments are for p = 1.5. In this case the order of derivative is α = 3 and the step size is dt = 10 -3 . For p = 1.3 the step size is dt = 5 • 10 -4 .

Decay profile. In Fig. 1 we demonstrate the validity of the analytic solution of the decay profile, Eq. (17b) in one-and two-dimensions. On the left, two eigenfunctions of the p-Laplacian are shown, (with eigenvalues λ = -0.017 (top) and λ = -0.0587 (bottom)). On the right, the corresponding theoretical and experimental decays are plotted for a randomly selected point x 0 . One can observe that the experimental results coincide well with the theoretical predictions.

Transform and Spectrum. We illustrate the application of the p-transform on eigenfunctions. Several eigenfunctions were generated (with different values of p), using the algorithm of [START_REF] Cohen | Energy dissipating flows for solving nonlinear eigenpair problems[END_REF], see Fig. 2. We expect that the spectrum of a numerical eigenfunction will be composed of essentially one dominant peak, approximating a delta response at t = 1 (p-2)λ . In Fig. 3 we show the decay profile for p-flow for the eigenfunction whose eigenvalue λ = -0.216. One can see in the different subfigures the decay profile

a) | f, u(t) | (b) f, d dt u(t) (c) f, d 2 dt 2 u(t) (d) f, d 3 dt 3 u(t)
Figure 3:

(a) f, d 3 dt 3 u(t)
Figure 4: and its time derivatives. As expected, the decay profile is a polynomial of second order, its first order derivative is linear, the second derivative is a step function and the third is an isolated peak. In Fig. 4 we show the third order derivative of the p-flow and again, as α = 3, we get an isolated peak. For p = 1.5, the extinction time is T = -2/λ, following [START_REF] Luchko | Fractional fourier transform and some of its applications[END_REF]. It can be seen that the experimental results coincide with the analytic solution i.e. the extinction times are 9.24 and 0.667, as expected. We would like to note that for the eigenfunction with eigenvalue λ = -3 we obtained numerically an almost exact eigenfunction and therefore the extinction time is predicted with very high accuracy. In general, these results fit the analytic solutions formulated in [START_REF] Li | Numerical methods for fractional calculus[END_REF].

(a) f, u(t) (b) f, d dt u(t) (c) f, d 2 dt 2 u(t) (d) f, D α u(t)
In Fig. 5 we demonstrate the transform for p = 1.3, α = 2 3 7 . Here the transform is based on fractional derivatives. There are several ways to implement fractional derivatives [START_REF] Petráš | Fractional derivatives, fractional integrals, and fractional differential equations in matlab[END_REF]. The derivative operator is a semigroup, therefore, we choose to separate the fractional derivative to the natural and the fractional parts. For the natural part (in our case 2) we apply standard derivation and for the fractional remainder (in our case 3 7 ) we approximate it via Fast Fourier Transform (FFT). The eigenvalue is λ = -0.2876, yielding a theoretical extinction time of T theory = 4.967. The experimental result is very close where we get T = 4.9665.

Filtering. We implement the filtering process as suggested in Eqs. ( 25) and [START_REF] Wei | A p-laplace equation model for image denoising[END_REF]. Again, we use p-flow where p = 1.5, i.e. 1/(2 -p) is natural. In Fig. 6 we add two eigenfunctions with different eigenvalues. When the signal is composed of two independent eigenfunctions (see [START_REF] Gilboa | Semi-inner-products for convex functionals and their use in image decomposition[END_REF]) we can expect to obtain two distinct groups of peaks. And indeed we get that result in Fig. 6d. Filtering the red marked spectrum out of the signal results in recovering the eigenfunction with the lower eigenvalue. Naturally, subtracting this result from the initial image approximates well the second eigenfunction. Now, we repeat the last experiment but with noise added to an eigenfunction (see Fig. 7). Again the signal is separated into two groups in the time (corresponding to eigenvalue or scale) axis Fig. 7d. We recovered the eigenfunction by filtering out the higher "frequencies". The chosen part to be recovered is the red part in the Fig. 7d.

Conclusions

In this work we formulated necessary and sufficient conditions for general shape preserving flows. We have shown that homogeneous operators admit this condition when the flow is initialized with a nonlinear eigenfunction, with respect to the operator of the flow. Moreover, a closed form solution was formulated for any homogeneous operator and eigenfunction. Consequently, blowup and extinction times were easily derived.

Following these insights, we could extend the TV and one-homogeneous transforms to the p-Laplacian flow, when the flow becomes extinct in finite time, p ∈ (1, 2). It was shown that the flow's decay profile is a polynomial of order 1/(2 -p). Fractional derivatives were suggested to define the nonlinear p-transform. A full theory is still under-way, where currently we have developed results in the classical case of integer derivatives. Numerically, it was shown that the concept extends well to fractional derivatives. We have demonstrated our approach in filtering and reconstruction.

According to the norm equivalence we can write

d

(

  Figure Eigenfunction decay.On the left, eigenfunctions of the p-Laplacian are shown for p = 1.5 in one and two dimensions. On the right, u(t, x 0 ) vs. time and the theoretical plot a(t)f (x 0 ).

3 Figure 2 :

 32 Figure 2: Eigenfunctions with different values of λ and p.

(

Figure 5 :

 5 Figure 5: A demonstration of the general case where 1/(2 -p) is not natural and accordingly the delta function appears at t = 4.9667.

A An alternative proof of Theorem 2, part 2

We compute the decay profile using the first part of the proof of Theorem 2.

The time derivative of P (u(t)), at t = 0, can be expressed as,

Let us recall that

d n dt n P (u). Since u(t) is infinitely differentiable, we can use a Taylor series to express u(t),

Using Eq. ( 27) and the initial condition u(0) = f we get

Assigning g (1) = p -1 yields [START_REF] Li | Numerical methods for fractional calculus[END_REF].

B Finite extinction time of the p-flow

The finite extinction time was proved in Theorem 2 for an eigenfunction as initial condition. In the following proposition we prove that the p-flow vanishes in finite time for arbitrary initial conditions. We show here a different approach based on [START_REF] Burger | Spectral decompositions using one-homogeneous functionals[END_REF] (Proposition 3.12).

Proposition 2 (Finite extinction time in the semi-discrete setting). Let u(t) be the solution to a semi-discrete scale-space flow (6) (continuous in time and discrete in space) with an initial condition u(0) = f . We assume f does not contain any null-space element of P . Then there exists T > 0 such that u(t) = 0 ∀t ≥ T .

Proof.


Pascale Minet
email: pascale.minet@inria.fr

Zied Soua
email: zied.soua@gmail.com

Ines Khoufi
email: ines.khoufi@gmail.com

An Adaptive Schedule for TSCH Networks in the Industry 4.0

Keywords: TSCH network, Scheduling, Network building, Data gathering, Latency, Adaptivity, Industry 4.0

TSCH (Time Slotted Channel Hopping) networks are excellent candidates to support applications of industry 4.0 with latency, energy and reliability requirements. However, such applications have to face topology, traffic or even application changes. How can TSCH networks adapt to such changes? We propose a solution ranging from network construction to data gathering. We show how the schedule can adapt to changing topology, traffic and application conditions. We compute the theoretical bounds with regard to key performance indicators and compare them with the values obtained by NS3 simulation.

I. CONTEXT AND MOTIVATION

The expected fourth industrial revolution, more known as Industry 4.0 [START_REF] Kagermann | Recommendations for implementing the strategic initiative Industrie 4.0: Final report of the Industrie 4.0 working group[END_REF], aims at producing higher quality products in a more efficient and flexible way. It relies on the Internet of Things that will enable a smarter monitoring and more selforganizing capabilities than traditional factories.

To produce better quality products and improve monitoring in Industry 4.0 [START_REF] Hermann | Design principles for industrie 4.0 scenarios[END_REF], strong requirements in terms of latency, robustness and power autonomy have to be met by the networks supporting the Industry 4.0 applications. The wireless TSCH (Time Slotted Channel Hopping) network specified in the e amendment [START_REF] Ieee Sa | IEEE Standard for Local and metropolitan area networks-Part 15.4: Low-Rate Wireless Personal Area Networks (LR-WPANs) -Amendment 1: MAC sublayer[END_REF] of the IEEE 802.15.4 standard [START_REF]IEEE Standard for Local and metropolitan area networks-Part 15.4: Low-Rate Wireless Personal Area Networks (LR-WPANs)[END_REF] has many appealing properties. Its schedule of multichannel slotted data transmissions ensures the absence of collisions. Because there is no retransmission due to collisions, communication is faster. Since the devices save energy each time they do not take part in a transmission, the power autonomy of nodes is prolonged. Furthermore, channel hopping enables to mitigate multipath fading and interferences.

To increase the flexibility and the self-organizing capacities required by Industry 4.0, the networks have to be able to adapt to changes. These changes may concern the application itself, the network topology by adding or removing devices, the traffic generated by increasing or decreasing the device sampling frequency, for instance. That is why the flexibility of the schedule ruling all networks communications is needed. In this paper, we show how a TSCH network can adapt to such changes. Data gathering transmissions are performed according to a schedule. This schedule is computed by the network coordinator (i.e. CPAN), then installed and activated on all network devices.

In the state of the art, many researchers have investigated the schedule of network communications to meet some requirements in terms of latency [START_REF] Yang | Analysis of Low Latency TSCH Networks for Physical Event Detection[END_REF], [START_REF] Rekik | A Performance Analysis of Orchestra Scheduling for TSCH network[END_REF], energy efficiency [START_REF] Vilajosana | A Realistic Energy Consumption Model for TSCH Networks[END_REF], [START_REF] Choudhury | Beacon Synchronization and Duty-Cycling in IEEE 802.15.4 Cluster-Tree Networks: A Review[END_REF], reliability [START_REF] Minet | Increasing reliability of a TSCH network for the industry 4.0[END_REF], or scalability [START_REF] Zats | Scalability of Time Synchronized Wireless Sensor Networking[END_REF], etc.. The solutions proposed are either centralized like in [START_REF] Khoufi | Scheduling transmissions with latency constraints in an IEEE 802.15.4e TSCH network[END_REF], [START_REF] Soua | MODESA: An optimized multichannel slot assignment for raw data convergecast in wireless sensor networks[END_REF], distributed like in [START_REF]Wave: A distributed scheduling algorithm for convergecast in ieee 802.15.4e tsch networks[END_REF], [START_REF] Accettura | Decentralized traffic aware scheduling in 6tisch networks: Design and experimental evaluation[END_REF] or autonomous like in [START_REF] Rekik | Autonomous and traffic-aware scheduling for tsch networks[END_REF]. The performances achieved strongly depend on the assumptions made concerning the network usage conditions (e.g. traffic load distribution, number of network devices). For instance, centralized solutions are able to compute the schedule minimizing the latency but they require a certain overhead to collect all the information needed to compute the schedule and then to install it on all the devices. As a consequence, they do not scale, unlike distributed solutions. Distributed solutions limit the messages exchanged but are usually unable to provide the optimal schedule and even in some cases are unable to bound how far they are from the optimal. Autonomous solutions can be considered as the extreme case of distributed ones: no message is exchanged between devices, the schedule is built from an implicit knowledge shared between nodes (e.g. knowledge of the cells implicitly assigned to a device for transmission to its parent). The main advantage of such solutions is their great flexibility: they are able to quickly adapt to traffic changes. However, these solutions work well only for very light network loads, for which the collision probability remains very low.

Very few researchers have studied how to schedule the transmissions in a network from the network creation. Devices associate to the TSCH network after having heard a beacon. Since the IEEE802.15.4 standard does not specify how the beacons are advertised, several solutions exist like random advertisement defined in [START_REF] Guglielmo | A performance analysis of the network formation process[END_REF], random vertical and random horizontal presented in [START_REF] Vogli | Fast Join and Synchronization Schema in the IEEE 802.15.4e MAC[END_REF], or the best solutions called MBS [START_REF] Guglielmo | A Model-based Beacon Scheduling algorithm for IEEE 802.15.4e TSCH networks[END_REF] and EDBA [START_REF] Khoufi | An enhanced deterministic beacon advertising algorithm for building TSCH networks[END_REF]. The merits of this paper are to show how an adaptive schedule may be installed and activated on all devices, even in case of high traffic, starting from network construction. It presents both theoretical results and simulation results obtained with NS3 [20]. This paper is organized as follows. Section II deals with the construction of a TSCH network and highlights the main concepts used to minimize the time needed to obtain an operational TSCH network. Section III focuses on data gathering and more precisely the schedule used to transmit data generated throughout the TSCH network to the sink. We show how the schedule adapts to topology, traffic or application changes. In Section IV simulation results obtained with NS3 are reported and compared with theoretical bounds. Finally Section V concludes this paper.

II. NETWORK CONSTRUCTION

A. Schedule and Slotframes

Transmissions within a TSCH network are done according to a global schedule built from the superposition of periodic slotframes. In this paper we distinguish three types of slotframes, ordered by decreasing priority:

• Beacon slotframe consisting of a sequence of beacon slots followed by a sequence of free slots (i.e. empty slots). Each beacon slot is assigned to a unique transmitter device, all other devices are in receive mode. • Data slotframe consisting of a sequence of data slots followed by a sequence of free slots. Each data slot is assigned to a unique transmitter device, only the receivers specified are in receive mode, whereas the other devices are sleeping to save energy. • Shared slotframe consisting of a sequence of shared slots followed by a sequence of free slots. In a shared slot, any device is allowed to transmit, whereas all devices having nothing to transmit are in receive mode. The priorities are such that a Beacon slot wins over a Data slot which wins over a Shared slot which wins over a free slot. To decrease the complexity, the Beacon slotframe size is assumed to be a multiple of both the Data slotframe size and the Shared slotframe size.

As an example, the global schedule resulting from the three Slotframes Beacon, Data and Shared is illustrated in Figure 1. A red slot is a Beacon slot, a green slot is a Data slot and a yellow slot is a Shared slot, whereas a slot with an X-cross denotes a free slot. Since according to the priorities, a Beacon slot wins over a Data slot which wins over a Shared slot which wins over a free slot, we get the following result: slots 0 to 7 are beacon slots; slot 8 is a Shared slot; slots 9 to 16 are Data slots, whereas slot 17 is a free one. Slots 18 to 25 are Data slots and slot 26 is a Shared slot. Then the schedule is repeated from slot 0.

B. Node Association

To join the TSCH network, a device randomly selects a physical channel used by this network and listens to a beacon advertising this network. Since the physical channel on which the beacon is broadcast changes at each beacon slot due to channel hopping, the joining device will eventually hear a beacon sent by one of its neighbors. Upon receipt of a valid beacon, this device gets synchronized with the TSCH network. It then requests its association to the network.

To minimize the overhead induced by the use of long addresses coded over 6 bytes for network devices, short addresses coded over 2 bytes are used. To ensure the uniqueness of short addresses, they are assigned by the CPAN (network coordinator) to each joining device when it is associated. As a consequence, the standardized association procedure is modified to become multihop: the association demand is forwarded from the joining node to the CPAN, and the associated response is forwarded from the CPAN to the joining node. In addition, a unique relative beacon index is assigned to the joining device by the CPAN. This beacon index is used by the joining device to compute the beacon slot in which it will broadcast its beacon.

C. Beacon Broadcast

Beacons are broadcast according to an ordered sequence given by the beacon indexes of devices: first the CPAN, then the nodes at depth 1 in the order specified by their index, then the nodes at depth 2 in the order specified by their index, and so on. Since on the one hand, the beacon index ensures that whatever the depth of any device N , N is the unique device to have this index at this depth, and on the other hand, beacons are broadcast according to the order given by the relative beacon index received from the CPAN, denoted BeaconIndex(N ) for node N . Hence, beacons are transmitted without collisions. The beacon slot offset used by a device N to broadcast its beacon is given by Equation 1:

BeaconSlot(N ) = DevLess(Depth(N ))+BeaconIndex(N )
(1) where Depth(N ) denotes the distance expressed in the number of hops from N to the CPAN, and DevLess(Depth(N )) is the number of devices at a depth strictly less than Depth(N ).

This broadcast of beacons ensures that any node has received the beacon from its parent before broadcasting its own beacon. This property enables a consistent information dissemination throughout the network.

III. SCHEDULING DATA GATHERING TRANSMISSIONS

A. Computation of the Data schedule

The CPAN is in charge of computing the schedule for data gathering. For that purpose, the CPAN needs to know for each TSCH device:

• its neighbors with their depth, • its traffic. When the CPAN has received the traffic and topology notifications from all TSCH devices, it computes the schedule for data gathering. In order to decrease the latency of data gathering, the schedule tries to minimize the total number of slots assigned for data gathering. However, to reduce complexity and to avoid frequent computations of the schedule due to topology changes only, spatial reuse is not applied in cells, where a cell is uniquely defined by the pair (Slot Offset, Channel Offset). No spatial reuse means that a cell is granted to only one transmitter. To compute a schedule with no spatial reuse, it is sufficient to know the traffic and the parent of each node. This decreases the overhead, since any node needs to send only its potential parents instead of all its neighbors to the sink.

In the performance evaluation done, the schedule is computed on the CPAN by the debt-based scheduler providing no spatial reuse. Each TSCH device N , except the CPAN, has a debt that is equal to: Debt(N ) = 0 if N has no data message to transmit = Depth(N ) × Rem(N ) otherwise. where Rem(N ) is the remaining number of data messages N has to transmit. The debt-based scheduler schedules first the transmission of the device with the highest debt. To reduce latency, we assume that data generated in the network are timestamped when they are generated. To minimize the maximum latency that can be achieved by data gathering, Assumption 1 is used:

Assumption 1: Each device transmits first the message with the smallest timestamp. This assumption ensures that when a device has messages to send that are generated locally or have been received from its children, it always selects the message containing the oldest data. In addition, if several devices have the same debt, thedebt-based scheduler selects first the message with the smallest timestamp.

An example of a Data schedule is depicted in Figure 2, where the TSCH network uses 5 channels and the sink has 4 radio interfaces. 24 slots are used to transmit Data. The topology is that depicted in Figure 3 and traffic generation is given in Section IV-A.

B. Schedule Installation

The intuitive idea is to take advantage of the ordered sequence of beacon broadcasts starting with the CPAN and progressively involving all devices by increasing depth, to broadcast the schedule to all devices. Since the Data schedule size may be large, it will be fragmented if needed. Each beacon contains a fragment of the schedule. For that purpose, Information Elements are used according to the TSCH specification. The size of the fragment is bounded by the available space in a beacon, denoted BeaconSpace. The number of fragments F rag to broadcast a schedule of size ScheduleSize is given by:

F rag = ScheduleSize BeaconSpace (2)

C. Schedule Activation

All devices simultaneously activate the new schedule at the Absolute Slot Number (ASN) corresponding to the next beacon slot assigned to the CPAN. Hence, all transmissions, including data, beacons and control messages are consistently done according to this active schedule.

D. Schedule Update

As soon as a device observes a change, it notifies the CPAN, using Shared slots. Notice that since this device may be several hops away from the CPAN, the notification should be forwarded if needed. Upon receipt of a notification, the CPAN waits for a given amount of time in order to collect other notifications, because usually, several devices notify the same change or changes sharing the same cause. Then the CPAN takes the decision to:

• either update the current Data schedule by modifying a limited number of cells. The schedule gives for each cell of the Data slotframe the transmitter and the receiver(s). • or compute a new Data schedule. The new schedule or the update of the current one is broadcast in the next beacon of the CPAN and forwarded to all network devices according to the ordered sequence defined in Subsection II-C. All network devices will then activate this new schedule or the updated one in the next beacon slot of the CPAN.

IV. PERFORMANCE EVALUATION

We used the NS3 simulation tool to evaluate the performances of the solution presented in Sections II and III. For that purpose, we developed and integrated the modules implementing the solution proposed.

A. Simulation parameters

The network is made up of the CPAN denoted as device 0 and 32 other devices. We consider three topologies:

• The Regular topology, where all devices are regularly spread on two circles of range R and 2R, centered on the CPAN, where R denotes the radio range. This configuration is used as a reference to point out the impact of device density on network key performance indicators. • The Irregular1 topology, where all devices are uniformly distributed in three sectors around the CPAN. • The Irregular2 topology, depicted by Figure 3, is similar to the previous one, except that devices are distributed in 5 sectors around the CPAN and two-hop devices are more concentrated, resulting in a higher number of children per parent. The simulation parameters used for the performance evaluation with NS3 are summarized in Table I. We distinguish between common parameters that are used for the three topologies considered, from those strongly depending on the configuration such as for example the number of Data slots. Devices 1 to 14, 18, 31 and 32 generate 400 samples per second, whereas the other devices generate 4 samples per second. The parent-child relationship for the Irregular2 topology is given by Table II.

B. Simulation results and theoretical bounds

In this subsection, we evaluate different times that are key performance indicators of the TSCH network. For each time evaluated, we give both an upper theoretical bound and the value obtained for the configuration described in subsection IV-A. These times are related to network building as well as data gathering and adaptivity to topology/traffic/application changes. For this evaluation, we assume that there is no wireless link/device failure. For network building, we also assume that all devices are simultaneously powered on, leading to a bottleneck on network association: all devices wanting to simultaneously join the network.

1) Network building: The average number of neighbors for the three topologies studied is given by Table III. • receive their first beacon on the channel they randomly chose. This time is measured from the time the device is powered on and averaged. For one-hop devices, this delay is shorter with the Regular and Irregular2 topologies that have the smallest average number of neighbors 1-hop away from the CPAN (see the last column of Table III).

For two-hop devices, the delay is higher with the Regular topology that has the smallest average number of neighbors (see the third column of Table III). This can be explained by the fact that when the average number of neighbors increases, the probability of receiving a beacon on the channel selected by the joining device increases which leads to a smaller waiting time before the first beacon receipt. Furthermore, since only associated devices are allowed to send beacons, two-hop devices have to wait the association of their parent, explaining their higher time. • receive their association confirmation. This time is measured from the time of the first association request of the device considered and averaged. Since the association of two-hop devices is multihop, two-hop devices have a longer association time than one-hop devices. • send their first beacon. This time is measured from the receipt of the association confirmation of the device considered and averaged. The maximum time taken to send their beacon is upper bounded by the duration of two Beacon Slotframes, as verified by the values obtained by simulation with NS3. Finally, the delay needed by network building is depicted by Figure 6. It is the time elapsed between the transmission by the CPAN of its first beacon and the receipt by the last device of its association confirmation. The greatest time consuming step is the association of network devices. All three topologies lead to similar results. The average number of association requests per node is depicted by Figure 6, where one-hop devices send similar numbers of requests in the three topologies considered, whereas two-hop devices send the highest number of requests for the Irregular2 topology, that has the highest average number of neighbors for two-hop devices (see the third column of Table III).

2) Data gathering: The maximum time to install the schedule on all devices, denoted Install, is defined as the time elapsed between the CPAN broadcast of the beacon containing the first fragment of the schedule, and the receipt by the last device of the last fragment of the schedule. We where BeaconSlot(P arent(LastDev)) denotes the beacon slot assigned to the parent of the last device receiving the beacon and T ransDelay is the transmission delay of the last schedule fragment.

As an example, we measured the Install time for the Data schedule depicted in Figure 2. For that schedule, 11 fragments are needed, taking into account a BeaconSpace enabling to transmit at most 10 data links (i.e. cells). We got 4.4947 s which corresponds to the theoretical bound computed previously.: (10 * (43+19) * 10+0.0047)ms = 4.4947/s. This value could be optimized by a better coding of the cell with its transmitter and receiver. The coding of Information Elements used to insert a schedule fragment in the beacon message could also be optimized for a better efficiency.

The maximum time to activate the schedule on all devices, denoted Activate, is defined as the time elapsed between the CPAN broadcast of the first fragment of the schedule and the broadcast of the first CPAN beacon following the broadcast of the last fragment; We then have:

Activate = F rag BeaconF rameSize SlotDuration (3)
As an example, we measured the Activate time for the Data schedule depicted in Figure 2. We got 0.4730 s, which corresponds to the number of fragments 11 times the duration of the Beacon slotframe containing 43 slots of 10 ms.

The maximum latency for data gathering, denoted Latency, is defined as the time elapsed between data generation on a device node and its delivery to the sink. The worst case occurs when data are generated while beacons are transmitted, the device has to wait for the next Data slotframe that is not evicted by beacons. In the worst case and with assumption 1, it takes a number of slots equal to the number of Data slots in the Data slotframe to reach the CPAN. Hence, the upper bound:

Latency =((1 + BeaconU sed DataF rameSize ) DataF rameSize + DataU sed) SlotDuration (5)
where BeaconU sed is equal to the total number of devices in the TSCH network, each of them has its own beacon slot, DataU sed is the number of slots used in the Data schedule.

With the parameters given in Table I, the theoretical upper bound for the latency is equal to 0.99s. Simulation results corroborate this value. Hence, the time needed to adapt the schedule to topology or traffic changes can be bounded and the value obtained remains acceptable by most applications.

V. CONCLUSION

In this paper, we have shown how a TSCH network is autonomously built, supports data gathering and is able to adapt to changes in network topology, traffic and application requirements. The solution proposed preserves the merits of TSCH network, that can be listed hereafter. The time-slotted multichannel medium access enables parallel transmissions on several channels, leading to shorter latency and higher throughputs. In addition, channel hopping mitigates interference and multipath effects. Furthermore, since transmissions are scheduled, a conflict-free schedule is computed by the network coordinator (i.e. the CPAN). Hence, no collision occurs during data gathering. The absence of collision leads to a higher throughput, because there is no retransmission due to collisions. It also preserves nodes power autonomy.

This simple solution based on the coexistence of several slotframes. We distinguish three slotframes, that are the Beacon Slotframe, the Data Slotframe and the Shared Slotframe. The network schedule corresponds to the superposition of the three schedules given by each slotframe, where the slotframe with the highest priority wins.

This solution ensures a collision-free dissemination over the whole network. Beacons are broadcast in sequence by increasing depth of devices. This broadcast is also used to disseminate Data Schedules (new schedule or update).

In addition, this solution is adaptive. Topology, traffic or application changes are notified to the CPAN. Depending on the changes notified, the CPAN updates the current schedule or recomputes a new one. Shared slots are used to cope with unexpected events.

Simulation results obtained with NS3 confirm the theoretical upper bounds computed for network construction and data gathering. Hence, TSCH networks are able to adapt to traffic or topology changes in a reasonable time which is a strong requirement of Industry 4.0 applications. In some further work, we will study how to improve this delay to support the most demanding applications.

Fig. 1 .

 1 Fig. 1. The resulting schedule.

Fig. 2 .

 2 Fig. 2. An example of Data Schedule with 24 slots and 5 channels.

Fig. 3 .

 3 Fig. 3. The Irregular2 topology considered: 33 nodes spread in 5 sectors.

Figures 4 and 5

 5 Figures 4 and 5 depict the average delay obtained on 50 simulation runs by one-hop devices and two-hop devices, respectively, to:

Fig. 4 .

 4 Fig. 4. Delays for one-hop devices.

Fig. 5 .

 5 Fig. 5. Delays for two-hop devices.

Fig. 6 .

 6 Fig. 6. Delays for network building.

  have: Install = ((F rag -1)BeaconF rameSize + BeaconSlot(P arent(LastDev)))SlotDuration + T ransDelay

Fig. 7 .

 7 Fig. 7. Average number of association requests per device.

TABLE I SIMULATION

 I PARAMETERS.

TABLE II PARENT

 II -CHILD RELATIONSHIP FOR THE IRREGULAR2 TOPOLOGY.

	Parent Children
	0	1 to 22
	14	23, 24, 27, 28, 30
	18	25, 29, 31, 32
	17	26

TABLE III AVERAGE

 III NUMBER OF NEIGHBORS FOR THE THREE TOPOLOGIES.

			Avg # of neighbors
	Topology	for 1-hop dev. for 2-hop dev.	1-hop from the CPAN
				for 1-hop dev.
	Regular	11.5	7.2	11
	Irregular1 14.31	8.5	10.95
	Irregular2 17.13	19.5	12.36