Khadidja Chaib Draa

Ali Zemouche

Marouane Alma

Holger Voos

Mohamed Darouach

K Chaib Draa
email: kh.chaibdraa@gmail.com

Nonlinear observer-based control with application to an anaerobic digestion process

Keywords: Anaerobic digestion, LMI approach, Observer Design, Reference trajectory tracking

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction and Preliminaries

Introduction

Anaerobic Digestion (AD) is a biological process through which the organic waste is converted into a mixture of gaseous that is called biogas. The latter can be valorized as electricity, heat, biofuel, or it can be injected into the natural gas grid. Considering electricity and natural gas grids, the potential of biogas is substantial, and allows for the integration of higher shares of variable renewable energy sources in the electricity grid. Hybrid systems are currently being internationally developed. They combine biogas plants with weather dependent technologies for balancing the grids [START_REF] Scarlat | Biogas: Developments and perspectives in europe[END_REF]. However, due to the complexity of the AD process and fluctuations of input influents [START_REF] Jimenez | Instrumentation and control of anaerobic digestion processes: a review and some research challenges[END_REF], continuous monitoring and design of sophisticated control strategies are required for the safety of the process, and improving efficiency.

With the aim to enhance the AD process performances, different control strategies have been proposed in the literature. Often, the adopted control strategy depends on the model complexity, the available measurements and the design criteria (pollutant minimization, product maximisation or digester stabilisation). With respect to the proposed control for reduced analytical AD models (two-steps "acidogenesis-methanogenesis" models), we may cite the control of bicarbonates alkalinity concentration in the digester through the addition of stimulating substrates to the process [START_REF] Hess | Modélisation de la Qualité du Biogas Produit par un Fermenteur Méthanogène et sa Stratégie de Régulation en Vue de sa Valorisation[END_REF], [START_REF] Palacios-Ruiz | Regulation of volatile fatty acids and total alkalinity in anaerobic digesters[END_REF] and [START_REF] Méndez-Acosta | A robust control scheme to improve the stability of anaerobic digestion processes[END_REF]. In [START_REF] Hess | Modélisation de la Qualité du Biogas Produit par un Fermenteur Méthanogène et sa Stratégie de Régulation en Vue de sa Valorisation[END_REF], the linearizing control was employed to enhance the biogas quality, while in [START_REF] Palacios-Ruiz | Regulation of volatile fatty acids and total alkalinity in anaerobic digesters[END_REF] and [START_REF] Méndez-Acosta | A robust control scheme to improve the stability of anaerobic digestion processes[END_REF], the input to output linearizing control was used to stabilize the digester. In the first and second control strategies, the magnitude of the added input was assumed very small so the general structure of the model does not change. In other words, the added input affects only the state for which it has been added (for example, added alkalinity affects only dynamics of the alkalinity concentration). Such an assumption makes the control simpler to design, especially the input to output linearizing control, which requires a nonlinear transformation of the system. Indeed, the mentioned assumption relaxes complexity of the required transformation. However, neglecting the effect of an additional input in the model dynamics is not very convenient. In parallel, we may also cite the Model Predictive Control (MPC) that has been proposed in [START_REF] Mauky | Model predictive control for demand-driven biogas production in full scale[END_REF] to control the biogas production for a demand-driven electricity production. The idea is to optimize the plant feeding according to a fluctuating timetable of energy demand. The control was applied to a full-scale research plant, and has shown satisfactory results. However, the analytical stability proof of the closed loop system is yet difficult to prove.

Combining the ideas from [START_REF] Hess | Modélisation de la Qualité du Biogas Produit par un Fermenteur Méthanogène et sa Stratégie de Régulation en Vue de sa Valorisation[END_REF], [START_REF] Palacios-Ruiz | Regulation of volatile fatty acids and total alkalinity in anaerobic digesters[END_REF], [START_REF] Méndez-Acosta | A robust control scheme to improve the stability of anaerobic digestion processes[END_REF] and [START_REF] Mauky | Model predictive control for demand-driven biogas production in full scale[END_REF], where the alkalinity addition is employed to stabilize the reactor and enhance the biogas quality, and the plant feeding is used to optimize the production, we will propose a control strategy to track an admissible reference trajectory. This trajectory mimics a desired biogas production that corresponds to the operator objectives. Moreover, the control strategy takes into account all the process dynamics including those that are not accessible for measurement. Indeed, it is well known that the key parameters of the AD process, such as bacteria, are difficult and costly to measure; and this is the reason why state observers have been designed repeatedly in the literature. The asymptotic observer [START_REF] Bastin | On-line Estimation and Adaptive Control of Bioreactors, Process Measurement and Control[END_REF] is quite simple to design, and does not require knowledge of some specific nonlinear functions (kinetic functions). Nevertheless, it is very sensitive to the model uncertainties, and its convergence rate depends on the operational conditions. Thus, it has been extended to interval observers [START_REF] Bernard | Closed loop observers bundle for uncertain biotechnological models[END_REF], which have the advantage of using reliable measurements that are nonlinear functions of the state vector. The interval observers estimate the intervals where the system state is laying when the model is subject to large uncertainties. However, generally the rate of convergence is partially tunable, and it is not easy to exploit the intervals for control [START_REF] Chaib-Draa | An LMI-based H ∞ discrete-time non linear state observer design for an anaerobic digestion model[END_REF]. Furthermore, the Kalman filter [START_REF] Haugen | State estimation and model-based control of a pilot anaerobic digestion reactor[END_REF], [START_REF] Rocha-Cozatl | State and input estimation of an anaerobic digestion reactor using a continuous-discrete unknown input observer[END_REF] has been proposed to enhance the convergence rate of estimation error towards zero. Although the Kalman filter has shown suitable results in different chemical applications, unfortunately, the convergence of estimation errors to zero is not guaranteed. Within this context, we may cite the high gain observer [START_REF] Gauthier | A simple observer for nonlinear systems, applications to bioreactors[END_REF], [START_REF] Farza | Nonlinear observers for parameter estimation in bioprocesses[END_REF], which has a fast convergence rate to the model state variables; but its synthesis is complex, and it is very sensitive to noise [START_REF] Lombardi | Implementation of observer for on-line estimation of concentration in continuous-stirred membrane bioreactor: Application to the fermentation of lactose[END_REF]. In this paper, we will use the LMI-based nonlinear state observer proposed in [START_REF] Chaib-Draa | LMI-based H ∞ nonlinear state observer design for anaerobic digestion model[END_REF] due to its robustness, systematic implementation and fast convergence. Then, the estimated states are used in a feedback control, which accounts for all the process dynamics. The synthesis of the observer-based controller parameters is done by solving two different LMI conditions. From feasibility point of view, this kind of separation principle for nonlinear Lipschitz systems (proposed in this paper) provides enhanced LMIs, when compared to the existing LMI techniques in the literature. These LMIs ensure the exponential convergence of the estimation error towards zero, and guarantee that the tracking error satisfies a given H ∞ criterion. The proof of convergence is done by using rigorous mathematical arguments that are easy to follow. To get enhanced LMI conditions, we used the famous Young's relation in a convenient way, and the standard Lipschitz condition is exploited in a different form to take into account all the components of the nonlinearity of the system.

The remainder of the paper is organized as follows. In Section 2, we pose the problem of observer-based reference tracking control. Then, in Section 3, we provide new LMI conditions which ensure the exponential convergence to zero of the estimation error and the H ∞ asymptotic stability of the tracking error. Furthermore, in Section 4, we give a detailed description on how to apply the designed control to the AD process model. We also run a numerical simulation to emphasize the effectiveness of the proposed control scheme. Finally, we conclude the paper in Section 5.

Notation and Preliminaries

This section is devoted to the notation used throughout this paper and some useful lemmas.

Notation

The following notation will be used throughout this paper:

• ( ) is used for the blocks induced by symmetry;

• A T represents the transposed matrix of A;

• I r represents the identity matrix of dimension r;

• for a square matrix S, S > 0 (S < 0) means that this matrix is positive definite (negative definite);

• the set Co(x, y) = {λx + (1 -λ)y, 0 ≤ λ ≤ 1} is the convex hull of {x, y};

• e s (i) = 0, ..., 0, i th 1 , 0, ..., 0 s components T ∈ R s , s ≥ 1 is a vector of the canonical basis of R s .

Useful lemmas

We present hereafter two lemmas, which will be used to get tractable and less conservative LMI conditions. Lemma 1.1 ( [16]). Let ϕ: R n → R q be a differentiable function. Then, the following items are equivalent:

• ϕ is a globally γ ϕ -Lipschitz function;

• For all x, y ∈ R n , there exist z i ∈ Co(x, y), z i = x, z i = y, scalar constants a ij , b ij , and functions ψ ij : R n → R such that:

ϕ(x) -ϕ(y) = q,n i,j=1 ψ ij (z i )H ij x -y (1)
and

a ij ≤ ψ ij z i ≤ b ij , (2)
where

ψ ij (z i ) ∂ϕ i ∂x j (z i ), H ij = e q (i)e T n (j).
Notice that we have introduced this lemma in order to simplify the presentation of the design methodology. Indeed, throughout this paper, we will exploit (1)-( 2) instead of a direct use of the Lipschitz property.

Lemma 1.2 ([16]

). Let X and Y be two matrices of appropriate dimensions. Then, for any given symmetric positive definite matrix S of appropriate dimension, the following inequality holds:

X T Y + Y T X ≤ 1 2 X + SY T S -1 X + SY . (3)
This lemma will be very useful to enhance the feasibility of the LMI conditions.

Formulation of the Observer-Based Tracking Problem

One class of control design problems consists of planning and following a reference trajectory while the system variables are not fully measurable. Among these systems we find the AD process, where the bacteria that are key parameters of the process are not accessible for measurement. To overcome this problem, a state observer can be included in the control design, allowing reconstruction of the state vector. However, the choice of the state observer is crucial since its dynamics affect the stability of the closed loop system.

Concerning the two-step models of the AD process, we propose to use the nonlinear observer designed in [START_REF] Chaib-Draa | LMI-based H ∞ nonlinear state observer design for anaerobic digestion model[END_REF], due to its exponential convergence and systematic implementation. Moreover, to keep the results general and convenient for other nonlinear systems belonging to the same class of systems as the AD process, we represent the AD model by the following system ẋ = A(ρ u )x + Gγ(x)

+ Bu y = Cx (4)
where the state vector x ∈ R n , the input u ∈ R q and the linear output measurements y ∈ R p . The parameter ρ u ∈ R s is an L ∞ bounded and known parameter and the affine matrix A(ρ u ) is expressed under the form

A(ρ u ) = A 0 + s j=1 ρ j A j (5)
with

ρ min ≤ ρ u ≤ ρ max (6)
which means that the parameter ρ u belongs to a bounded convex set for which the set of 2 s vertices can be defined by:

V ρ u = ∈ R s : j ∈ {ρ j,min , ρ j,max } . (7)
The matrices

A i ∈ R n×n , G ∈ R n×m and C ∈ R p×n are constant.
The nonlinear function γ : R n -→ R m is assumed to be globally Lipschitz and can always be written under the detailed form:

γ(x) = m i=1 γ i ( ϑ i H i x) (8)
where

H i ∈ R n i ×n .
Remark 1. The parameter ρ u depends only on the control input vector u. Nevertheless, it can depend also on the output y and any other exogenous variable. However, without loss of generality, we used only u, which is motivated by the anaerobic digestion model considered in this paper.

In order to reconstruct the state vector of the system (4), we use the following state observer:

ẋ = A(ρ u )x + G m i=1 γ i ( θi ) + Bu + L(ρ u ) y -C x (9) with θi = H i x + K i (ρ u )(y -C x), (10)
and

L(ρ u ) = L 0 + s j=1 ρ u j L j , K i (ρ u ) = K i0 + s j=1 ρ u j K ij . ( 11
)
where x is the estimate of the system state x. The matrices L i ∈ R n×p and K ij ∈ R n i ×p are the observer gains, to be computed so that the estimation error e = x -x [START_REF] Gauthier | A simple observer for nonlinear systems, applications to bioreactors[END_REF] turns to be exponentially stable. Since γ(.) is globally Lipschitz, then according to Lemma 1.1 there exist z i ∈ Co(ϑ i , θi ), functions φ ij : R n i -→ R, and constants a ij , b ij , so that

G(γ(x) -γ(x)) = m,n i i,j=1 φ ij (z i )H ij ϑ i -θi (13)
and

a ij ≤ φ ij (z i ) ≤ b ij , (14)
where

φ ij (z i ) = ∂γ i ∂ϑ j i (z i ), H ij = Ge T m (i)e n i (j) (15)
For shortness, we set φ ij φ ij (z i ). Without loss of generality, we assume that a ij = 0 for all i = 1, . . . , m and j = 1, . . . , n i . For more details about this, we refer the reader to [START_REF] Arcak | Observer-based control of systems with slope-restricted nonlinearities[END_REF].

Since

ϑ i -θi = (H i -K i (ρ u )C
) e, then dynamics of the estimation error are given by

ė = A L (ρ u ) + m,n i i,j=1 φ ij H ij H i -K i (ρ u )C e (16)
with

A L (ρ u ) = A(ρ u ) -L(ρ u )C (17)
As already explained, dynamics of the estimation error affect the stability of the closed loop system. Thus, before discussing the stability analysis of the chosen observer and the computation of its parameters, we first present the observer based trajectory tracking control. Then, we discuss the stability conditions of the closed loop system composed of the system, the observer and the controller.

A feasible trajectory for the system (4) is a pair (x d , u d ) that satisfies the differential equation and generates the desired trajectory:

ẋd = A(ρ u d )x d + Gγ(x d ) + Bu d y d = Cx d (18)
where x d and u d represent the desired state and the desired input, respectively. Assuming that it is possible to find a feasible trajectory for the system, then it is possible to search for controllers u = f (x, x d , u d ) that track the desired trajectory. Besides, due to the absence of key measurements, we use the state estimate instead. That is the tracking control reads

u = -K(ρ u d )(x -x d ) + u d (19)
where

K(ρ u d ) = K 0 + s j=1 ρ u d j K j (20)
Let us define the tracking error by

x = x -x d (21)
Its dynamic can be easily obtained as

ẋ = A(ρ u ) -BK(ρ u d ) + m,n i i,j=1 ϕ ij (t)H ij H i x + BK(ρ u d )e + A(ρ u ) -A(ρ u d ) x d ω(t) (22)
where

ϕ ij ∂γ i ∂ϑ j i (ν i ), with ν i ∈ Co(x, x d ) and
ϕ ij ≤ ϕ ij ≤ ϕ ij (23)
The aim consists of finding the controller and observer gain matrices, so that the tracking error x satisfies the following

H ∞ criterion x L n 2 ≤ µ ω 2 L n 2 + η x0 e 0 2 (24)
where µ > 0 is the gain from ω to x, and η > 0. In the next section, we provide the stability conditions that satisfy this objective.

New LMI Procedure to Design the Observer-Based Controller

In this section we present a kind of separation principle for nonlinear systems. Since the dynamics [START_REF] Zemouche | Circle criterion-based H ∞ observer design for Lipschitz and monotonic nonlinear systems-Enhanced LMI conditions and constructive discussions[END_REF] do not depend on the reference tracking error x and the functions φ ij are bounded, then we can study the convergence of the estimation error e separately, and use it in dynamics of the tracking error as a bounded disturbance exponentially converging towards zero. The following theorem provides the synthesis conditions expressed in term of LMIs.

Design of the observer gains

We first state a theorem which provides the LMI conditions ensuring the exponential convergence of the estimation error towards zero. Theorem 3.1. For a given positive scalar β, if there exist symmetric positive definite matrices Q, S i , i = 1, . . . , n, and matrices X i , X ij , i = 1, . . . , n and j = 1, . . . , s, of appropriate dimensions, such that the following LMI conditions (25) are fulfilled for system (4):

     A Q, X, + βQ Π Π 1 . . . Π m ( ) -ΛS      ≤ 0, ∀ ∈ V ρ (25)
then the observer (9)-( 11) with gains

L j = Q -1 X T j , K ij = S -1 i X T ij (26)
guarantees the exponential convergence of the estimation error e(t), following the inequality:

e(t) ≤ λ max (Q) λ min (Q) e 0 e -β 2 t , (27)
where

A Q, X, = A T 0 Q + QA 0 -C T X 0 -X T 0 C + s j=1 j A T j Q + QA j -C T X j -X T j C (28)
Π i = M 1 i Q, S i . . . M n i i Q, S i (29)
M j i Q, S i = QH ij + H T i S i -C X i0 + s j=1 j X ij (30) Λ = block-diag Λ 1 , ..., Λ m (31)
Λ i = block-diag 2 φi1 I n i , . . . , 2 φin i I n i (32) S = block-diag S 1 , . . . , S m (33)
S i = block-diag S i , . . . , S i . ( 34
)
Proof. Consider the quadratic Lyapunov function

V (e) = e T Qe, Q = Q T > 0. (35)
Then, its derivative V (e) along the trajectories ( 16) is given by

V (e) = e T   A L (ρ) + m,n i i,j=1 φ ij H ij H K i T Q + Q A L (ρ) + m,n i i,j=1 φ ij H ij H K i   e (36)
with

A L = A(ρ) -L(ρ)C, H K i = H i -K i (ρ)C ( 37
)
The derivative of Lyapunov function (36) is negative if

A T L (ρ)Q + QA L (ρ) + m,n i i,j=1 φ ij    QH ij X T ij H K i Y i +Y T i X ij    < 0 (38)
From Lemma 1.2, we know that for any symmetric positive definite matrices S ij , we have

X T ij Y i + Y T i X ij ≤ 1 2 X ij + S ij Y i T S -1 ij X ij + S ij Y i M j i Q,S ij (39)
Since the matrix Y i does not depend on the index j and depends on the same K i (ρ k ), then to get an LMI, we need to put

S ij = S i , ∀(i, j) (40)
Consequently, from (23) and the fact that without loss of generality a ij = 0, inequality (38) is satisfied if

A T L (ρ)P + P A L (ρ) + m,n i i,j=1 M j i Q, S i T - 2 b ij S i -1 M j i Q, S i ≤ 0. ( 41
)
Therefore, from Schur lemma and with the change of variables

X j = L T j Q and X ij = K T ij S i , inequality (41) is equivalent to   A T L (ρ)P + P A L (ρ) Π 1 . . . Π m ( ) -ΛS   ≤ 0. (42)
Since (42) is affine in ρ u , then from the convexity principle [START_REF] Boyd | Linear matrix inequalities in system and control theory[END_REF] and the notations (28)-(34) we deduce that the LMIs (25) guarantees the inequality V (e) + βV (e) ≤ 0,

which leads to

V (e) ≤ V (e 0 )e -βt . (44)
From the fact that V (e) satisfies the inequality

λ min (Q) e(t) 2 ≤ V (e(t)) ≤ λ max (Q) e 0 2 e -βt (45)
then we get the exponential convergence of the estimation error as follows:

e(t) ≤ λ max (Q) λ min (Q) e 0 e -β 2 t , (46)
where β 2 is the convergence rate of the error e(t). This ends the proof of Theorem 3.1.

Design of the controller parameters

Hereafter, we present a second theorem, which provides new LMI conditions constituting the second step of the proposed observer-based controller synthesis methodology. Theorem 3.2. If there exist symmetric positive definite matrices P, Z ij , i, j = 1, . . . , n, and matrices Y i of appropriate dimensions, such that the convex optimization problem (47) is solvable:

min(µ) subject to (48) (47)       Θ Σ Σ 1 . . . Σ m ( ) -ΛZ       ≤ 0, ∀ ∈ V ρ ( 48
)
then the controller ( 19)-(20) with gains

K j = Y j P -1 , j = 1, . . . , s, (49)
guarantees the following H ∞ performance criterion on the tracking error x:

x L n 2 ≤ µ ω K(ρ u d )e 2 L n 2 + λ max P x0 2 , (50)
where

Θ =   Θ 11 P 0 ( ) -I n   , Θ 11 = A P, Y, I n B ( ) -µI n+q A P, Y, = PA T 0 + A 0 P -Y 0 B -BY 0 (51) + s j=1 j PA T j + A 0 P -Y j B -BY j Σ i = N 1 i P, Z i1 . . . N n i i P, Z in i (52) N j i P, Z ij =   PH T i 0 0   +   H ij 0 0   Z ij (53) Λ = block-diag Λ 1 , ..., Λ m (54)
Λ i = block-diag 2 φi1 I n i , . . . , 2 φin i I n i (55) Z = block-diag Z 1 , . . . , Z m (56)
Z i = block-diag Z i1 , . . . , Z in i . ( 57
)
Proof. Consider the Lyapunov function

V (x) = xT P -1 x. ( 58
)
By calculating the derivative of V (x) along the trajectories of (22), we get

V + xT x -µω T ω = x ω P -1 D + D T P -1 + I n P -1 B ( ) -µI n+q x ω T (59)
with

D = A(ρ u d ) -BK(ρ u d ) + m,n i i,j=1 ϕ ij (t)H ij H i , (60)
B = [I n B] and ω = ω K(ρ u d )e . (61)
Inequality (59) holds if the following matrix is positive definite

P -1 D + D T P -1 + I n P -1 B ( ) -µI n+q < 0, (62)
which is equivalent to

DP + PD T + P 2 B ( ) -µI n+q < 0 (63)
by pre-and post-multiplying (62) by

P 0 0 I n .
Using the Schur lemma, it follows that inequality (63) is equivalent to

  DP + PD T B P ( ) -µI n+q 0 ( ) 0 -I n   < 0, (64)
which can be rewritten under the form:

  (A(ρ u d ) -BK(ρ u d )) P + P (A(ρ u d ) -BK(ρ u d )) T B P ( ) -µI n+q 0 ( ) 0 -I n   + m,n i i,j=1 ϕ ij (t)            X T i       PH T i 0 0       Y ij H T ij 0 0 +Y T ij X i            < 0. (65)
From Lemma 1.2, for any symmetric positive definite matrices Z ij , the inequality (65) holds if the following one is fulfilled:

  (A(ρ u d ) -BK(ρ u d )) P + P (A(ρ u d ) -BK(ρ u d )) T B P ( ) -µI n+q 0 ( ) 0 -I n   + m,n i i,j=1 ϕ ij (t) 2 X i + Z ij Y ij T Z -1 ij X i + Z ij Y ij N j i P,Z ij < 0. (66)
Therefore, using the notation of Theorem 3.2, applying the Schur lemma, and from the convexity principle, the inequality (66) is equivalent to (48). This means that

V (x) + xT x -µω T ω ≤ 0, (67)
which leads to (50) by integrating (67) from 0 to +∞. Indeed, by integrating (67), we get which is equivalent to (50).This ends the proof.

V (+∞) -V (0) + +∞ 0 xT (t)x(t)dt -µ +∞ 0 ωT (t)ω(t)dt ≤ 0 Since V (+∞) > 0 and V (0) ≤ λ max P x0 2 ,

Fulfilment of the H ∞ criterion (24)

Now we complete the convergence analysis part of the proposed H ∞ observer-based tracking design methodology. The last and final step of the convergence proof can be summarized in the next proposition. Proposition 3.1. If there exist symmetric positive definite matrices P, Q, Z ij , S i , i, j = 1, . . . , n, matrices Y i , X i , X ij , i = 1, . . . , n; j = 1, . . . , s of appropriate dimensions, and a positive scalar β, such that the following two items hold:

i) The LMI conditions (25) are fulfilled;

ii) the convex optimization problem (47) is solvable, then, the observer-based controller (19)-(20) with gains given by (26) and (49), guaratees the H ∞ criterion (24) on x with µ returned by (47) and η given by

η max λ max (P), 1 β λ max (Q)λ max (K (ρ u d )K(ρ u d )) λ min (Q) . ( 68
)
Proof. Since LMIs (25) hold, then the estimation error satisfies (27). Then, we have

ω K(ρ u d )e 2 L n 2 ≤ ω 2 L n 2 + λ max K (ρ u d )K(ρ u d ) e 2 L n 2 ≤ ω 2 L n 2 + λ max K (ρ u d )K(ρ u d ) λ max (Q) λ min (Q) e 0 2 +∞ 0 e -βτ dτ = ω 2 L n 2 + λ max K (ρ u d )K(ρ u d ) β λ max (Q) λ min (Q) e 0 2 . ( 69
)
On the other hand, the solvability of the convex optimization problem (47) means that the performance criterion(50) holds. It follows that

x L n 2 ≤ µ ω 2 L n 2 + max λ max (P), λ max (Q)λ max (K (ρ u d )K(ρ u d )) βλ min (Q) x0 e 0 2 . ( 70
)
This ends the proof of Proposition 3.1.

Remark 2. The proposed LMI-based design technique is easy to implement numerically. Indeed, once the parameters of a given model are known, it suffices to solve the proposed LMIs to get the observer-based controller gains ensuring the H ∞ criterion. To solve such LMIs, we use some available numerical software toolboxes like Matlab LMI Toolbox by using YALMIP. As for the computational aspect and the implementation complexity, there are no obstacles for real-time applications because the observer-based controller parameters are computed offline by solving the LMIs.

Application to an Anaerobic Digestion Process

Mathematical Model of the AD Process

AD modelling has been widely investigated in the literature. Often, the resulted model is driven by the application objectives, waste nature and its composition, the available data and its reliability. In this paper, being motivated by the control and observer design for the AD processes, we will use the same model considered in [START_REF] Chaib-Draa | LMI-based H ∞ nonlinear state observer design for anaerobic digestion model[END_REF] which is based on the AM2 model [START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF] with the inclusion of one additional input representing the addition of stimulating alkalinity as shown in Figure 1. The resulted model structure reads

               ẋ1 = -k 1 µ 1 (x 1 )x 2 + u 1 S 1in -ux 1 ẋ2 = (µ 1 (x 1 ) -αu)x 2 ẋ3 = k 2 µ 1 (x 1 )x 2 -k 3 µ 2 (x 3 )x 4 + u 1 S 2in -ux 3 ẋ4 = (µ 2 (x 3 ) -αu)x 4 ẋ5 = k 4 µ 1 (x 1 )x 2 + k 5 µ 2 (x 3 )x 4 + u 1 C in -ux 5 -q c (x) ẋ6 = u 1 Z in + u 2 Z ad -ux 6 (71) in F 1 Feeding Tank in in in in Z C S S , , , 2 1 ad Z n i F 2 V out F v 6 5 4 3 2 1 , , , , , x x x x x x m c q q ,

Anaerobic Digester

Added Alkalinity Tank

y 1 = q c (x) y 2 = [x 1 , x 3 , x 6 ] T (72) with µ 1 (x 1 ) = µ 1 x 1 x 1 + k s 1 , µ 2 (x 3 ) = µ 2 x 3 x 3 + k s2 + x 2 3 k i2 ( 73
)
co 2 = x 5 + x 3 -x 6 , φ = co 2 + K H P T + k 6 k La µ 2 (x 3 )x 4
(74)

q c (x) = k La [co 2 -K H P c (x)], q m (x) = k 6 µ 2 (x 3 )x 4 (75) P c (x) = φ -φ 2 -4K H P T co 2 2K H (76)
where x 1 (g/l) is the organic substrate concentration to be consumed by the acidogenic bacteria x 2 (g/l) for growth and production of volatile fatty acids x 3 (mmol/l) (assumed to behave like pure acetate), x 4 (g/l) is the methanogenic bacteria concentration, x 5 (mmol/l) represents the inorganic carbon concentration and x 6 (mmol/l) the alkalinity concentration.

The related fed concentrations to the digester S 1in , S 2in , C in and Z in are supposed to be known and constant. The control inputs are u 1 = F 1in v (1/day) and u 2 = F 2in v (1/day), where F 1in is the waste feeding rate and F 2in is the input flow rate of the added alkalinity (Z ad ) to the digester. Since the later volume (v) is constant, then the output flow rate u = u 1 + u 2 (see Figure 1). The produced biogas is assumed to be composed of methane q m (x) and co 2 (q c (x)) gaseous. The later partial pressure is computed by P c (x). The rest of the used parameters in the model are defined in Table 1. The AD model ( 71)-(72) can be easily written under the form (4) using the following parameters

s = 1, n = 6, m = 2, n 1 = 2, n 2 = 2, n 3 = 4, ρ = u, A 0 = 0 n×n , A 1 = -block-diag(1, α, 1, α, 1, 1), G =   -k 1 1 k 2 0 k 4 0 0 0 -k 3 1 k 5 0 0 0 0 0 -1 0   T , γ(x) = µ 1 (x 1 )x 2 , µ 2 (x 3 )x 4 , q c (x) T , B = S 1in 0 S 2in 0 C in Z in 0 0 0 0 0 Z ad T , C =   1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1   , H 1 = 1 0 0 0 0 0 0 1 0 0 0 0 , H 2 = 0 0 1 0 0 0 0 0 0 1 0 0 , H 3 =     0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1     .

Simulation Results

In this section, we present a numerical example where we apply the proposed observer based tracking control to the AD process. The purpose is to track an admissible reference reflecting an enhanced quality of biogas. Moreover, to be as close as possible to real life experiments, we apply an additive zero-centered Gaussian noise to the performed measurements.

To run the simulation, we use the parameter values given in Table 1, and we take S 1in = 16 g/l, S 2in = 170 mmol/l, C in = 76.15 mmol/l, Z in = 200 mmol/l, Z ad = 700 mmol/l. Also, we initialize the system and the observer at x(0) = [1.8, 0.4, 12, 0.7, 53.48, 55] T and x0 = [1.8, 0.6, 12, 0.3, 45, 55] T , respectively. Regarding the desired reference, it corresponds to the steady state x d = [1.95, 0.60, 5.4, 1.38, 242.8, 240.34] T and the admissible reference input u d = [0.4966, 0.0436] T (we note that the pair (x d , u d ) satisfies the equations ( 4)). Besides, to solve the LMI conditions given by Theorem 3.2, we choose ρ min = 0.1 (1/day) and ρ max = 0.8 (1/day).

After solving the LMI conditions (25) and the optimization problem (47), which have been found feasible for β = 0.06, and by using the LMI Toolbox of Matlab and the solver SeDuMi, we have obtained the results depicted in Figures 2345678910. In the first 8 figures, the blue curve refers to the system state and the produced biogas quality without applying control. The red curve represents the desired reference, the green curve refers to the trajectory when applying control, and the black curve represents estimate of the controlled state. As it can be seen from these figures, despite the large initial estimation error and the corrupted measurements, the observer converges to the simulated system.

Moreover, we can see clearly from the results that despite the attached noise to measurements, the closed loop system tracks the desired trajectory. Indeed, we can see from Figure 8 that the biogas quality is enhanced by almost 10% when applying the designed control. This is a promising result from energy point of view of the process.

Last but not least, we want to emphasise that the controller behaviour, see Figures 9 and10, remains in an acceptable range of values. Indeed, in practice, for the well functioning of the digester and from the physical meaning of the state variables, variation of the control inputs is limited to avoid bacteria washout. This is, actually, what is obtained with the proposed observed based control strategy. Organic substrate concentration x1 (g/l) Remark 3. In a general situation, the upper and lower bounds of the parameter ρ u are given by calculating the minimum and maximum values of the parameter. It is well know that given a bounded function, it is always possible to calculate the minimum and maximum values of such a function. These values are fixed (and constant) and depend on the parameter. On the other hand, in the case of the considered anaerobic digestion model, the parameter ρ u depends on the input vector. Then, the values ρ min and ρ max are computed numerically by using the minimum and maximum values of the state vectors generated by the model. They can also depend on the saturated value of the control input if the saturation issue is considered. Consequently, the solutions of the LMIs depend on the saturated value. The LMIs may be feasible for a value and unfeasible for larger values (from the convexity principle). However, the technical analysis of this issue is one of the future work we aim to consider in a thorough way, from analytical and numerical points of view. We also aim to consider model uncertainties and noises with extended comparisons with other techniques from the literature.

Remark 4. Note that including a comparison of the results with another technique from the literature would be interesting and improve the paper. However, detailed comparisons are part of our future work. We aim to consider model uncertainties and process and measurement noises. On the other hand, some advantages of the proposed method compared to the MPC are obvious. MPC can work on thousands of different processes, especially for biogas. For instance, in [START_REF] Mauky | Model predictive control for demand-driven biogas production in full scale[END_REF], the authors proposed a model predictive control for demand-driven biogas production, which works successfully. Their method may perform quite as well as the proposed LMI-based method. Nevertheless, they incur higher computational effort, because the MPC requires online computations, and there are some other drawbacks when compared to the proposed LMI-based method:

• Although MPC works on several processes, it has some drawbacks related to the stability when we face general nonlinear systems. On the other hand, the proposed LMI-based techniques overcome some MPC drawbacks, and guarantee global convergence for nonlinear systems.

• The computational complexity of MPC for real-time applications is very high compared to the proposed methodology where the LMIs are solved offline.

• Finding the MPC parameters ensuring local convergence is not always easy. For each model, we need to proceed with different MPC parameters that may be difficult to find for some complex models. On the other hand, with our LMI design method, it is sufficient to write the model under the considered structure and compute the bounds of the partial derivatives of the nonlinearities to solve the LMIs.

Remark 5. It should be noticed that the proposed method does not solve all the AD process issues. Motivated by the selected AD model, the aim of the paper is to consider the problem of observer-based trajectory tracking for a class of nonlinear systems using an LMI framework. Hence, new LMI conditions are proposed, and the result is applied to an AD process. Indeed, the paper does not consider robustness and disturbance rejections and startups/shutdowns, which are important issues in real-life applications. The theoretical analysis of these issues is one of the future work we aim to consider in a deepen way, from analytical and numerical points of view, with extended comparisons with other techniques from the literature.

Conclusion

An observer-based reference tracking control for the AD process has been proposed in this paper. The control scheme is composed of an exponential nonlinear state observer and a feedback control. Actually, the state observer has been included in the control design to remedy the lack of some key measurements. Regarding the control law, feedback control has been used because it takes into account all dynamics of the process model.

In order to deal with the stability of the closed loop system, we have presented a kind of separation results. Indeed, relying on the triangular form of the system composed of dynamics of the estimation and tracking errors, we could establish the stability of each separately. Using an adequate reformulation of the Young's inequality and the Lipschitz property, we have provided the stability conditions in the form of enhanced end easily tractable LMIs. Firstly, to ensure the exponential stability of the estimation error. Then, to guarantee the H ∞ asymptotic stability of the tracking error.

We want to highlight that in order to extend the use of our technique and make it applicable for other systems belonging to the same class of systems as the AD process, we have presented the results in a general way. Furthermore, we have provided numerical simulations to illustrate the effectiveness of the proposed approach. In view of the simulation results, we target in the near future to extend the design methodology for saturation constraints on the control inputs. Moreover, as prospects, we will consider the disturbance rejections, startups/shutdowns, and parametric uncertainties, which are important issues in real-life applications. Extended comparisons to other available techniques in the literature would also increase the significance of the proposed design methodology.

Figure 1 :

 1 Figure 1: Controlled anaerobic digestion process.

Figure 3 :

 3 Figure 3: Acidogenic bacteria concentration x 2 (g/l).

Figure 4 :

 4 Figure 4: Acetate concentration x 3 (mmol/l).

Figure 5 :Figure 6 :Figure 7 :

 567 Figure 5: Methanogenic bacteria concentration x 4 (g/l).

Figure 10 :

 10 Figure 10: Control input u 2 (1/day).

  then we have

	+∞	+∞
	xT (t)x(t)dt ≤ µ	ωT (t)ω(t)dt + λ max P x0	2
	0	0

Table 1 :

 1 Model Parameters[START_REF] Bernard | Dynamical model development and parameter identification for an anaerobic wastewater treatment process[END_REF]

	Acronyms Definition	Units	Value
	α	Proportion of dilution rate for bacteria	mmol/l	0.5
	k 1	Yield for substrate (x 1 ) degradation	g/(g of x 2 )	42.1
	k 2	Yield for VFA (x 3 ) production	mmol/(g of x 2 ) 116.5
	k 3	Yield for VFA consumption	mmol/(g of x 4 ) 268
	k 4	Yield for co 2 production	mmol/g	50.6
	k 5	Yield for co 2 production	mmol/g	343.6
	k 6	Yield for ch 4 production	mmol/g	453
	µ 1	Maximum acidogenic bacteria (x 2 ) growth rate	1/day	1.25
	µ 2	Maximum methanogenic bacteria (x 4 ) growth rate 1/day	0.74
	k s1	Half saturation constant associated with x 1	g/l	7.1
	k s2	Half saturation constant associated with x 3	mmol/l	9.28
	k i2	Inhibition constant associated with x 3	mmol/l	256
	k b	Acidity constant of bicarbonate	mol/l	6.5 10 -7
	K H	Henry's constant	mmol/(l.atm)	27
	P T	Total preasure	atm	1.013
	k La	Liquid/gas transfer constant	1/day	19.8
	4.2. Correspondance between (4) and (71)-(72)

  Figure 2: Organic substrate concentration x 1 (g/l).

		3.5
																			x 1
		3				1.98 2 2.02													x 1c x 1d
						1.92 1.94 1.96													x1
		2.5					19.95	20	20.05	20.1	20.15	20.2	20.25	20.3	20.35	20.4
		2
		1.5
	Acidogenic bacteria concentration (g/l)	0.35 0.4 0.45 0.55 0.6 0.65 0.7 0.5	0	5	0.605 0.6055 0.606 0.6065	27.8	28	10	28.2	28.4		15 28.6	28.8		29	20 29.2	29.4	25	x 2d x 2c x2 x 2	30
										Time (days)


Karim Benhenni

Ali Hajj Hassan

Yingcai Su

Local polynomial estimation of regression operators from functional data with correlated errors

Keywords: El Niño time series, Fréchet derivatives, Functional fixed-design data, Local polynomial estimation, Nonparametric regression operator, Short and long memory processes

This article considers the problem of nonparametric estimation of the regression operator r in a functional regression model Y = r(x) + ε with a scalar response Y, a functional explanatory variable x, and a second order stationary error process ε. We construct a local polynomial estimator of r together with its Fréchet derivatives from functional data with correlated errors. The convergence in mean squared error of the constructed estimator is studied for both short and long range dependent error processes. Simulation studies on the performance of the proposed estimator are conducted, and applications to independent data and El Niño time series data are given.

Introduction

Data in the form of a collection of continuous function curves become increasingly common in practice. Accordingly, the interest in functional data analysis (FDA) has been growing steadily over the past few decades. The books [START_REF] Bosq | Linear Processes in Function Spaces: Theory and Application[END_REF][START_REF] Ferraty | Nonparametric Analysis for Functional Data: Theory and Practice[END_REF][START_REF] Ramsay | Applied Functional Data Analysis: Methods and Case Studies[END_REF][START_REF] Ramsay | Functional Data Analysis[END_REF], and the overview paper [START_REF] Cuevas | A partial overview of the theory of statistics with functional data[END_REF] provide a general coverage of FDA theory and methodology. More recent contributions to this field can be found in the book [START_REF] Aneiros | Contributions to Statistics: Functional Statistics and Related Fields[END_REF]. FDA-based statistical techniques find applications in a wide range of areas such as environmental, medical, meteorological sciences and statistical control processes [START_REF] Sancho | Variability analysis by statistical control process and functional data analysis: Case of study applied to power system harmonics assessment[END_REF]. Interestingly, FDA techniques were even applied to problems of radon entry rate estimation [START_REF] Brabec | Simplified radon entry rate estimation methodology from tracer and radon concentrations based on functional data analysis approach[END_REF] and identification of abnormal nitrogen oxides emission levels [START_REF] Febrero | Outlier detection in functional data by depth measures with application to identify abnormal NO x levels[END_REF].

In the finite-dimensional framework, there is a vast literature on local polynomial regression estimation under the assumption of independent random errors. In this regard, asymptotic expressions for the bias and variance of the corresponding estimators can be found in [START_REF] Fan | Local Polynomial Modeling and Its Applications[END_REF] and [START_REF] Wand | Kernel Smoothing[END_REF]. In the case of regression estimation with dependent errors, local polynomial estimation is studied under mixing conditions in [START_REF] Masry | Local polynomial estimation of regression functions for mixing processes[END_REF], association in [START_REF] Masry | Local polynomial regression estimation with correlated errors[END_REF], and (stationary) short range dependence in [START_REF] Francisco-Fernández | Local polynomial regression estimation with correlated errors[END_REF] and [START_REF] Pérez-González | Asymptotic properties of local polynomial regression with missing data and correlated errors[END_REF]. For more general stationary error processes (correlated errors) and local polynomial estimators based on repeated independent observations, the rates of convergence of mean squared errors (MSE) and asymptotic distributions were obtained in [START_REF] Benhenni | Local polynomial estimation of the mean function and its derivatives based on functional data and regular designs[END_REF].

In the infinite-dimensional framework, in the case of dependent data and random functional designs, consistent estimators of regression operators with scalar responses under certain mixing conditions were obtained in [START_REF] Ferraty | Nonparametric analysis for functional data, with application in regression method[END_REF]. There are many practical situations in which functional designs are fixed and the errors are correlated. For functional fixed designs, the asymptotic properties of the Nadaraya-Watson kernel estimator in terms of MSE consistency and asymptotic normality were obtained in [START_REF] Benhenni | Estimation of the regression operator from functional fixed-design with correlated errors[END_REF]. Local linear regression for functional data and scalar response with independent errors was introduced in [START_REF] Baillo | Locally modeled regression and functional data[END_REF] and [START_REF] Boj | Distance-based local linear regression for functional predictor[END_REF]. The MSE consistency of the local linear regression estimator when the covariate takes values in a Hilbert space, namely L 2 [0, 1], was derived in [START_REF] Baillo | Locally modeled regression and functional data[END_REF]. An easier computational expression for the local estimator was constructed and its almost complete convergence, along with the corresponding rates, was shown in [START_REF] Barrientos-Marin | Local linear regression for functional predictor and scalar response[END_REF]. A more general approach for constructing the local linear estimator was taken in [START_REF] Berlinet | Local linear regression for functional data[END_REF], where a bound on the MSE was also given. Other related problems such as conditional quantile regression or conditional distribution estimation based on local linear methods can be found in [START_REF] Rachdi | Theoretical and practical aspects of the quadratic error in the local linear estimation of the conditional density for functional data[END_REF] and [START_REF] Zoulikha | Functional quantile regression: Local linear modelisation[END_REF].

Here, we extend these results to the local polynomial estimation framework and study the effect of short and long memory error processes, as defined in [START_REF] Beran | Statistics for Long-Memory Processes[END_REF], on the rate of convergence of the MSE. More precisely, we construct a truncated local polynomial estimator of the regression operator together with its Fréchet derivatives based on functional data with correlated errors, and then establish the asymptotic behavior, including the rate of convergence, of the corresponding MSE. The procedure proposed in this article for the estimation of Fréchet derivatives of the regression operator is nonparametric. To the best of our knowledge, the only earlier work involving estimation of functional derivatives is due to Hall et al. [START_REF] Hall | Estimation of functional derivatives[END_REF]. In the construction and analysis of the proposed estimators, they used the first order derivative of the regression operator taken at the orthonormal eigenfunctions of the estimated covariance kernel.

The nonparametric k-nearest-neighbor (kNN) smoothing technique was used by Kara et al. [START_REF] Kara | Data-driven kNN estimation in nonparametric functional data analysis[END_REF] in the estimation of a regression operator for iid samples. In their approach, the value of the parameter k is fixed but the smoothing parameter, the bandwidth, is data dependent. It is of interest to see how the kNN method works out in the setting of dependent samples. A nonparametric approach with data-driven smoothing parameters to the regression operator estimation was also investigated in [START_REF] Chagny | Adaptive estimation in the functional nonparametric regression model[END_REF] for iid errors. It seems worthwhile to investigate possible extensions of the approach of this article to a semi-parametric set-up, functional single index model as investigated in [START_REF] Ait-Saïdi | Cross-validated estimation in the single-function index model[END_REF] and [START_REF] Goia | A partitioned single functional index model[END_REF], and also to multiple index functional regression models as in [START_REF] Chen | Single and multiple index function regression models with nonparametric link[END_REF].

In nonparametric time series modeling, to overcome the difficulty of using an overwhelmingly high number of past observations, Aneiros and Vieu [START_REF] Aneiros | Nonparametric time series prediction: A semi-functional partial linear modeling[END_REF][START_REF] Aneiros | Partial linear modeling with multi-functional covariates[END_REF] proposed a functional approach, namely the semi-functional partial linear regression (SFPLR) model. The work reported in this article may be extended to their semi-parametric situations, when the response, explanatory variables, dependence of data, and covariates are properly identified or introduced. That being said, SFPLR modeling is interesting in its own right and certainly a promising approach for the myriad of real life functional time series data. The results of this article are potentially useful for linear and nonlinear inverse problems in geophysical sciences [START_REF] Mcgillivray | Methods for calculating Fréchet derivatives and sensitivities for the non-linear inverse problem: A comparative study[END_REF] as well. In our discussion and analysis, short and long memory dependent errors are assumed. It is worth noting that the idea of functional canonical correlation coefficients as in [START_REF] He | Functional canonical analysis for square integrable stochastic processes[END_REF] provides another way to characterize the dependence structure between and within functional samples. All of these possible extensions or further studies are beyond the scope of this article and thus will not be discussed further. However, readers are referred to [START_REF] Geenens | A Nonparametric Functional Method for Signature Recognition[END_REF] and [START_REF] Ling | Nonparametric modeling for functional data: Selected survey and tracks for future[END_REF] for recent surveys on the topic of the nonparametric FDA.

The organization of this paper is as follows. In Section 2, we state the functional regression model and define the aforementioned nonparametric estimator. In Section 3, we set up the appropriate assumptions and give the MSE convergence of the local estimator of the regression operator and its Fréchet derivatives along with the corresponding rates. Sections 4 and 5 are devoted to various simulation studies for both short and long memory error processes and also to an application to spectrometric data and El Niño time series data. Finally, a summary and some further comments are given in Section 6, and the proofs of the main results are presented in the Appendix. Moreover, the codes and data for reproducing the simulations and data analysis are included in the Online Supplement.

Estimation of the regression operator and its Fréchet derivatives

We consider a regression model with regression operator r defined, for all i ∈ {1, . . . , n}, by

Y i = r(x i ) + ε i ,
where Y 1 , . . . , Y n are observations of the response variable Y, x 1 , . . . , x n the explanatory data of functional fixed-design type, and ε 1 , . . . , ε n correlated random errors. The response Y is a real-valued random variable such that E |Y| < ∞. The explanatory variable x is assumed to belong to a metric space (H, d) equipped with the metric d. The random error ε is assumed to be a centered second order stationary process with autocovariance function defined, for all i, j ∈ Z, by E(ε i ε j ) = ρ ε (ij), and unit variance, viz. E(ε 2 i ) = ρ ε (0) = 1. In this section, a local polynomial estimator of the regression operator r, together with its Fréchet derivatives, is proposed. In the next section, we investigate its asymptotic performance for both short and long range dependence (or memory) error processes. For this, we recall briefly the definitions of short and long memory processes.

For some stochastic processes, such as the popular autoregressive and moving average models, the autocovariance function decays exponentially and is summable, viz.

+∞ τ=-∞ |ρ ε (τ)| < ∞.
These random processes are called short range dependent (SRD) or short range correlated stochastic processes. Some standard time series models, such as the popular Fractional ARIMA processes or Fractional Difference processes, assume the non-summability of the autocovariance sequence {ρ ε (i) : i ∈ Z}, which captures the intuition underlying long range dependence (LRD) or long range memory, viz.

+∞ τ=-∞ |ρ ε (τ)| = ∞.
For instance, a stationary process {ε t : t ∈ Z} is a long memory process if for some real γ ∈ (0, 1], and some positive constant C, its autocovariance function satisfies, as j → ∞,

|ρ ε ( j)| ∼ C | j| -γ ;
(1) see [START_REF] Beran | Statistics for Long-Memory Processes[END_REF] for more details.

Let U be an open subset of H.

If f : U → R is p + 1 times differentiable on U, the Taylor expansion of f around x ∈ U is f (y) = p j=0 1 j! D j f (x) d j (x, y) + R p (y),
where D j f (x) denote the jth Fréchet derivative of f at x, with x and y being points of U and d the metric of H. The remainder is given by

R p (y) = 1 p! 1 0 D p+1 f {x + td(x, y)} d p+1 (x, y)(1 -t) p dt = o{d p (x, y)}
as d(x, y) → 0; see [START_REF] Cartan | Cours de calcul différentiel[END_REF] and [START_REF] Zeidler | Nonlinear Functional Analysis and Its Applications I: Fixed-Point Theorems[END_REF] for more details.

When the observations Y 1 , . . . , Y n are possibly dependent, our goal is to estimate the regression operator r and its Fréchet derivatives r j = D j r. To this end, we introduce a functional truncated version of the local polynomial kernel estimate. The estimators of r 0 (x), . . . , r p (x) based on the observations Y 1 , . . . , Y n and the functional fixed design x 1 , . . . , x n are obtained for x ∈ U, by minimizing the criterion

n i=1          Y i - p j=0 β j (x)d j (x i , x)          2 K h {d(x i , x)} with K h (•) = h -1 K(•/h
) and β j (x) = r j (x)/ j!, for all j ∈ {0, . . . , p}, where the kernel K is a real-valued function defined on R + with a compact support and h = h(n) is the bandwidth such that h ∈ R + and h → 0 as n → ∞.

Let Y n = (Y 1 , . . . , Y n ) and denote the canonical basis of R p+1 by (e k ) p k=0 , i.e., the (k + 1)st component of the column vector e k is 1 and the rest of the components are zeros. Define the n × (p + 1) matrix X n by

X n =             1 d 1 (x 1 , x) . . . d p (x 1 , x) . . . . . . . . . . . . 1 d 1 (x n , x) . . . d p (x n , x)            
and the diagonal matrix

W n = diag [K h {d(x 1 , x)}, . . . , K h {d(x n , x)}].
Then the estimators r j (x) are given, for all j ∈ {0, . . . , p}, by

r j (x) = j! e j β n (x), (2)
where

β n (x) = (X n W n X n ) -1 X n W n Y n .
The truncated functional Nadaraya-Watson estimator of the regression operator r can be obtained by taking j = 0 and p = 0, viz.

r NW (x) =

x i ∈B(x,h) W i (x) Y i
where B(x, h) = {y : d(x, y) ≤ h} is the ball with center x and radius h in the metric space (H, d), and

W i (x) = ∆ i (x)/ x i ∈B(x,h) ∆ i (x) with ∆ i (x) = K h {d(x, x i )}.

Asymptotic performance

Notations and assumptions

The following notations and assumptions are required for the statements of the main results. The kernel K is a function with compact support on [0, 1] and has a continuous derivative in (0, 1) such that K(1) > 0 and K (s) ≤ 0.

(

) 3
Thus the kernel K is a decreasing function that gives more weight to the argument close to zero. For all j ∈ {0, . . . , p} and a bandwidth h > 0, let

ϕ j (h) = 1 n n i=1 d j (x i , x) h j 1 B(x,h) (x i )
where 1 D (x) stands for the indicator function of a set D in H, and for all k, ∈ {0, . . . , p} and s, s ∈ (0, 1],

µ k, (hs, hs ) = 1 n 2 1≤i, j≤n i j d k (x i , x) (hs) k d (x j , x) (hs ) ρ (i -j)1 B(x,hs) (x i )1 B(x,hs ) (x j ).
In the functional random design, the small ball probabilities are determined by the distribution function F x (t) = Pr{d(X, x) ≤ t}, and E{d j (X, x)} represents the jth moment of the random variable d(X, x). In the functional fixed design, this is equivalent to the empirical jth moment frequency ϕ j (t). Likewise, µ k, (t, s) represents the empirical (k, )th joint weighted moment frequency. Then we assume that ϕ j and µ k, satisfy a certain regular variation condition at the origin, i.e., as h → 0,

τ h j (s) = ϕ j (hs)/ϕ j (h) → τ j (s) (4)
and

ψ h k, (s, s ) = µ k, (hs, hs )/µ k, (h, h) → ψ k, (s, s ). ( 5
)
When Fréchet derivatives of the regression operator are estimated, assumptions (4) and ( 5) are natural extensions of the assumptions considered in [START_REF] Benhenni | Estimation of the regression operator from functional fixed-design with correlated errors[END_REF] and [START_REF] Ferraty | Nonparametric regression on functional data: Inference and practical aspects[END_REF], which correspond to the case j = 0 of this article. Finally, we list the following constants that are related to K, the τ j s, and ψs; they will be used in the next section. For ν ∈ {1, 2},

C ν j = K ν (1) - 1 0 τ j (s)s j (K ν ) (s)ds, C 1,1 k, = K 2 (1) -K(1) 1 0 ψ k, (s, 1)K (s)ds -K(1) 1 0 ψ k, (1, s)K (s)ds + 1 0 1 0 ψ k, (s, t)K (s)K (t)dsdt. (6)

Asymptotic results

For the asymptotic expressions of the bias and the variance of the local polynomial estimators, we have the following results.

Theorem 1. Assume that the regression operator r is p + 1 times continuously differentiable around x. If assumptions (3) and (4) are satisfied, then for all j ∈ {0, . . . , p}, we have

E{r j (x)} -r j (x) = j! r p+1 (x) (p + 1)! h p+1-j e j Σ -1 h μh {1 + o(1)},
where

Σ h = (ϕ k+ (h)C 1 k+ ) k, ∈{0,...,p} is a (p + 1) × (p + 1) matrix, μh = (ϕ p+1 (h)C 1 p+1 , . . . , ϕ 2p+1 (h)C 1 2p+1
) is a column vector of size p + 1, and the constants C 1 j are as in (6).

Theorem 2. Assume that the regression operator r is p + 1 times continuously differentiable around x. If assumptions (3), (4), and (5) are satisfied, then for all j ∈ {0, . . . , p}, we have

var{ r j (x)} = 1 n e j H -1 Σ -1 h (Σ (2) h + nΣ (1,1) h )Σ -1 h H -1 e j {1 + o(1)},
where

H = diag(1, h, . . . , h p ), the matrices Σ (2) h = (ϕ k+ (h)C 2 k+ ) k, ∈{0,...,p} and Σ (1,1) h = (µ k, (h, h)C 1,1 k, ) k, ∈{0,.
..,p} are of size (p + 1) × (p + 1), and the constants C 2 j and C 1,1 k, are as in [START_REF] Barrientos-Marin | Local linear regression for functional predictor and scalar response[END_REF]. We can obtain simpler expressions of the asymptotic bias and variance terms under some assumptions on the empirical frequencies.

Corollary 1. Assume that there exist some positive constants α 0 , . . . , α p and γ k, , for all k, ∈ {0, . . . , p} such that, as h → 0, ϕ j (h)/ϕ 0 (h) → α j and µ k, (h, h)/µ 0,0 (h, h) → γ k, , then under assumptions (3), (4), and (5), we have

E{r j (x)} -r j (x) = j! r p+1 (x) (p + 1)! h p+1-j e j M -1 µ{1 + o(1)} and var{ r j (x)} = 1 nh 2 j ϕ 0 (h) e j M -1 M (2) M -1 e j {1 + o(1)} + µ 0,0 (h, h) h 2 j ϕ 2 0 (h) e j M -1 M (1,1) M -1 e j {1 + o(1)}, where µ = (α p+1 C 1 p+1 , . . . , α 2p+1 C 1 2p+1 ) is a column vector of size p + 1, M = (α k+ C 1 k+ ) k, ∈{0,...,p} , M (2) = (α k+ C 2 k+ ) k, ∈{0,...,p} , and M (1,1) = (γ k, C 1,1 k, ) k, ∈{0,.
..,p} are matrices of size (p + 1) × (p + 1). Remark 1. The asymptotic MSE of the Nadaraya-Watson estimator rNW can be deduced from Theorems 1 and 2 by taking j = 0 and p = 0 (local constant)

MSE(r NW ) = O(h 2 ) + 1 nϕ 0 (h) C 2 0 (C 1 0 ) 2 + µ 0,0 (h, h) ϕ 2 0 (h) C 1,1 0,0 (C 1 0 ) 2 {1 + o(1)},
where ϕ 0 (h) is the empirical frequency of curves belonging to the closed ball of center x and radius h. In this case, the dominant term in the variance is of order O{µ 0,0 (h)/ϕ 2 0 (h)} for a long memory error process, and O[1/{nϕ 0 (h)}] for a short memory error process.

We also deduce from Theorems 1, 2 and Corollary 1, the corresponding rates of convergence of the local estimator of the regression operator and its Fréchet derivatives when the errors are iid.

Corollary 2. If Y 1 , . . . , Y n (or ε 1 , . . . , ε n ) are independent and identically distributed, then under assumptions (3), (4), and (5), and if nh 2 j ϕ 0 (h) → ∞ as n → ∞, we have

MSE{ r j (x)} = O(h 2(p+1-j) ) + 1 nh 2 j ϕ 0 (h) e j M -1 M (2) M -1 e j {1 + o(1)}.
In the case of long and short memory error processes, we deduce the following rates of convergence of the local estimator of r and its Fréchet derivatives. Corollary 3. In the case of a long memory, under Assumptions (1), (3), (4), and (5) we have that if

{nϕ 0 (h)} γ h 2 j → ∞ as n → ∞, then MSE{ r j (x)} = O(h 2(p+1-j) ) + O 1 {nϕ 0 (h)} γ h 2 j .
Corollary 4. In the case of a short memory, under Assumptions (3), (4), and (5)

we have that if nϕ 0 (h)h 2 j → ∞ as n → ∞, then MSE{ r j (x)} = O(h 2(p+1-j) ) + O 1 nϕ 0 (h)h 2 j .
It should be noted that, as in the finite-dimensional framework, the performance in terms of the rate of convergence of the MSE declines when higher Fréchet derivatives of the regression operator are estimated. The bias term is of smaller order for higher-order polynomial approximation. The rate of convergence of the variance does not depend on the order p and hence can only be affected by the order of approximation through the asymptotic constant.

It can also be observed from these results that if the local polynomials are used for estimating any order Fréchet derivative r 1 , . . . , r p , the rate of convergence coming from the variance component is lower in the case of the long memory error process than in the case of the short memory error process.

Simulation studies

This section gives some simulations which show the advantage of using the local polynomial estimates of the regression operator r in the case of the functional fixed design model with correlated errors. To this end, we simulate two types of smooth curves: parabolic curves and sinusoidal curves. A sample of 100 curves of each type is plotted in Figure 1.

For both cases, we generate different curves of sample size n = 200; each curve is discretized with 100 equidistant values in [0, 1]. The response variable is obtained from the model defined, for all k ∈ {1, . . . , n}, by

y k = r(x k ) + k .
Recall that the error process is a second order stationary process. As a short memory process, we consider an Ornstein-Uhlenbeck process (OU) which satisfies the stochastic differential equation

d t = (θ 1 -θ 2 t ) dt + θ 3 dW t ,
where θ 1 , θ 2 and θ 3 are three deterministic parameters and W t denotes the Wiener process. Here, we take θ 1 = 0, θ 2 = 0.75 and θ 3 = 0.5. For a long memory error process, we use an autoregressive fractionally integrated moving average ARFIMA(p, d, q) with parameters p = q = 0 and fractional differencing parameter d = 0. The three estimators are constructed from the Epanechnikov kernel. The explanatory variable x belongs to the metric space L 2 [0, 1] equipped with the metric d: for any two curves x k and x k ,

d(x k , x k ) = 1 0 {x k (t) -x k (t)} 2 dt.
The selection of a proper smoothing parameter (bandwidth) is a crucial step for the estimation of the regression operator. We use a local bandwidth choice which is known to outperform the global one; see [START_REF] Benhenni | Local smoothing regression with functional data[END_REF]. The method for selecting the bandwidth h at the curve x consists of choosing h x to minimize, over a particular subset H n of smoothing parameters to be specified, the local cross-validation criterion

LCV x (h) = 1 n n i=1 {Y i -r i h (x i )} 2 W n,x (x i ),
where r i h is the leave-one-curve-out estimate of r. The locally selected bandwidth is

h x = arg min h∈H n LCV x (h).
We take, for all i ∈ {1, . . . , n}, the local weights

W n,x (x i ) = 1 if d(x i , x) < h, 0 otherwise.
In practice, one replaces the minimization of LCV x over the infinite set (0, ∞) (i.e., h > 0) with a minimization problem over a finite set. To this end, we consider the local bandwidth h k (x) such that the ball centered at x of radius h k (x) contains exactly k neighbors (i.e., k curves around x). As a direct consequence, we have for all x ∈ H, W 1,x (x i ) + • • • + W n,x (x i ) = k when h = h k (x). In this way, we can look for the optimal local bandwidth among those which consider respectively one neighbor, two neighbors, etc. The optimal retained number k opt of neighbors allows to build the optimal local bandwidths h k opt (x) at any x ∈ H. Finally, the performance of an estimator r is evaluated by computing the relative mean squared error, viz.

RMSE(r) = n k=1 {y k -r(x k )} 2 / n k=1 (y k -y) 2 .

Parabolic curves and linear operator

We first generate 200 parabolic smooth curves by setting, for all k ∈ {1, . . . , 200} and t ∈ [0, 1], x k (t) = kt 2 . We take the following linear regression operator

r(x) = 1 0 x(t) dt.
Figure 2 gives the plot of the y k s versus the r(x k )s for the two error processes when SNR = 10. The simulation results are summarized in Table 1. The values of RMSE of each estimator increase when the SNR decreases. It appears that our estimators (linear and quadratic) outperforms the Nadaraya-Watson estimator regardless of the type of error process and the value of SNR that are being used. For instance, when SNR = 10, there is a reduction of the RMSE, 26% for OU and 33% for ARFIMA (respectively 56% for OU and 54% for ARFIMA), of the local linear estimator (respectively the local quadratic estimator) over the NW estimator. Lastly, Figures 34give a graphical comparison by displaying the estimated responses in terms of the observed responses for the two types of error processes when SNR = 10.

Sinusoidal curves and nonlinear operator

The second simulation study consists in creating a set of 200 sinusoidal smooth curves by setting, for all k ∈ {1, . . . , 200} and t ∈ [0, 1], x k (t) = √ k sin{4(kt)}. We take the following nonlinear regression operator:

r(x) = 1 0 dt 1 + x 2 (t)
.

Figure 5 gives the plot of the y k s versus the r(x k )s for the two error processes. Once again, the simulation results summarized in Table 2 show the superiority of the local linear and quadratic estimators over the Nadaraya-Watson estimator. For instance, when SNR = 10 there is a significant reduction of the RMSE, 81% for OU and 51% for ARFIMA (respectively 90% for OU and 57% for ARFIMA), of the local linear estimator (respectively the local quadratic estimator) over the NW estimator. As in the parabolic curves situation, the performance of the estimators is affected by the nature of the error process. To illustrate this point, we displayed in Figures 67the observed responses in terms of their estimations when SNR = 10. Moreover, it seems from these figures that the linear and the quadratic estimators are more robust with respect to some high values (outliers) than the classical NW estimator. Finally, The shape of the curves has an impact on the quality of the estimation. The estimates obtained are less efficient than those of parabolic curves.

Application to real data

This section aims to confirm the superiority of the local linear and quadratic estimators over the (local constant) Nadaraya-Watson one using two well-known functional data sets: the spectrometric curves (case of independent errors) and El Niño time series (case of dependent errors).

Spectrometric data

For the ith unit among 215 pieces of finely chopped meat, we observe one spectrometric curve x i which corresponds to the absorbency measured at 100 wavelengths, and its fat content y i is obtained by an analytical chemical process.

We also compare the performance of the three estimators by computing the relative mean squared errors of the 215 estimated values ŷi . The corresponding results are shown in Table 3. Furthermore, we compare the estimation performance of the three estimators constructed from the semi-metric based on the second derivative because of the smoothness of the spectrometric curves; see Figure 8.

When the L 2 metric is used, there is an important gain of 74% (respectively 68%) in RMSE of linear estimator (respectively quadratic estimator) compared to the NW one. The use of the semi-metric based on the second derivative significantly improves the quality of the estimation of each estimator, but the linear and quadratic estimators outperform the NW estimator. Moreover, it is seen that for this data set, there is no need to use a higher-order approximation since the local linear estimator already has a very good performance. All these facts are illustrated in Figure 9, where the estimated fat values of each estimator are plotted versus the observed ones with respectively L 2 metric and second derivative semi-metric.

El Niño data

The second application concerns the El Niño time series which gives monthly sea surface temperatures over 54 years, from June, 1950 to May, 2004. This series was considered by [START_REF] Ferraty | Conditional quantiles for dependent functional data with application to the climatic El Niño phenomenon[END_REF] for dependent functional data. The data are recorded as a sequence of 648 real values as follows. Let Z(1), . . . , Z(648) be the El Niño time series. Then as in [START_REF] Ferraty | Conditional quantiles for dependent functional data with application to the climatic El Niño phenomenon[END_REF], we can construct subsequences by setting z i (t) = Z{12(i -1) + t} for all t ∈ {1, . . . , 12} and i ∈ {1, . . . , 54}. Because the climatic phenomenon is changing continuously over time, it is natural to consider each annual curve as a continuous curve {z i (t) : t ∈ [0; 12)} observed at 12 discretized points. Therefore, this time series can be viewed as a sample of 54 dependent functional data.

In order to study the performance of our estimators we will ignore the 54th year and will predict it from the 53 previous ones. That is, for all i ∈ {1, . . . , 52} and any fixed δ ∈ {1, . . . , 12}, we take Y i (δ) = z i+1 (δ). In this way, we build a sample of 52 dependent pairs (Y 1 (δ), z 1 ), . . . , (Y 52 (δ), z 52 ) where the Y i (δ)s are real-valued random variables and the z i s are functional data.

Our main goal here is to predict, for each of the 12 months, the sea surface temperature of the 54th year by means of the Niño data of the 53 previous years. The corresponding estimates based on L2 metric are numerically illustrated by computing the mean squared error While the performance of the local linear estimate is good, there is no real improvement of these local polynomial estimates over the classical estimates for prediction purposes. Moreover, it is seen for this data set that the higherorder approximation does not bring a major improvement over the local constant or linear estimates. This may be due to the fact that the bandwidths used in the three local estimates, local constant (p = 0), linear (p = 1) and quadratic (p = 2) are very different and this has a significant effect on the bias and the variance. Besides, the local polynomial estimate depends very much on the inverse of a matrix and this may be subject to high numerical instability when the sample size n and the order p increase.

Other semi-metrics could be used, such as the one based on principal component analysis. However, the main goal here is not really the choice of the semi-metric but a comparison in terms of performance of the three estimators corresponding to p ∈ {0, 1, 2}.

Further comments

We established some asymptotic results for estimating the regression operator and its Fréchet derivatives from a polynomial approximation of any order in a Banach space equipped with a metric. We showed how these asymptotic performances are affected by the order of the local polynomial estimation. When the curves live in a semi-metric space with a semi-metric based on the derivatives of the curves, namely for every integer ≥ 2,

x ( ) k (t) -x ( ) k (t) and C 1,1 k, (h) = K 2 (1) -K(1) 1 0 ψ k, (hs, h)K (s)ds -K(1) 1 0 ψ k, (h, hs)K (s)ds + 1 0 1 0 ψ k, (hs, ht)K (s)K (t)dsdt → C 1,1 k, ,
where C ν j and C 1,1 k, are as in [START_REF] Barrientos-Marin | Local linear regression for functional predictor and scalar response[END_REF].

Proof of Theorem 1. Write r n = (r(x 1 ), . . . , r(x n )) and define the (p + 1) × (p + 1) matrix S n = X n W n X n with (k, )th element given by

s k+ = n i=1 d k+ (x i , x)K h {d(x i , x)}.
It follows from the definition of the estimator (2) that E{ β n (x)} = S -1 n X n W n r n . The Taylor expansion of r(x i ) up to order p + 1 yields

r(x i ) = r(x) + d(x i , x)r 1 (x) + • • • + d p+1 (x i , x) (p + 1)! r p+1 (x) + o{d p+1 (x i , x)}
and thus

r n = X n β(x) + {β p+1 + o(1)}             d p+1 (x 1 , x) . . . d p+1 (x n , x)             .
As a result, the bias in the estimation of β(x) is

E{ β n (x)} -β(x) = {1 + o(1)}β p+1 S -1 n c n ,
where c n = (s p+1 , . . . , s 2p+1 ) . We approximate the elements of the matrix S n as follows.

For each i ∈ D x (h), we have

K h {d(x i , x)} = K(1)1 B(x,h) (x i ) - 1 0 1 B(x,hs) (x i )K (s)ds.
Then

s k+ = K(1)nh k+ ϕ k+ (h) -nh k+ 1 0 ϕ k+ (hs)s k+ K (s)ds = nh k+ ϕ k+ (h) K(1) - 1 0 τ h k+ (s)s k+ K (s)ds .
From assumption (4) with ν = 1, we then have, for all k, ∈ {0, . . . , p},

s k+ = nh k+ ϕ k+ (h)C 1 k+ (h){1 + o(1)}. It follows that S n = n{1 + o(1)}HΣ h H. Now, c n =             s p+1 . . . s 2p+1             =              nh p+1 ϕ p+1 (h)C 1 p+1 (h){1 + o(1)} . . . nh 2p+1 ϕ 2p+1 (h)C 1 2p+1 (h){1 + o(1)}              = nh p+1 H μh {1 + o(1)}, so that S -1 n c n = h p+1 H -1 Σ -1 h μh {1 + o(1)} and hence E{ β n (x)} -β(x) = β p+1 h p+1 H -1 Σ -1 h μh {1 + o(1)},
from which the asymptotic bias expression of r j (x), for all j ∈ {0, . . . , p}, in Theorem 1 is obtained.

Proof of Theorem 2. Define the n × n matrix V n = {ρ ε (ij)} and the (p + 1)

× (p + 1) matrix S * n = X n W n V n W n X n . Noting that var(Y n ) = V n , we can see from the definition of the estimator (2) that var{ β n (x)} = S -1 n S * n S -1 n ,
where the (k, )th element of the matrix S * n is

s * k, = n i=1 n j=1 d k (x i , x)d (x j , x)∆ i (x)∆ j (x)ρ (i -j) = n i=1 d k+ (x i , x)∆ 2 i (x)ρ (0) + 1≤i, j≤n i j d k (x i , x)d (x j , x)∆ i (x)∆ j (x)ρ (i -j).
For the first term, following the same steps as for deriving the asymptotic approximation of the elements s k+ of the matrix S n , we obtain

n i=1 d k+ (x i , x)∆ 2 i (x)ρ (0) = nh k+ ϕ k+ (h)C 2 k+ (h)ρ (0){1 + o(1)}.
For the second term, let D x (h) = {i : d(x, x i ) ≤ h} = {i : x i ∈ B(x, h)}. Then using the assumption (5), we deduce that 1≤i, j≤n i j Given that D x (hs) ⊂ D x (h) for all s ∈ [0, 1], we have D x (h) ∩ D x (hs) = D x (hs), it follows that 1≤i, j≤n i j d k (x i , x)d (x j , x)∆ i (x)∆ j (x)ρ (ij) = K 2 (1)n 2 h k+ µ k, (h, h) -K(1)n 2 h k+

d k (x i , x)d (x j , x)∆ i (x)∆ j (x)ρ (i -j) = K 2 (1) i, j∈Dx (h)
= 1 n H -1 Σ -1 h H -1 nHΣ (2)
h Hρ (0) + n 2 HΣ (1,1) h H

1 n H -1 Σ -1 h H -1 {1 + o(1)} = 1 n H -1 Σ -1 h Σ (2)
h ρ (0) + nΣ (1,1) h

Σ -1 h H -1 {1 + o(1)}.
This concludes the argument.

Proof of Corollary 3. Let n x (h) = card {i : x i ∈ B(x, h)} = card {D x (h)}, so that ϕ 0 (h) = n x (h)/n. Then under assumptions (1), ( 3), (4), and (5), we have, for some constant C, Given that the error process is of short memory, we have n k=1 |ρ (k)| = O(1). It then follows that µ 0,0 (h, h) = O{n x (h)/n 2 } = O{n -1 ϕ 0 (h)}.

µ 0,0 (h, h) = 1 n 2 1≤i, j≤n i j |ρ (i -j)| 1 B(x,h) (x i )1 B(x,h) (x j ) ≤ 1 n 2 i, j∈Dx (h) i j C |i -j| -γ ≤ 2C n 2 i, j∈Dx (h) i< j |i -j| -γ ≤ 2C n 2 n x (h) k=1 {n x (h) -k}k -γ ≤ 2Cn x (h) 2-γ n 2 × 1 n x (h) n x (h) k=1 1 - k n x (h) k n x (h) -γ n x (h) 2-γ
Therefore, var{ r j (x)} = O 1 nh 2 j ϕ 0 (h) , which concludes the argument.

3 .σ 2 r

 32 The standard deviation σ of the error process is controlled by considering various values of the Signal-to-Noise Ratio (SNR), namely: 2, 3, 5, 10 and 20, whereSNR = k )r(x)} 2and r(x) denotes the average of the r(x k )s.Our main goal here is to compare the performance of the following estimators of the regression operator r: a) the Nadaraya-Watson (NW) estimator obtained for j = 0 and p = 0; b) the linear estimator obtained for j = 0 and p = 1; c) the quadratic estimator obtained for j = 0 and p = 2.

Fig. 1 .

 1 Fig. 1. A sample of 100 curves of parabolic and sinusoidal curves.

Fig. 2 .

 2 Fig. 2. Parabolic curves: y k s versus r(x k )s with SNR = 10.

Fig. 3 .

 3 Fig. 3. Parabolic curves: SNR = 10, OU error process: y k s versus their estimations.

Fig. 4 .

 4 Fig. 4. Parabolic curves: SNR = 10, ARFIMA error process: y k s versus their estimations.

Fig. 5 .

 5 Fig. 5. Sinusoidal curves: y k s versus r(x k )s with SNR = 10.

Fig. 6 .

 6 Fig. 6. Sinusoidal curves: SNR = 10, OU error process: y k s versus their estimations.

Fig. 7 .

 7 Fig. 7. Sinusoidal curves: SNR = 10, ARFIMA error process: y k s versus their estimations.

Fig. 8 .

 8 Fig. 8. The spectrometric curves.

Fig. 9 .

 9 Fig. 9. Spectrometric data: Estimated versus observed fat values for the three estimators with L 2 metric (top panels) and second derivative semimetric (bottom panels).

12 δ=1

 12 {ŷ 53 (δ)y 53 } 2 over the 12 estimated temperatures. They are as follows: NW Linear Quadratic 0.233 0.617 2.624

  x i , x)d (x j , x)ρ (ij) x i , x)d (x j , x)ρ (ij) x i , x)d (x j , x)ρ (ij) x i , x)d (x j , x)ρ (ij).

1 0µ-K( 1 )n 2 h k+ 1 0h k+ 1 0 1 0

 11111 k, (hs, h)K (s)dsµ k, (h, hs)K (s)ds + n 2 µ k, (hs, ht)K (s)K (t)dsdt = n 2 h k+ µ k, (h, h){K 2 (1) -K(1) (s, t)K (s)dsK (t)dt}{1 + o(1)} = n 2 h k+ µ k, (h, h)C 1,1 k, (h){1 + o(1)}.Therefore, combining the two asymptotic expressions, we haves * k, = nh k+ ϕ k+ (h)C 2 k+ (h)ρ (0){1 + o(1)} + n 2 h k+ µ k, (h, h)C 1,1 k, (h){1 + o(1)}.It follows that S * n = nHΣ(2) h Hρ (0){1 + o(1)} + n 2 HΣ (1,1) h H{1 + o(1)} and var{ β(x)} = S -1 n S * n S -1 n

n 2 1 0( 1 - 2 j 2 i

 1122 u)u -γ du ϕ 0 (h) 2 {nϕ 0 (h)} γ where the notation a n b n is used to indicate that a n = O(b n ) or b n = O(a n ) as n → ∞. Under the assumptions of Corollary 1, we deduce that var{ r j (x)} = O 1 nh 2 j ϕ 0 (h) {nϕ 0 (h)} γ = O 1 h 2 j {nϕ 0 (h)} γ since nϕ 0 (h) → ∞ as n → ∞ and γ ∈ (0, 1).Proof of Corollary 4. From assumptions (3), (4) and (5) and using the same notations as in the previous proof, we have, for some constant C,µ 0,0 (h, h) = 1 n 2 1≤i, j≤n i j |ρ (ij)| 1 B(x,h) (x i )1 B(x,h) (x j ) ≤ 2 n , j∈Dx (h) i≤ j |ρ (ij)| ≤ 2 n 2 n x (h) k=1 {n x (h) -k} |ρ (k)| ≤ n x (h) n 2 n k=1 |ρ (k)| .

Table 1

 1 Parabolic curves: RMSE of each estimator in terms of error processes, ARFIMA and OU, for different values of SNR.

	Error Process SNR	NW	Linear	Quadratic
		20	0.00162 0.00114	0.00063
		10	0.00310 0.00230	0.00136
	OU	5	0.00604 0.00471	0.00294
		3	0.00998 0.00797	0.00523
		2	0.01488 0.01213	0.00814
		20	0.02805 0.01926	0.01310
		10	0.05581 0.03729	0.02569
	ARFIMA	5	0.10551 0.06778	0.04731
		3	0.15977 0.09971	0.07024
		2	0.21188 0.13002	0.09211

Table 2

 2 Sinusoidal curves: RMSE of each estimator in terms of error processes, ARFIMA and OU, for different values of SNR.

	Error Process SNR	NW	Linear	Quadratic
		20	0.10374 0.01430	0.01016
		10	0.11012 0.02055	0.01147
	OU	5	0.12091 0.03200	0.01763
		3	0.13366 0.04515	0.02892
		2	0.14819 0.06169	0.04375
		20	0.12086 0.04064	0.03324
		10	0.15040 0.07399	0.06397
	ARFIMA	5	0.20712 0.13691	0.11923
		3	0.27356 0.20925	0.18117
		2	0.34200 0.28051	0.24283

Table 3

 3 Spectrometric data: RMSE of each estimator with L 2 metric and semi-metric based on second derivative.

	Distance	NW	Linear	Quadratic
	L 2	0.24886 0.06371	0.07866
	Second derivative 0.01292 0.00659	0.00635

dt,we conducted simulations studies with = 2 using the same type of curves and regression operators and arrived at the same conclusions as in the metric space considered in Section 4 ( = 0). However, the theoretical results in this work cannot be generalized to a semi-metric space because the polynomial expansion of the regression operator is not available. The case with random designs and correlated errors involves more technical difficulties and should be treated separately in possible further work. For instance, we can refer to[START_REF] Benhenni | Regression models with correlated errors based on functional random design[END_REF] when the Nadaraya-Watson estimator of r is used. Extensions of this work to semi-parametric functional models could also be considered in future research.

Acknowledgments

We would like to thank three anonymous reviewers, Prof. Philippe Vieu (the Managing Editor), and Prof. Christian Genest (the Editor-in-Chief) for their constructive comments and valuable suggestions which have led to an improved version of this article.

Appendix

In the following proofs, we will use the fact that from assumptions (3), ( 4) and ( 5), for any ν ∈ {1, 2} and as h → 0,


Abderazek Talbi

Pierre-Laurent Coddet

Malek Tabbal

Anne-Lise L Thomann

Eric Millon

Arnaud Stolz

Chantal Boulmer-Leborgne

Gerard M O'connor

Nadjib Semmar
email: nadjib.semmar@univ-orleans.fr

Comparative study of laser induced periodic surface structures formed on pulsed laser deposited and magnetron sputtered titanium oxide films

Keywords: LIPSS, Nanostructuring, Pulsed laser deposition, Magnetron sputtering, Titanium oxide thin films, Femtosecond laser

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Laser irradiation of materials offers a direct one-step route to induce micro and nanostructures on a variety of solid surfaces in view of tailoring their properties for specific technological applications [START_REF] Bonse | Laser-Induced Periodic Surface Structures -A Scientific Evergreen[END_REF][START_REF] Ahn | Low damage electrical modification of 4H-SiC via ultrafast laser irradiation[END_REF][START_REF] Vorobyev | Direct femtosecond laser surface nano/microstructuring and its applications[END_REF][START_REF] Reif | On the physics of selforganized nanostructure formation upon femtosecond laser ablation[END_REF][START_REF] Gurevich | Laser Induced Periodic Surface Structures induced by surface plasmons coupled via roughness[END_REF]. The control of the characteristics (shape, size, orientation and periodicity) of such Laser Induced Periodic Surface Structures (LIPSS) has been boosted by the advent of ultra-short lasers (particularly femto-second lasers). In addition to the Low Spatial Frequency LIPSS (LSFL) whose periodicity is typically very close to the laser wavelength, femto-second lasers also permit the formation of High Spatial Frequency LIPSS (HSFL). These latter have a periodicity that is significantly lower than that of LSFL thus widening the range of different topologies that can be formed at the surface of a material through laser irradiation.

Furthermore, there has been growing interest in the investigation of LIPSS formation on materials in thin film form [START_REF] Bobzin | Laserstructured high performance MS coatings[END_REF][START_REF] Rebollar | Laser induced periodic surface structures on polymer films: From fundamentals to applications[END_REF] (as opposed to bulk material) in view of the ever-growing applications of thin films and coatings in a wide variety of industrial domains such as aerospace, medicine, electronics and optics. For thin films, LIPSS formation brings additional challenges [START_REF] Zamfirescu | The role of the substrate material type in formation of laser induced periodical surface structures on ZnO thin films[END_REF][START_REF] Hashida | Crystal structures on a copper thin film with a surface of periodic self-organized nanostructures induced by femtosecond laser pulses[END_REF] that include the effect of interfacial stresses between film and substrate, cavitation, physical properties of the underlying substrate as well as the crystalline properties of the films that often depends on its initial elaboration.

In this paper, we present a comparative study of the formation of LIPSS using a femtosecond laser on thin films titanium oxide (TiOx) synthesized by Pulsed Laser Deposition (PLD) and Magnetron Sputtering (MS). This material was selected because of its wide applications in devices for harvesting energy (photo-voltaics, photo-catalyst) as well as its good stability at high temperature and its low negative environmental impact [START_REF] Zamfirescu | The role of the substrate material type in formation of laser induced periodical surface structures on ZnO thin films[END_REF]. The effect of laser fluence and laser scanning parameters on the properties of the LIPSS were investigated using two different femtosecond laser sources and compared for both types of films (i.e. PLD and MS).

Experimental

Two deposition techniques, namely Pulsed Laser Deposition (PLD) and Magnetron Sputtering (MS) were used to grow titanium oxide thin films with typical thicknesses ranging between 300 and 500 nm on 1.2 μm thick SiO2 layers thermally grown on Si (100) substrates.

For the PLD experiments, a KrF excimer laser beam (248 nm, 10 Hz) was focused onto a stoichiometric ceramic TiO2 target placed at a distance of 50 mm from a heated substrate holder, under a 10 -6 mbar vacuum base pressure. Film growth was performed on the conventionally cleaned substrates at a laser fluence of 4 J/cm 2 whereas the substrate was set at 700•C in order to induce a better film crystallinity. Under these experimental growth conditions, the homogeneous films were synthesized over an area of 1x1 cm 2 . They were found to be substoichiometric with a ratio of O/Ti ranging from 1.6 to 1.8, as determined by Rutherford Backscattering Spectrometry (RBS) [START_REF] Boulbar | Pure and Nb-doped TiO1.5 films grown by pulsed-laser deposition for transparent p-n homojunctions[END_REF]. For the MS deposition runs, a 10 cm diameter and 4.0 mm thick titanium metallic target (Purity 99.995 %) was magnetron sputtered in a pulsed DC mode using a mixture of Ar and O2 gas. Details of the sputtering system are available elsewhere [START_REF] Coddet | Characteristics and properties of a magnetron sputtered gadolinia-doped ceria barrier layer for solid oxide electrochemical cells[END_REF]. Film growth was thus achieved on unheated rotating substrate holder placed at a distance of ~127 mm from the target at a total pressure of 1 Pa (base pressure of ~10 -5 Pa). The flow rates of argon and oxygen were fixed at 20 and 1.5 sccm, respectively. Under these conditions, homogenous (over an area of 5x5 cm 2 ) titanium oxide thin films with a Ti/O ratio of 1.8, is also determined by RBS.

The first set of experiments to micro/nanostructure the PLD and MS grown TiOx thin films was performed using a 100 fs linearly polarized Ti: Sapphire laser emitting at 800 nm and operating at a repetition rate of 1 kHz. The laser spot diameter was focused to 350 µm on the surface of the films. The third harmonic at 266 nm, generated by tripling the frequency of the laser beam utilizing a BBO crystal, was used for these experiments, with a pulse duration of 100 fs. Treatment of large surfaces was accomplished by fixing the samples and a mirror reflecting the laser beam onto the sample was controlled by an x-y-translation device taking into account the laser pulse number by adjusting the translation step and working at a scan speed of 50 mm/s.

More detailed information about these experimental parameters can be found in Talbi et al. [START_REF] Talbi | Nanostructuring of titanium oxide thin film by UV femtosecond laser beam: From one spot to large surfaces[END_REF].

Additional laser irradiation experiments were also performed on the MS films only using a Yb:YKW 500 fs linearly polarized laser emitting at 1030 nm that can operate at a repetition rate up to 300 kHz. An irradiation configuration similar to the one using for the Ti-sapphire laser was used, but at a scan speed of 4 mm/s and a repetition rate of 100 kHz. The beam waist in this case was determined to be 32 µm, as estimated using the Liu method [START_REF] Liu | Simple technique for measurements of pulsed Gaussian-beam spot sizes[END_REF].

Results and discussion

Figure 1 represents SEM images of TiO1.8 thin films deposited on SiO2 using (a) PLD and (b) MS. The surface morphology of the PLD film appears very smooth with a barely distinguishable nano-grained structure. This contrasts with a very obvious granular structure of the MS films and a rougher surface than their PLD counterparts. Cross-sectional SEM images (not shown) indicate a columnar growth in both cases, so the difference in surface morphology could be due to surface mobility effects of the incident species during film growth. Indeed, the kinetic energy species in the laser ablation plume are typically greater than that of sputtering species (in the absence of electrical bias of the substrate, as is the case here). Furthermore, the PLD films were grown at a substrate temperature of 700 o C whereas the MS films were not heated during growth (a substrate temperature increase of no more than 50 o C would be expected though). This difference in growth temperature could also explain the difference in surface roughness between the PLD and the MS films.

=

Where d is the spot diameter and s is the step of translation along the x-axis or y-axis [START_REF] Talbi | Nanostructuring of titanium oxide thin film by UV femtosecond laser beam: From one spot to large surfaces[END_REF]. The SEM images for both samples, whether deposited by PLD or MS, show that laser irradiation under these conditions lead to the formation of regularly spaced dots with two periods at 266 and 160 nm, the latter being nearly half of the laser wavelength. Typical heights of these dots were estimated at 40 nm by preliminary Atomic Force Microscopy (AFM) measurements. These periods could be considered as LSFL and HSFL respectively, but the patterns are wavy and do not strictly follow a direction that is either parallel or perpendicular to the beam polarization. This effect is probably due to the overlap of the successive beams. These patterns are similar to those obtained under static beam conditions (one laser spot without scanning) and were mainly attributed to a self-organization process. In such a process, the melting of the surface layers can be induced by the ultra-short laser pulse at UV wavelength, where material absorption is high and depth penetration is low. Melted material would then agglomerate in regular dots as a result of surface energy minimization [START_REF] Huynh | Evidence of liquid phase during laser-induced periodic surface structures formation induced by accumulative ultraviolet picosecond laser beam[END_REF].

In a previous paper [START_REF] Talbi | Nanostructuring of titanium oxide thin film by UV femtosecond laser beam: From one spot to large surfaces[END_REF], it has been shown varying laser fluence and number of pulses can lead to a variety of different surface morphologies ranging from dots to LSFL, HSFL as well as intermediate morphologies. Thus, PLD and MS films were irradiated under different conditions, as seen in Figure 3, namely at a fluence of 25 mJ/cm 2 and a total number of pulses of 3300. The obtained morphology consists of LSFL with a period of ~ 240 nm and an orientation predominantly parallel to the beam polarization. Careful observation of the LSFL shows that there are some locations where a "linking" between the LSFL occurs. The presence of such links or bridges supports the self-organization model [START_REF] Reif | On the physics of selforganized nanostructure formation upon femtosecond laser ablation[END_REF] where there is evolution of the surface from regularly spaced dots to LSFL with increasing laser fluence. Comparing the images of the PLD films (Fig. 3-a) and those of the MS films (Fig. 3-b), it is also clear that the formed patterns are quite similar quantitatively (very close spatial period of the LSFL-and HSFL-LIPSS as determined above) and qualitatively (shape/type of the LIPSS, dots or lines, depending on laser irradiation conditions). Further increase in laser fluence, up to 30 mJ/cm 2 for a number of pulses of 13000 leads to a destruction of the previously described nanostructures through a dominant ablation-evaporation process, again for both PLD and MS films. However, despite these very similar results obtained within the experimental conditions used in this work, the delamination of the MS thin films is noticed in some zones of the treated surface for a laser dose of 25 mJ/cm² and 3300 pulses (as seen in Figure 4). Such delamination effects were not observed for the PLD films. The delamination of a thin film under laser irradiation can be a result of ablation phenomena induced through a thermo-mechanical process [START_REF] Lee | Micro-scale patterning of indium tin oxide film by spatially modulated pulsed Nd:YAG laser beam[END_REF][START_REF] Farid | Stress assisted selective ablation of ITO thin film by picosecond laser[END_REF]. It is clear from Fig. 4 that, after film delamination, the substrate surface appears very clean without any damage and the delaminated film exhibits the regular LSFL pattern, inferring that the delamination occurs after or simultaneously to the LSFL formation. It is also noteworthy through observation of the delaminated layers in Fig. 4 that the amplitude of the LIPSS pattern is much less than the thickness of the film. We can thus conclude that de-wetting of the substrate essentially plays no role in the LIPSS formation. Since, no such delamination was observed in the case of PLD thin films, it can be concluded that the adhesion and thermo-mechanical properties of the latter films on the substrate are better than their MS counterparts. This could be attributed to the nature of the PLD process whereby the highly energetic particles present in the ablation plume and impinging on the substrate and on the film as it grows, can enhance surface adhesion on the substrate. Furthermore, high internal compressive stress is expected to be present in the MS films because they were grown at room temperature, which is not the case of the PLD films. Yet, the delamination of the MS films occurs only after laser irradiation. It is widely accepted in literature that, the cracking, fragmentation and delamination of material occurs due to the thermo-elastic stresses induced by laser irradiation [START_REF] Lee | Micro-scale patterning of indium tin oxide film by spatially modulated pulsed Nd:YAG laser beam[END_REF][START_REF] Farid | Stress assisted selective ablation of ITO thin film by picosecond laser[END_REF][START_REF] Connor | Laser patterning of very thin indium tin oxide thin films on PET substrates[END_REF]. As result of the femtosecond pulse, the laser energy absorption causes a confined heat in the irradiated volume. This thermal confinement leads to an increasing temperature within the irradiated zone leading to thermoelastic stress [START_REF] Farid | Stress assisted selective ablation of ITO thin film by picosecond laser[END_REF]. Due to the repeated rapid heating/cooling processes, a very high density of such thermo-elastic stress can be generated which causes the fracture by fatigue and delamination of the layers.

In order to further investigate the effect of laser fluence on the delamination of MS films and with the aim of nano-structuring large areas of the films (2.5x2.5 cm 2 ), MS films were irradiated using the 500 fs laser emitting at 1030 nm expecting a lower peak power. The repetition rate, the scan speed and the spacing between lines were fixed to respectively 100 kHz, 4 mm/s and 5 μm leading to a number of pulses of 24000. This much higher frequency of the laser (100 kHz) compared to the previous one (1 kHz) allowed irradiation of large area within a reasonable period of time (tens of minutes). Figure 5 shows the effect of surface morphology of MS films following irradiation at different values of laser fluence, ranging from 55 to 220 mJ/cm 2 . These fluences exceed the surface damage threshold which was estimated at 36 mJ/cm² after 10000 pulses under static laser irradiation i.e. with no beam scanning. The damage threshold was deduced by SEM observation of the damage zone expansion as the number of pulses is increased under static laser irradiation. A somewhat unexpected behavior of surface morphology evolution is observed, with increasing laser fluence. Indeed, for fluences ranging from 55 to 85 mJ/cm², film delamination is severe but the effect appears to be drastically reduced by for a fluence of 110 mJ/cm² and higher, and a micro-cracked surface morphology is observed.

A schematic showing the different damage threshold fluences for a Gaussian laser beam profile (as is the case for the 500 fs laser) is illustrated in Fig 6 . This schematic can be used to explain the effects seen in Fig. 5, through the following mechanisms. The surface of the film can be damaged through two main processes: (1) Cracks and delamination processes induced by thermo-elastic stresses (2) Ablation, evaporation and melting processes induced by the lattice heating due to heat accumulation. Assuming that, after multiple shots, the cracking (stress) threshold is lower than thermal threshold (melting), F the applied fluence is greater than the cracking (stress) threshold but lower than the ablation (evaporation and melting) threshold, laser irradiation can lead to significant film damage manifested through cracking and delamination processes in form of flakes (in our case, this would correspond to laser fluences lower than 110 mJ/cm²). On the other hand, if the applied fluence is above the melting and cracking threshold fluence but does not exceed the ablation (evaporation) threshold, the laser irradiation would induce thermo-elastic stresses leading to the generation of cracks, while simultaneously triggering a melting process. This latter can lead to stress relaxation and induce "repair' of the surface damage induced by cracks and delamination. This proposed mechanism as well as LIPSS formation under femto-second laser irradiation of MS films at a wavelength of 1030 nm will be the subject of future work.

Conclusion

In this paper, we present a comparative study of LIPSS formation on sub-stoichiometric titanium oxide thin films (TiOx, with 1.6 ≤ x≤1.8) that were grown using two different techniques, namely pulsed laser deposition (PLD) and magnetron sputtering (MS). Irradiation was performed using femtosecond lasers beams, emitting at wavelengths 266 and 1030 nm, and at different number of laser pulses and laser fluences. Conditions for the formation of LIPSS were determined for both types of films, at 266 nm. A variety of nanostructures such as LSFL, HSFL, spikes, dots were detected and their shapes and sizes were found to strongly depend on laser fluence and number of pulses, as essentially dictated by free-surface energy minimization.

LIPSS appeared to follow similar formation pattern for both PLD and MS films, except for some delamination detected for the MS films. This was attributed to better adhesion and lower internal stress of the PLD films compared to their MS counterparts. Large surface nanostructuring experiments performed at 1030 nm showed less delamination with the increase of the laser fluence, a phenomenon attributed to the competition between the cracking and the melting processes of the layers.

Figures captions:

Figure 2

 2 Figure 2 show SEM images of the (a) PLD and (b) MS TiOx films after irradiation with

1 -

 1 Figure 1-SEM images of TiOx thin films deposited by (a) PLD and (b) MS.

Figure 2 -

 2 Figure 2-SEM images of (a) PLD and (b) MS TiOx films irradiated by 100 fs laser pulses at a wavelength of 266 nm with F = 15 mJ/cm² and N = 13000.

Figure 3 -

 3 Figure 3-SEM images of (a) PLD and (b) MS TiOx films irradiated by 100 fs laser pulses at a wavelength of 266 nm with F = 25 mJ/cm² and N = 33000.

Figure 4 .

 4 Figure 4. SEM image of a delamination observed on a MS by film irradiated 100 fs laser pulses at a wavelength of 266 nm.

Figure 5 .

 5 Figure 5. SEM images of TiOx MS films irradiated at 1030 nm at different laser fluences for N= 24000.

Figure 6 .

 6 Figure 6. Schematic showing the different damage threshold fluences for a Gaussian laser beam profile.

  Figure 1

Figure 2 (

 2 Figure 2

  Figure 3

Figure 4 .

 4 Figure 4.

  Figure 5.

Figure 6 .

 6 Figure 6.

Acknowledgements

This project was supported through a French National project `'Projet Investissement d'Avenir -Tours2015' in collaboration with STMicroelectronics in Tours. MT would like to thank the University of Orléans and the SAFAR program for a financial support during the completion of this work.


est destinée au dépôt

  Interaction entre la mémoire de travail et l'attention sélective dans le domaine auditif non-verbalSalomé BlainTo cite this version:Salomé Blain. Interaction entre la mémoire de travail et l'attention sélective dans le domaine auditif non-verbal. Neuroscience. Université de Lyon, 2020. English. NNT : 2020LYSE1179. tel-03392189