Keywords: Differential equations, Differentiators, Parameter estimation

In this article, we propose a robust method to compute the output of a fractional linear system defined through a linear fractional differential equation (FDE) with timevarying coefficients, where the input can be noisy. We firstly introduce an estimator of the fractional derivative of an unknown signal, which is defined by an integral formula obtained by calculating the fractional derivative of a truncated Jacobi polynomial series expansion. We then approximate the FDE by applying to each fractional derivative this formal algebraic integral estimator. Consequently, the fractional derivatives of the solution are applied on the used Jacobi polynomials and then we need to identify the unknown coefficients of the truncated series expansion of the solution. Modulating functions method is used to estimate these coefficients by solving a linear system issued from the approximated FDE and some initial conditions. A numerical result is given to confirm the reliability of the proposed method.

Fractional order differentiation by integration: an application to fractional linear systems D.Y. Liu * T.M. Laleg-Kirati * O. Gibaru * * , * * *

INTRODUCTION

Fractional calculus were introduced in many fields of science and engineering long time ago. It was first developed by mathematicians in the middle of the ninetieth century [START_REF] Ross | The development of fractional calculus 1695-1900[END_REF]. During the past decades, fractional calculus has gained great interest in several applications [START_REF] Oldham | The Fractional Calculus[END_REF], for instance in control engineering [START_REF] Oustaloup | The CRONE control of resonant plants: application to a flexible transmission[END_REF], signal processing [START_REF] Magin | On the fractional signals and systems[END_REF], and finance [START_REF] Mainardi | Fractional calculus and continuous-time finance II: the waiting-time distribution[END_REF], etc. Unlike classical differential equations, fractional differential equations can better describe some natural phenomena and especially some dynamical systems [START_REF] Ciesielski | Numerical simulations of anomalous diffusion[END_REF][START_REF] Metzler | The random walk's guide to anomalous diffusion: a fractional dynamics aproach[END_REF]. This is due essentially to their memory and hereditary properties. Many studies have been done in the existence and uniqueness of solutions to fractional differential equations (see, e.g. [START_REF] Podlubny | Fractional Differential Equations[END_REF][START_REF] Kilbas | Theory and Applications of Fractional Differential Equations[END_REF][START_REF] Deng | Existence and uniqueness of solutions of initial value problems for nonlinear fractional differential equations[END_REF]).

In the case of disturbed source, the behavior of solution has been studied in [START_REF] Diethelm | The Analysis of Fractional Differential Equations[END_REF]. However, it is in general difficult to find analytic solutions to fractional differential equations, reliable approximations and numerical techniques are then essential. One can distinguish two main categories of methods for solving fractional differential equations: frequency-domain methods and time-domain methods. Frequency-domain methods include the popular Laplace transform method [START_REF] Podlubny | Fractional Differential Equations[END_REF] and Fourier transform method [START_REF] Kilbas | Theory and Applications of Fractional Differential Equations[END_REF]. They usually provide analytical solutions for fractional differential equations. Concerning time-domain methods, an increasing number of numerical schemes are being developed for both linear and nonlinear fractional differential equations. In nonlinear case, these methods include Lubich's difference methods, Diethelm's method based on quadrature, Adams-Bashforth-Moulton method, Homotopy-perturbation method, Adomian's decomposition method (see [START_REF] Magin | On the fractional signals and systems[END_REF]).

During the last few years, many authors applied function approximation theory to solve linear fractional differential equations with constant coefficients. Solutions were approximated by truncated series expansions involving different bases, such as orthogonal polynomials [START_REF] Doha | Efficient Chebyshev spectral methods for solving multiterm fractional orders differential equations[END_REF] and wavelets [START_REF] Rehman | The Legendre wavelet method for solving fractional differential equations[END_REF][START_REF] Li | Numerical solution of fractional differential equations using cubic B-spline wavelet collocation method[END_REF], with unknown coefficients. Thus, the problem of solving a fractional differential equation was transformed into a problem of identifying the coefficients of an equation. Then, different methods were considered to solve this problem by solving a linear system of algebraic equations. These methods include operational matrix method [START_REF] Rehman | The Legendre wavelet method for solving fractional differential equations[END_REF], collocation method [START_REF] Li | Numerical solution of fractional differential equations using cubic B-spline wavelet collocation method[END_REF], and spectral tau method [START_REF] Doha | Efficient Chebyshev spectral methods for solving multiterm fractional orders differential equations[END_REF]. Unlike the spectral tau method used in [START_REF] Doha | Efficient Chebyshev spectral methods for solving multiterm fractional orders differential equations[END_REF], the operational matrix method in general uses polynomials to approximate the solutions of fractional differential equations, and the fractional derivatives or the fractional integrals of these solutions. But the fractional derivatives and the fractional integrals of a classical polynomial are not polynomials. Moreover, the collocation method in-terpolates source functions at collocation points which are usually not irregularly spaced. This is a constraint for many applications. Also, the method is not robust against noisy source functions. Very recently, Legendre polynomials and spectral tau method were used in [START_REF] Bhrawy | A quadrature tau method for variable coefficients fractional differential equations[END_REF] to solve fractional differential equations with time-varying coefficients. In this paper, we use Jacobi polynomials and the modulating function method to extend the methods used in [START_REF] Bhrawy | A quadrature tau method for variable coefficients fractional differential equations[END_REF].

There are practical situations where the source functions are measured and corrupted by an additive noise. However, none of the above methods has considered this case. In this paper, we are interested in the fractional order differentiation by integration with Jacobi polynomials. This method consists in estimating the fractional derivative of an unknown signal by an integral formula involving Jacobi polynomials and the noisy observation of the unknown signal [START_REF] Liu | Fractional order differentiation by integration with Jacobi polynomials[END_REF][START_REF] Liu | Fractional order differentiation by integration and error analysis in noisy environment[END_REF]. It generalizes the method of differentiation by integration [START_REF] Mboup | Numerical differentiation with annihilators in noisy environment[END_REF][START_REF] Liu | Differentiation by integration with Jacobi polynomials[END_REF] from integer order to fractional order. Let us recall that the method of differentiation by integration with Jacobi polynomials is the generalization of the Lanczos generalized derivative [START_REF] Lanczos | Applied Analysis[END_REF] (p. 324) (1956) in noisy case. The differentiators proposed in [START_REF] Mboup | Numerical differentiation with annihilators in noisy environment[END_REF] were originally introduced using the algebraic parametric techniques [START_REF] Fliess | Compression différentielle de transitoires bruités[END_REF][START_REF] Fliess | An algebraic framework for linear identification[END_REF][START_REF] Mboup | Parameter estimation for signals described by differential equations[END_REF][START_REF] Liu | An error analysis in the algebraic estimation of a noisy sinusoidal[END_REF][START_REF] Liu | Nonasymptotic fractional order differentiators via an algebraic parametric method[END_REF][START_REF] Ushirobira | Algebraic parameter estimation of a biased sinusoidal waveform signal from noisy data[END_REF][START_REF] Perruquetti | An algebraic approach for Humane posture estimation in the sagital plane using accelerometer noisy signal[END_REF] which exhibited good robustness properties with respect to corrupting noises without the need of knowing their statistical properties [START_REF] Fliess | Analyse non standard du bruit[END_REF].

One of the methods that have been proposed to solve parameter identification problems is the modulating function method. The latter was pioneered by Shinbrot [START_REF] Shinbrot | On the analysis of linear and nonlinear systems[END_REF] in 1957 to solve the parameter identification problem by using integral transformation. The idea was motivated by Laplace and Fourier transforms. Essentially, the use of modulating functions allows to transform a differential equation, involving input-output noisy signals on a specified time interval, into a sequence of algebraic equations. Moreover, it avoids the derivatives of the noisy signals and annihilates the effects of initial conditions, thus allows the direct use of noisy signals in the integral expression. Hence, this method takes advantage from the low-pass filtering property of modulating functions integrals [START_REF] Liu | Parameters estimation of a noisy sinusoidal signal with timevarying amplitude[END_REF][START_REF] Fliess | Analyse non standard du bruit[END_REF]. It was used to parameter identification for nonlinear systems, time-varying systems and noisy sinusoidal signals [START_REF] Preising | Theory and application of the modulating function method I: Review and theory of the method and theory of the spline type modulating functions[END_REF][START_REF] Ungarala | Time-varying system identification using modulating functions and spline models with application to bio-processes[END_REF][START_REF] Liu | Parameters estimation of a noisy sinusoidal signal with timevarying amplitude[END_REF].

The aim of this paper is to propose a robust method to find the output of fractional linear systems, defined through linear fractional differential equations with timevarying coefficients, where the input can be noisy. For this purpose, we use the fractional order differentiation by integration with Jacobi polynomial and the modulating functions method. The calculus are made from the noisy input and some initial conditions. This paper is organized as follows: in Section 2, we recall the method of fractional order differentiation by integration. It is used in Section 3 to solve fractional differential equations together with the modulating function method. In Section 4, we apply the proposed method to an example. Finally, some conclusions and perspectives are given in Section 5.

FRACTIONAL ORDER DIFFERENTIATION BY INTEGRATION

In this section, we recall the main idea behind fractional order differentiation by integration [START_REF] Liu | Fractional order differentiation by integration with Jacobi polynomials[END_REF][START_REF] Liu | Fractional order differentiation by integration and error analysis in noisy environment[END_REF], which generalizes the method of differentiation by integration with Jacobi polynomials [START_REF] Mboup | Numerical differentiation with annihilators in noisy environment[END_REF][START_REF] Liu | Differentiation by integration with Jacobi polynomials[END_REF] from the integer order to the fractional order. Before doing so, we recall the definitions and some useful properties of the Caputo fractional derivative and the Jacobi orthogonal polynomials, respectively.

Caputo fractional derivative

Unlike classical integer order derivatives, there are several definitions for fractional derivatives which are in general not equivalent with each others [START_REF] Podlubny | Fractional Differential Equations[END_REF][START_REF] Kilbas | Theory and Applications of Fractional Differential Equations[END_REF]. In this paper, we use the Caputo fractional derivative.

Let f ∈ C l (R) with l ∈ N * , where C l (R) denotes the set of the l-times continuously differentiable functions defined on R. Then, the Caputo fractional derivative (see [START_REF] Podlubny | Fractional Differential Equations[END_REF] p. 79) of f is defined as follows:

∀ t ∈ R * + , D α t f (•) := 1 Γ(l -α) t 0 (t -τ ) l-α-1 f (l) (τ ) dτ, (1)
where 0 ≤ l -1 < α ≤ l, and Γ(z) = ∞ 0 exp(-x) x z-1 dx is the Gamma function (see [START_REF] Abramowitz | editeurs. Handbook of mathematical functions[END_REF] p. 255). Hence, the Caputo fractional derivative of an n th order polynomial f (t) = t n can be given as follows (see [START_REF] Podlubny | Fractional Differential Equations[END_REF] p. 72):

∀ t ∈ R + , D α t f (•) =    0, if n < α, Γ(n + 1) Γ(n + 1 -α) t n-α , if n ≥ α. (2)
In the following are some useful proprieties of the Caputo fractional derivative: ∀ t ∈ R * + , • Linearity (see [START_REF] Podlubny | Fractional Differential Equations[END_REF] p. 91):

D α t {λ 1 f 1 (•) + λ 2 f 2 (•)} = λ 1 D α t f 1 (•) + λ 2 D α t f 2 (•), (3)
• Scale change (see [START_REF] Oldham | The Fractional Calculus[END_REF] p. 76):

D α λ t f (•) = 1 λ α D α t f (λ •), (4)
where

λ ∈ R * + , λ 1 , λ 2 ∈ R and f, f 1 , f 2 ∈ C l (R).

Jacobi orthogonal polynomials

The n th order shifted Jacobi orthogonal polynomial (see [START_REF] Abramowitz | editeurs. Handbook of mathematical functions[END_REF] p. 775) can be defined on [0, 1] as follows:

P (µ,κ) n (τ ) := n j=0 n + µ j n + κ n -j (τ -1) n-j τ j , (5)
where µ, κ ∈] -1, +∞[, and

n+µ j = Γ(n+µ+1) Γ(n+µ-j+1) Γ(j+1) , n+κ n-j = Γ(n+κ+1)
Γ(κ+j+1) Γ(n-j+1) . Let g 1 and g 2 be two continuous functions defined on [0, 1], then we define the following scalar product (see [START_REF] Abramowitz | editeurs. Handbook of mathematical functions[END_REF] p. 774):

g 1 (•), g 2 (•) µ,κ = 1 0 w µ,κ (τ ) g 1 (τ ) g 2 (τ ) dτ, (6)
where w µ,κ (τ ) = (1τ ) µ τ κ is the associated weighted function. Thus, the norm of the Jacobi polynomial P (µ,κ) n is given by: P

(µ,κ) n 2 µ,κ = Γ(µ+n+1) Γ(κ+n+1) Γ(µ+κ+n+1) Γ(n+1) (2n+µ+κ+1) .
Then, by applying the linearity (3) to (5) and using (2), the Caputo fractional derivative of the Jacobi polynomial can be given in the following lemma. Lemma 1. The α th (α ∈ R + ) order derivative of the n th order Jacobi orthogonal polynomial P (µ,κ) n defined in ( 5) is given as follows:

∀ τ ∈ [0, 1], D α τ P (µ,κ) n (•) =      0, if n < α, n j=0 min(n-j,n-l) i=0 c µ,κ,n,j,i Γ(n -i + 1) Γ(n -i + 1 -α) τ n-i-α , if n ≥ α,
where c µ,κ,n,j,i = (-1) i n+µ j n+κ n-j n-j i

, and l ∈ N * with 0 ≤ l -1 < α ≤ l.

Fractional Jacobi differentiator

Let y ∈ C l (R), and y ̟ be a noisy observation of y on an interval

I = [0, h] ⊂ R + : ∀ t ∈ I, y ̟ (t) = y(t) + ̟(t), (7)
where ̟ is an integrable noise1 . We are going to estimate the α th order derivative of y using the observation y ̟ .

Let us ignore the noise for a moment. By taking the Jacobi orthogonal series expansion of y(h•) ([36] p. 6), we have:

∀ ξ ∈ [0, 1], y(hξ) = +∞ i=0 P (µ,κ) i (•), y(h•) µ,κ P (µ,κ) i 2 µ,κ P (µ,κ) i (ξ). (8)
We take the truncated Jacobi orthogonal series expansion of y(h•) to approximate y on I by the following N th (N ∈ N) order polynomial:

∀ ξ ∈ [0, 1], D (0)
h,µ,κ,N y(hξ

) := N i=0 λ (µ,κ) i P (µ,κ) i (ξ), (9)
where λ

(µ,κ) i = P (µ,κ) i (•),y(h•) µ,κ P (µ,κ) i 2 µ,κ
. Then, we calculate the fractional derivative of the polynomial D

h,µ,κ,N y(•) such that we can approximate the fractional derivative of y. We denote the α th order derivative of

D (0) h,µ,κ,N y(•) by D (α)
h,µ,κ,N y(•). By applying the scale change property (4) and the linearity (3), to [START_REF] Kilbas | Theory and Applications of Fractional Differential Equations[END_REF], we obtain:

∀ ξ ∈ [0, 1], D α hξ y(•) ≈D (α) h,µ,κ,N y(hξ) :=D α hξ D (0) h,µ,κ,N y(•) = 1 h α D α ξ D (0) h,µ,κ,N y(h•) = 1 h α N i=0 λ (µ,κ) i D α ξ P (µ,κ) i (•). (10)
Consequently, by substituting y by its noisy observation y ̟ in [START_REF] Deng | Existence and uniqueness of solutions of initial value problems for nonlinear fractional differential equations[END_REF] we obtain the following fractional order differentiator:

∀ ξ ∈ [0, 1], D α hξ y(•) ≈ D (α) h,µ,κ,N y ̟ (hξ) := 1 h α N i=0 λ(µ,κ) i D α ξ P (µ,κ) i (•), (11)
where λ(µ,κ)

i = P (µ,κ) i (•),y ̟ (h•) µ,κ P (µ,κ) i 2 µ,κ
. By expanding the scalar product •, • µ,κ , we get the following integral form for

D (α) h,µ,κ,N y ̟ (hξ): ∀ ξ ∈ [0, 1], D (α) h,µ,κ,N y ̟ (hξ) = 1 h α 1 0 Q µ,κ,α,N (τ, ξ) y ̟ (hτ ) dτ, (12)
where Q µ,κ,α,N (τ, ξ) = w µ,κ (τ )

N i=0 P (µ,κ) i (τ ) P (µ,κ) i 2 µ,κ D α ξ P (µ,κ) i (•).
According to [START_REF] Diethelm | The Analysis of Fractional Differential Equations[END_REF], the differentiator

D (α)
h,µ,κ,N y ̟ (hξ) uses an integral involving Jacobi polynomials to estimate the fractional order derivative of a noisy signal. We call this differentiator fractional Jacobi differentiator as in [START_REF] Liu | Fractional order differentiation by integration with Jacobi polynomials[END_REF][START_REF] Liu | Fractional order differentiation by integration and error analysis in noisy environment[END_REF], and this method fractional order differentiation by integration.

The estimation error for the fractional Jacobi differentiator in the noisy case can be divided into two sources: the truncated term error which is due to the truncated terms in the Jacobi orthogonal series expansion of y(h•), and the noise error contribution which is due to the noise ̟. These errors have been studied in [START_REF] Liu | Fractional order differentiation by integration and error analysis in noisy environment[END_REF]. Moreover, it has been also showed in [START_REF] Liu | Fractional order differentiation by integration and error analysis in noisy environment[END_REF] how to choose the design parameters N , κ and µ for the fractional Jacobi differentiator so as to reduce these errors.

APPLICATION TO FRACTIONAL LINEAR SYSTEMS

Fractional linear systems

In this section, we consider a fractional linear system defined by a Fractional Differential Equation (FDE) in the following form:

∀ t ∈ I = [0, h] ⊂ R + , L i=0 a i (t) D αi t y(•) = M j=0 b j (t) D βj t u(•), ( 13
)
where y is the output, u is the input, L, M ∈ N, a i = 0, b j = 0 are time-varying coefficients and assumed to be known. We suppose that α i ,

β j ∈ R + are such that 0 ≤ α 0 < α 1 < • • • < α L and 0 ≤ β 1 < β 2 < • • • < β M .
Let u ̟ be a noisy observation of u on I. The goal is to find an approximation of the output y by solving the fractional differential equation defined in [START_REF] Doha | Efficient Chebyshev spectral methods for solving multiterm fractional orders differential equations[END_REF] with the noisy input u ̟ and the following initial conditions:

y (i) (0) = d i ∈ R, for i = 0, . . . , ⌈α L ⌉ -1, (14)
where ⌈α L ⌉ denotes the smallest integer greater than or equal to α L . The existence and uniqueness of solution of the initial value problem ( 13)-( 14) have been studied in [START_REF] Deng | Existence and uniqueness of solutions of initial value problems for nonlinear fractional differential equations[END_REF] in noise-free case. We are going to numerically solve this problem by using the fractional Jacobi differentiator and the modulation functions method in the noisy case.

Application of fractional Jacobi differentiator

In order to solve the initial value problem ( 13)-( 14), we use the polynomial defined in [START_REF] Kilbas | Theory and Applications of Fractional Differential Equations[END_REF] to approximate the solution y. Moreover, the fractional order derivatives of y can be approximated by using the fractional Jacobi differentiator given in [START_REF] Deng | Existence and uniqueness of solutions of initial value problems for nonlinear fractional differential equations[END_REF]. On the one hand, the left side of the FDE (13) can be approximated as follows:

∀ ξ ∈ [0, 1], L i=0 a i (hξ) D αi hξ y(•) ≈ L i=0 a i (hξ) D (αi) h,µ,κ,N y(hξ) = L i=0 a i (hξ) h αi N j=0 λ (µ,κ) j q (αi) µ,κ,j (ξ) = N j=0 λ (µ,κ) j L i=0 a i (hξ) h αi q (αi) µ,κ,j (ξ), (15)
where q

(αi) µ,κ,j (ξ) = D αi ξ P (µ,κ) j (•).
On the other hand, we denote the right side of the FDE [START_REF] Doha | Efficient Chebyshev spectral methods for solving multiterm fractional orders differential equations[END_REF] by

U (ξ) = M j=0 b j (hξ) D βj hξ u(•).
Then, U (•) can be estimated by using the fractional Jacobi differentiator given in [START_REF] Diethelm | The Analysis of Fractional Differential Equations[END_REF] and the noisy input u ̟ . Hence, the FDE ( 13) is approximated by the following equation:

∀ ξ ∈ [0, 1], N j=0 λ (µ,κ) j Q (µ,κ) j (ξ) = Ũ (ξ), (16) where Q (µ,κ) j (ξ)
= L i=0 a i (hξ) h αi q (αi) µ,κ,j (ξ)
, and Ũ (•) is an estimation of U (•). Consequently, the problem of solving the FDE ( 13) is transformed into a problem of identification of unknown coefficients λ (µ,κ) j

, for j = 0, 1, . . . , N . We are going to apply the modulating functions method to solve this problem.

Modulating function method for FDEs

Let us define the following subset of the fractional orders in (13): B := {β j ; β j ∈ N, j = 0, . . . , M }, and β * = max βj ∈B β j .

If B = ∅, then β * = 0. Then, we take a class of functions {g i } N +1 i=⌈αL⌉ which satisfy the following conditions:

g i ∈ C β * ([0, 1]), g (j)
i (0) = g

(j) i (1) = 0, ∀ j = 0, 1, . . . , β * -1. (17)
These functions are called modulating functions [START_REF] Preising | Theory and application of the modulating function method I: Review and theory of the method and theory of the spline type modulating functions[END_REF].

By multiplying equation ( 16) by each function g i and by integrating them on [0, 1], we then obtain the following system:

N j=0 λ (µ,κ) j 1 0 g i (ξ) Q (µ,κ) j (ξ) dξ = 1 0 g i (ξ) Ũ (ξ) dξ, (18)
for i = ⌈α L ⌉, . . . , N . Then, the right side of ( 18) is calculated in the following way: [START_REF] Doha | Efficient Chebyshev spectral methods for solving multiterm fractional orders differential equations[END_REF], then by applying β j times integration by parts and the property [START_REF] Liu | Fractional order differentiation by integration with Jacobi polynomials[END_REF] which annihilates the boundary values, we obtain:

• If β j ∈ R + \N in (13), then we estimate D βj hξ u(•) by the fractional Jacobi differentiator D (βj ) h,µ,κ,N u ̟ (hξ) in Ũ (ξ). • If β j ∈ N in
1 0 g i (ξ) b j (hξ) u (βj ) (hξ) dξ = (-1) βj h βj 1 0 d βj dξ βj {g i (ξ) b j (hξ)} u(hξ) dξ. (19)
Hence, by using the noisy input u ̟ we get:

1 0 g i (ξ) b j (hξ) u (βj ) (hξ) dξ ≈ (-1) βj h βj 1 0 d βj dξ βj {g i (ξ) b j (hξ)} u ̟ (hξ) dξ. (20)
Thus, we avoid the estimation of the integer order derivative of the noisy signal u ̟ in [START_REF] Liu | Fractional order differentiation by integration and error analysis in noisy environment[END_REF]. The estimation error can be reduced (see [START_REF] Liu | Error analysis of Jacobi derivative estimators for noisy signals[END_REF][START_REF] Liu | Convergence Rate of the Causal Jacobi Derivative Estimator[END_REF][START_REF] Liu | Error analysis for a class of numerical differentiator: application to state observation[END_REF] for more details on the error analysis for the integer order differentiation by integration method).

Consequently, due to the low-pass filtering property of the fractional Jacobi differentiator and the modulating function method, the noise effect in the right side of ( 18) can be reduced. Now, we take the values of the integer order derivatives of the approximation polynomial defined in ( 9) at ξ = 0 to estimate the initial conditions. By applying similar calculations than those developed for [START_REF] Deng | Existence and uniqueness of solutions of initial value problems for nonlinear fractional differential equations[END_REF] we then get:

N j=0 λ (µ,κ) j q (i) µ,κ,j (0) h i = d i , for i = 0, . . . , ⌈α L ⌉ -1. (21)
Consequently, by using ( 18) and ( 21), we can solve the following linear system to estimate the coefficients λ

(µ,κ) i : M (µ,κ) N    λ (µ,κ) 0 . . . λ (µ,κ) N    =    d 0 . . . d N    , (22)
where d i = 1 0 g i (ξ) Ũ (ξ) dξ, for i = ⌈α L ⌉, . . . , N , and

M (µ,κ) N (i, j) =        q (i) µ,κ,j (0) h i , for 0 ≤ i ≤ ⌈α L ⌉ -1, 0 ≤ j ≤ M, 1 0 g i (ξ) Q (µ,κ) j (ξ) dξ, for ⌈α L ⌉ ≤ i ≤ N, 0 ≤ j ≤ M.

NUMERICAL RESULTS

In order to demonstrate the reliability of our method, we consider the following FDE (see [START_REF] Bhrawy | A quadrature tau method for variable coefficients fractional differential equations[END_REF]):

     4 i=0 a i (t) D αi t y(•) = f (t), y(0) = 2, y ′ (0) = 0, (23)
where a 0 (t) = t

1 5 , a 1 (t) = t 1 4 , a 2 (t) = t 1 3 , a 3 (t) = t 1 2 , a 4 (t) ≡ 1, α 0 = 0, α 1 = 0.333, α 2 = 1, α 3 = 1.234, α 4 = 2, and y(t) = 2 -1 2 t 2 .
We need to choose an input such that f satisfies [START_REF] Fliess | An algebraic framework for linear identification[END_REF]. The choice of this input is not unique. In this example, we take u(•) = y(•) such that We can see this discrete noisy input in Figure 1.

f (t) =
Since the noisy input is given in the discrete case, we apply the trapezoidal rule to approximate the integrals in ( 12) and ( 18). This produces a numerical error in the estimated output. Since y and u are two polynomials of degree 2, we take N = 2 in the fractional Jacobi differentiators. Thus, there is no truncated term errors in the used fractional Jacobi differentiators. Moreover, we take µ = κ = 0 and the following modulating function: g 2 (ξ) = ξ 3 (1ξ) 7 for ξ ∈ [0, 1]. The estimation error in the estimated output can be divided into two sources: the numerical error and the noise error contribution. Since y is a polynomial, this estimation error can also be a polynomial. We can see the estimated output and the corresponding estimation error in Figure 2 and Figure 3.

CONCLUSION

In this paper, we have proposed a method to determine the output of a fractional linear system from some initial conditions, where the input can be noisy. The system is defined through a linear fractional differential equation with time-varying coefficients. An integration-based fractional order differentiation method has been used to estimate the fractional derivative of noisy signals. It has been shown that this method exhibits good robustness properties with respect to a corrupting noise. We have also shown that with this method, the problem of solving the linear fractional differential equation becomes a problem of identifying parameters. We have used modulating function method for the identification of the parameters which is suitable for time-varying linear systems involving input-output noisy signals. Numerical results confirms the efficiency of the proposed method. In order to reduce the estimation error with appropriate modulating functions, we will think to give some error analysis in a future work. Questions related to observability and parameter identification for a fractional linear system will also be considered with such approach.

4 j=0aFig. 1 .Fig. 2 .

 412 Fig.1. The noise-free input and its noisy observation.

Fig. 3 .

 3 Fig. 3. Absolute estimation error.

More generally, the noise is a stochastic process, which is integrable in the sense of convergence in mean square (see[START_REF] Liu | Error analysis of Jacobi derivative estimators for noisy signals[END_REF]).


Adrien Couetoux

Olivier Teytaud

Hassen Doghmen

Learning a Move-Generator for Upper Confidence Trees

We experiment the introduction of machine learning tools to improve Monte-Carlo Tree Search. More precisely, we propose the use of Direct Policy Search, a classical reinforcement learning paradigm, to learn the Monte-Carlo Move Generator. We experiment our algorithm on different forms of unit commitment problems, including experiments on a problem with both macrolevel and microlevel decisions.

Introduction

Monte-Carlo Tree Search (MCTS) [START_REF] Coulom | Efficient Selectivity and Backup Operators in Monte-Carlo Tree Search[END_REF] is a versatile algorithm for Markov Decision Processes (MDP) or games. It is elegant (can be described in a few lines), versatile (can be applied in many settings such as MDP, games, stochastic games), and can work with or without expert information in the playouts (hence a great success in general game playing [START_REF] Sharma | Knowledge generation for improving simulations in uct for general game playing[END_REF]). It uses the framework of bandits [START_REF] Kocsis | Bandit based Monte-Carlo planning[END_REF] which does not depend on a particular application. It is moderately efficient when the number of time steps is big, but surprisingly stable on high-dimensional problems.

It has been greatly improved by including Progressive Widening and Double Progressive Widening [START_REF] Coulom | Computing elo ratings of move patterns in the game of go[END_REF][START_REF] Chaslot | Progressive Strategies for Monte-Carlo Tree Search[END_REF], RAVE values [START_REF] Gelly | Combining online and offline knowledge in UCT[END_REF], Blind Values [START_REF] Couetoux | Improving the exploration in upper confidence trees[END_REF], and handcrafted Monte-Carlo moves [START_REF] Wang | Modifications of UCT and sequence-like simulations for Monte-Carlo Go[END_REF][START_REF] Lee | The Computational Intelligence of MoGo Revealed in Taiwan's Computer Go Tournaments[END_REF]. A crucial component is the Monte-Carlo move generator, also known as the playout generator.

In this paper, we focus on the addition of specialized Monte-Carlo moves, i.e. we modify default policy, to help dealing with stochastic planning problems. Finding a default policy that is optimal for all instances of a stochastic problem can be extremely difficult and time consuming. The solution we propose here is to apply a Direct Policy Search to the available default policy. This way, even an initially poor default policy can be improved to fit different instances of one stochastic planning problem.

In Section 2 we describe existing algorithms (Monte-Carlo Tree Search in Section 2.1, and existing algorithms for improving Monte-Carlo move generators in Section 2.2). We also introduce Direct Policy Search for improving a Monte-Carlo move generator in Section 2.3. In Section 3, we experiment our algorithms on three forms of a unit commitment problem (Section 3.1), and on an investment problems (Section 3.2).

In this section we will present the vanilla Monte-Carlo Tree Search algorithm (Section 2.1), Direct Policy Search, and Monte-Carlo move generators improvements (Sections 2.2 and 2.3).

Monte-Carlo Tree Search and Upper Confidence Trees

Monte-Carlo Tree Search (here presented in the framework of a MDP) consists in simulating plenty of series of decisions as long as we have time before choosing an action, and keeping statistics of these games. What follows, is the formal description of the state of the art continuous MCTS, i.e. MCTS with Double Progressive Widening (MCTS-DPW), as seen in [START_REF] Couetoux | Continuous Upper Confidence Trees[END_REF]. As the reader can see, it mainly requires two things: (i) a transition function, capable of simulating what happens when an action a is taken in state s, and returns a new state s ′ and a reward r. (ii) a default policy ϕ, capable of returning an action a, given a state s. When nothing is specified, it is assumed that this function returns a random action, following a random distribution that covers the entire set of feasible actions in state s.

MCTS algorithm with DPW and default policy ϕ

Input: a state S. Output: an action a. Initialize: ∀s, nbSims(s) = 0 while Time not elapsed do // starting a simulation. s = S. while s is not a terminal state and nbSims(s) > 0 do

Apply DPW in state s. Let s ′ be the state given by DPW. s = s ′ end while while s is not a terminal state // {happens when a non final and new state s is visited} do

Choose action a, according to ϕ Simulate action a; get a new state s ′ s = s ′ end while Get a reward r = Reward(s) // s is a final state, it has a reward. For all states s in the simulation above, let r nbV isits(s) (s) = r. end while Return the action which was simulated most often from S.

The algorithm therefore relies on a Monte-Carlo move generator, also called default policy, ϕ. The default policy can be a simple random uniform generator (when no expertise is available for making more reasonable simulations), but handcrafted functions can perform better.

Heuristics and Monte-Carlo move generators

Whereas in the 2-player case, it is known that making a Monte-Carlo generator stronger (stronger in the sense: as a stand-alone policy), does not necessarily make the MCTS built on top of it stronger (see [START_REF] Wang | Modifications of UCT and sequence-like simulations for Monte-Carlo Go[END_REF]), we conjecture that in the one-player case it is usually quite efficient.

The recent improvements in the world of Computer Go basically comes from improvements of the Monte-Carlo move generator, implemented so that the Monte-Carlo simulator does not contradict life&death known results; Zen, CrazyStone, Pachi, are examples of such strong programs, around 2 Dan for short time settings and 4 Dan for long time settings. Other tools have been proposed as generic solutions for learning Monte-Carlo move generators:

-Simulation balancing [START_REF] Silver | Monte-carlo simulation balancing[END_REF][START_REF] Huang | Monte-carlo simulation balancing in practice[END_REF] has been proposed for automatically learning the Monte-Carlo move generator in 2-player games.

-PoolRave [START_REF] Rimmel | Biasing Monte-Carlo Simulations through RAVE Values[END_REF], in which the Monte-Carlo move is replaced, with a fixed probability p ∈ (0, 1), by a move uniformly drawn among the c moves with best RAVE score in the last node of the simulation with at least k simulations.

-Contextual Monte-Carlo [START_REF] Rimmel | Multiple Overlapping Tiles for Contextual Monte Carlo Tree Search[END_REF] in which the Monte-Carlo move-generator is improved by online learning a tile-based value function.

These tools are efficient, but the main successes in Monte-Carlo Tree Search nonetheless come from handcrafted Monte-Carlo move generators. In this paper, we used a specialized Monte-Carlo move generator, chosen specifically on our main target problem, as well as a less specialized function. We improve them by Direct Policy Search (Section 2.3).

Direct Policy Search for generating Monte-Carlo Move Generators

Direct Policy Search (DPS) is an approach very different from Upper Confidence Tree; it is based on selecting a policy among a parametric family of policies by optimization of its parameters. The pseudo-code is as follows: Direct Policy Search is usually applied offline, i.e. a single θ is obtained once and for all. However, optimizing Θ to maximize θ → Simulate(s, M DP, π(θ)) specifically for the current state s for which we look for a decision is possible. We apply DPS and use the obtained policy π( θ) as a Monte-Carlo move generator in our MCTS.

The paper in [START_REF] Bengio | Using a financial training criterion rather than a prediction criterion[END_REF] proposes to apply DPS (the terminology in the paper is different, but it is essentially DPS) based on a heuristic function obtained by experts, by smoothing the heuristic and adding parameters in it (the smoothing is here for making the problem easier to optimize). This is our approach in the rest of this paper, except that we do not smooth policies as the randomized nature of our problems make the objective function smooth enough. We use selfadaptation [START_REF] Meyer-Nieberg | Self-adaptation in evolutionary algorithms[END_REF] as a noisy optimization algorithm. As a summary, our algorithm is as follows for choosing a move in state s within time t: Procedure OptimisticHeuristics(s, φ, t, M DP ) Input: a state s, a time t, a parametric family of policies φ θ . Output: an action a. Apply DPS with time budget t/2 for choosing θ (use warm start if possible) Apply MCTS with time budget t/2 for choosing action a.

Here, we compare the performances of different sequential decision making algorithms. Namely, we implemented vanilla MCTS, MCTS with a fixed default policy, MCTS with a default policy improved online by DPS, and DPS alone.

We made experiments on three different forms of the unit commitment problem, and on a more general energy management problem called bilevel.

Unit commitment problem

We work on a stock management problem, from [START_REF] Teytaud | Including Ontologies in Monte-Carlo Tree Search and Applicationsan Open Source Platform[END_REF].

The main points in the problem are that: (1) Unit Commitment problems can not be solved efficiently by traditional methods; these problems are usually simplified so that classical methods, like Bellman's stochastic dynamic programming, can be applied. The motivation of our work on Unit Commitment by Monte-Carlo Tree Search methods is that we want to work without simplifying too much the model. ( 2) Unit Commitment problems exist at many time scales (from milliseconds, up to years for hydroelectric stocks or tenths of years if investments are included) and many dimensionalities (from a few stocks to thousands of state variables), depending on the scope under analysis. We here work on small scale problems for the sake of statistical significance (working on our full problems requires by far too much time for reproducing runs tenths of times).

In this paper, we will consider three variants of the unit commitment problem. The significant difference between these three variants is the way the stocks are connected. In the first one, they are lined up on a one dimension chain (we will call it the one river problem). In the second one, they are linked so that they form a binary tree, with the root being the last stock that the water goes through (we will call it the "binary rivers problem"). Finally, the third one is simply a random arrangement of the stocks, with one single constraint: no cycles are allowed.

Two different heuristics for the Unit commitment problem. The expert parametrized heuristic that we use has been designed using knowledge about the problem, to make it particularly efficient on the one river variant of the unit commitment problem. On the other hand, the naive heuristic uses almost no knowledge about the problem. Given a state s, it requires the current time to go t, and the average demand at the current time step D avg . We provide below pseudo-codes of both heuristics.

Expert heuristic

Input: a state S, of dimension N . Parameter: a vector θ of dimension 3. Default value is [1, 0, 1]. Information required from the problem:

-D(t): expected electricity demand during time step t.

-D timeT oGo (t): expected total demand after time step t.

-T SA(s, timeT oGo): total stock available (this assume a 1 river structure).

-T I(timeT oGo, averageInf lows): expected total usable water from inflows (this assume a 1 river structure).

Output: an action a.

initialize:

total water available = T W A = (θ0 + θ1 × timeT oGo) × (T SA + T I) -x = production by hydroelectricity = 0 -increaseW ater ← true -S available = 0≤i≤N -1 si 2. while (increaseW ater and x < S available ) si being the current level of stock i do define the marginal cost mc of increasing water, approximated as

mc = IC(x, s, t, D(t)) + θ2 × LT C(x, T W A, t, D timeT oGo (t))
where:

-IC(x, s, t, D(t)) is negative; it is the marginal benefit associated to the reduction of thermal production. -LT C(. . . ) is the sum of thermal production cost, if expected total demand D timeT oGo , decreased by the total production from the water stocks if equally distributed on the time steps to go, is produced thermally.

if marginal cost mc > 0 then then increaseW ater ← f alse else

x ← x + 1. end if 3. end while 4. Compute q, the ratio min(0, x S available ) 5. Return the action vector a defined as follows: ∀0 ≤ i ≤ N -1, ai = q.Li Naive heuristic, polynomial with degree m Input: a state S, of dimension N . Parameter: a vector θ of dimension m + 1. Default values are [1, 0, . . . , 0] Information required from the problem: t the remaining time steps, and Davg the average demand after the current time step Output: an action a.

Compute total amount of water to use

Wuse = max(0, Davg.(θ0 + θ1t + • • • + θm+1t m )) 2.
Given S and the current level Li of each stock i, W available = 0≤i≤N -1 Li 3. Compute q, the ratio min(0, Wuse W available ) 4. Return the action vector a defined as follows:

∀0 ≤ i ≤ N -1, ai = q.Li
Experimental results on the Unit Commitment problem. We present here the results obtained on all three variants of the unit commitment problem. Each time, we compared the following algorithms: (i) vanilla MCTS, as presented in Section 2.1, (ii) MCTS-naive, a MCTS using the naive heuristic as a default policy,(iii) MCTS-expert, a MCTS using the expert heuristic as a default policy, (iv) MCTS-naive-DPS, a MCTS using the naive heuristic improved by DPS, (v) MCTS-expert-DPS, a MCTS using the expert heuristic improved by DPS and when relevant, (vi) the non tuned naive and expert heuristics. The x axis shows the time budget allocated per decision made, in logarithmic scale, and went from 0.01 second to 2.56 second. The y axis shows the average reward. Each average reward was computed using 1000 runs. Error bars show the 95% confidence intervals. The higher the reward, the better the algorithm performed. It should be noted that the rewards cannot be compared between different variants of the unit commitment problem. Indeed, only the connections between the stocks change, and changing this changes the amount of water effectively available. Our results for the 1-river problem and for the binary rivers problem are shown on the left side and the right side of Fig. 1, respectively. We did not plot the results of the naive heuristic, that scored -150000 and -380000 respectively (far below other methods), for the sake of readability. In both experiment, MCTS-naive-DPS and MCTS-expert-DPS outperform by a factor of at least 100 the third placed algorithm, MCTS-expert. MCTS-naive and MCTS vanilla share the fourth and fifth places. Our results on the random rivers problem are shown in Fig. 2. In this experiment, the most significant difference in the results is that MCTS-naive-DPS is about 10 times faster than MCTS-expert-DPS.

Over all three versions of the unit commitment problem, the most efficient and robust version has been MCTS-naive-DPS. Even on the one river problem, that the expert heuristic was particularly well tuned for, we could not see huge benefits from using it as a parametric function for DPS.

Experiments on the investment problem

In this section we experiment our algorithm on a problem (simplified from [START_REF] Teytaud | Including Ontologies in Monte-Carlo Tree Search and Applicationsan Open Source Platform[END_REF]) as follows: -At each time step, we decide investments; there is a limited amount of money to invest, and investments must be distributed over 7 different possible infrastructures. There are therefore 7 decision variables for each time step. -At each time step, a lower level problem (the management of the energy production system) is built and solved, and its cost is the cost of the current transition of the investment problems. -There are 10 time steps, the last one has a strong influence because it reflects the long-term. Our results compare the following strategies:

a heuristic which gives a constant ratio of the investment on each possible infrastructure (the parameters of the heuristic are this proportions); the default parametrization is the same ratio for all infrastructures; -DPS on a "sum of Gaussians" policy (parameters: positions of the Gaussians, widths, associated decisions; see [START_REF] Teytaud | Including Ontologies in Monte-Carlo Tree Search and Applicationsan Open Source Platform[END_REF]; -DPS on a "neural network" policy (parameters: weights, theresholds; see [START_REF] Teytaud | Including Ontologies in Monte-Carlo Tree Search and Applicationsan Open Source Platform[END_REF]); -DPS on a "sum of Gaussians" policy, added to the heuristic with default parametrization; -DPS on a "neural network" policy, added to the heuristic with default parametrization; -MCTS, on top of each of the above.

Results are presented in Fig. 3.

Conclusion

We combine DPS and MCTS. The DPS provides the Monte-Carlo simulator of the MCTS. The resulting algorithm, has no free parameter and outperforms by far the vanilla MCTS. We use human expertise at two levels: (i) For partial observation handling, i.e. the belief state estimation was handcrafted, so that the problem is essentially a MDP rather than a partially observable MDP. The details of this are beyond the scope of this paper. (ii) In the Monte-Carlo move generator, because in spite of nice and interesting efforts in the literature, no generic algorithm, in the current state of the art, can define a Monte-Carlo move generator as efficiently as a human expert (in the case of Go, but also in the case of unit commitment problems). Nonetheless our DPS could strongly improve the heuristic by optimizing its parameters. We agree with the traditional statement that MCTS is surprisingly efficient when no human expertise is available, but we clearly see that human expertise was an easy key for a speed-up 100, as well as human expertise is the key of recent progress in MCTS for the classical challenge of the game of Go.

Importantly, the need for human expertise is considerably reduced by the use of DPS for optimizing the heuristics, so that our results are a step towards generic MCTS tools.

Fig. 1 .

 1 Fig.1. Performances of different variants of MCTS on the 1 river unit commitment problem (left) and the binary rivers (right), with 7 stocks, 24 time steps. Y axis shows the reward (the higher the better).

Fig. 2 .

 2 Fig.2. Performances of different variants of MCTS on the randomly connected unit commitment problem (7 stocks, 24 time steps).Y axis shows the reward (the higher the better).

Fig. 3 .

 3 Fig. 3. Results on the energy investment problem. The five DPS curves (curves 1 to 5) are very close to each other; results are better than for the heuristic alone, and versions without the heuristic are almost the same as versions with the heuristic. The MCTS+DPS+neural network was the most efficient strategy, outperforming MCTS+DPS+neural network+heuristic. The sums of Gaussians require more time for learning, hence the poor results for moderate budgets.