Rafael Barreto Jijón

Carlos Canudas-De-Wit

Silviu-Iulian Niculescu

Jonathan Dumon

Adaptive Observer Design under Low Data Rate Transmission with Applications to Oil Well Drill-string

Keywords: Stick-Slip, Oil Well drill string, D-OSKIL, unknown parameter adaptive observer, time-variant, delay, stability

In oil well drilling operations, one of the important problem to deal with is represented by the necessity of suppressing harmful stick-slip oscillations. A control law named D-OSKIL mechanism uses the weight-on-the-bit force as a control variable to extinguish limit cycles. It uses the value of the bit angular velocity that is found through an unknown parameter observer by means of the measure of the table rotary angular speed. To improve this former estimation, we add the measurement of the angular velocity of the bit that, due to the technological constraints, arrives delayed. This new design leads us to the analysis of a time-varying delay system.

I. INTRODUCTION

Oil well drilling operations present a particular friction phenomenon called Stick-Slip Oscillations which consists in a sub-normal irregular rotation movement, e.g., the top of the drill string rotates with a constant angular velocity, meanwhile downhole, the bit (cutting device) angular velocity varies in a range from zero to six times the one measured at the surface. This effect appears mainly when the bit is in contact with rock formations [START_REF] Armstrong-Hélouvry | Stick-slip arising from stribeck friction[END_REF], [START_REF] Kyllingstad | A study of stick/slip motion of the bit[END_REF].

Due to the very large forces that are applied into the bit, the presence of stick-slip self-excited oscillations can cause irreversible damage to the equipment, as for example, shut-downs, operation delays, decreasing service life of drill strings and downhole equipment. Since avoiding this kind of oscillations can provide important savings in terms of exploitation time, spare parts costs and maintenance this task has become a challenge for drillers and scientists. Further details about stick-slip and drilling oscillations can be found in [START_REF] Zamanian | Stick-slip oscillations of drag bits by considering damping of drilling mud and active damping system[END_REF], [START_REF] Baumgart | Stick-slip and bit-bounce of deep-hole drillstrings[END_REF], [START_REF] Navarro-Lopez | Vibraciones mecánicas en una sarta de perforación: Problemas de control[END_REF] and [START_REF] Spanos | Oil well drilling: A vibrations perspective[END_REF].

A model for the drill string system and for the stickslip oscillations as well as an appropriate control law called Drilling-Oscillation Killer (D-OSKIL) have been proposed in [START_REF] Canudas-De Wit | D-OSKIL: A new mechanism for controlling Stick-Slip Oscilations in oil well drillstrings[END_REF]. The controller uses mainly a vertical force sometimes named "weight on the bit" (WoB) as an additional variable to eliminate stick-slip effects. Additionally, not all the values of the states of the system are available for control. Adaptive observer has been designed to provide an estimation of such R. Barreto Jijón is with the INPG, GIPSA-Lab, NeCS team. Grenoble, France. (e-mail: rafael.barreto.jijon@gmail.com) carlos.canudas-de-wit@gipsa-lab.inpg.fr)

S.-I. Niculescu is with the L2S, CNRS-Supelec, Gif-sur-Yvette, France.

(e-mail: Silviu.Niculescu@lss.supelec.fr)

J. Dumon is with the CNRS, GIPSA-Lab, NeCS team. Grenoble, France.

(e-mail: jonathan.dumon@gipsa-lab.inpg.fr)

states and, additionally, to estimate the friction coefficient. In this case, the measured variable is the rotary angular speed. Since the physics effects occurring downhole have no strong influence at surface due to the attenuation along the drill string, the measurement does not effectively reflect them. This means that the signal to noise ratio is small and it affects the quality of estimation. We are to improve the observer's behavior, using some coarse information coming from the drilling toll with a new measurement: the bit angular velocity. As explained in the forthcoming section, the technological constraints induce the presence of a transmission delay. It is worth to mention that such a delay introduces an additional difficulty in the control problem as pointed out in the literature in different other cases [START_REF] Gu | Stability of time-delay systems[END_REF], [START_REF] Niculescu | Delay effects on stability. A robust control approach[END_REF], [START_REF] Michiels | Stability and stabilization of timedelay systems. An eigenvalue based approach[END_REF] The paper is organized as follows: In section II, we briefly introduce the model of the dynamical system along with a description of the testbed where the observer is to be implemented as well as some operational conditions. Section III includes some discussions on the particular physical constraints as well as the technological feasibility of the implementation of the observer. In section IV we tackle the observer's development and its stability conditions and finally, some concluding remarks in section VI end the paper.

Throughout the paper, the following notations will be used: the space C v n,τ is a set defined by {φ ∈ C n,τ : φ c < v} where C n,τ = C ([-τ, 0], R n ) represents the Banach space of the definite piecewise-continuous vector functions mapping the interval [-τ, 0] into R n with the uniform convergence topology, and

φ c = sup -τ≤t≤0 φ (t)
represents the norm of the function φ ∈ C n,τ .

II. SYSTEM DESCRIPTION

A. Model

The system is modeled as two coupled masses as shown in Figure 1. J r and J b are two inertial masses locally damped by d r and d b . The inertias are coupled through an elastic shaft of stiffness k and damping c. Let us define ϕ r , ϕ b as the angular positions of the rotary and the bit respectively; φr , φb as their angular velocities, u(t) = WoB is the weight on the bit control signal, v(t) is the rotary table torque control signal used to regulate φr , µ is the friction coefficient; A, B, H, C o are model matrices given in (3), Ψ(t) = Ψ(u(t)) = Hu(t), x is the state vector and y o is the output variable. For a more detailed description, see, for instance, [START_REF] Canudas-De Wit | D-OSKIL: A new mechanism for controlling Stick-Slip Oscilations in oil well drillstrings[END_REF]. The model to be used is described by: ẋ(t) = Ax(t) + Bv(t) + Ψ(t)µ

(1)

y o = C o x = φr ,
where the state x = [x 1 x 2 x 3 ] T is defined as follows:

x 1 = ϕ r -ϕ b , x 2 = φr , x 3 = φb ( 2
)
and

A =    0 1 -1 -k J r -d r +c J r c J r k J b c J b -c+d b J b    , B =   0 1 J r 0   H =   0 0 -1 J b   , C o = 0 1 0 (3)

B. Testbed

In order to validate the proposed results, an appropriate testbed has been built. It is shown in Figure 2 as well as its schematic in Figure 3. It has a Host PC for development, compiling and user interface. A target PC is also used for real-time execution of the code like the observer and the D-OSKIL controller [START_REF] Lu | Experimental study of the D-OSKIL mechanism for controlling the stick-slip oscillations in a drilling labratory testbed[END_REF]. The rotary system is composed of a DC motor, a transmission box, the rotary table, the drill string, a bit and a quadrature encoder. The support platform consists of a DC motor to move the rotary system vertically, a specimen holder, a support structure and a tension force sensor. The acoustic signal is simulated as sonar pulses like a beeper and microphone pear as shown in Figure 3. Furthermore, the delay is added artificially by software.

We will define the following operation conditions: Well depth: d ∈ [0, 8000] meters, since drilling penetration speed is usually very small, we can consider it as negligible: ḋ ≈ 0. In practice, we have found that φb ∈ [0, 31] rad/s with typical value φb = 5 rad/s.

III. WIRELESS-TRANSMISSION TECHNOLOGIES

In this section, we focus on the way to bring the measurement from downhole to surface so we can use it for improving the observer's behavior.

There are mainly two types of transmission: through telemetry signals along the drilling fluid often referred to as mud-pressure pulses [START_REF] Jeffryes | Method and apparatus for enhanced acoustic mud-pulse telemetry[END_REF] and through acoustic waves along the drillstring [START_REF] Clayer | Procède de transmission acoustique de données de forage d'un puits[END_REF].

In most of the literature, electronic equipments are designed for data acquisition and to play the modulator roll. It should be implemented as an autonomous system energized either by a mud operated electrical turbine or by a battery pack [START_REF] Tubei | Mud-pulser telemetry system for down hole Measurement-While-Drilling[END_REF].

A. Mud-pulse telemetry

This technology uses the mud that goes through the drilling system as a transmission media. The data will be represented by pressure pulses. According to [START_REF] Tubei | Mud-pulser telemetry system for down hole Measurement-While-Drilling[END_REF], the pulser actuator (a stepper-motor-based device) and a main valve restricts the flow and creates some pressure-pulse sequence. A piezoelectric device captures these variations that are then analyzed by a micro-controller. Evidently, due to the irregular nature the mud flow, the low frequency vibrations produced by mud pumps and pulsation dampeners the signals are corrupted by noise. Furthermore, they have an important attenuation. Some characteristics to highlight are [START_REF] Liu | Transmission behavior of mud-pressure pulse along well bore[END_REF], [START_REF] Clayer | Procède de transmission acoustique de données de forage d'un puits[END_REF], [START_REF] Drumheller | The propagation of sound waves in drill strings[END_REF] its cost-effective data transfer, its very low bite rate (1 or 2 bits per second). Mud-pulse velocity declines with the disturbances of mud density, gas content and mud compressibility. It becomes more difficult with increasing well depth. Pulse waves travel through the borehole at 1200 meters per second [START_REF] Liu | Transmission behavior of mud-pressure pulse along well bore[END_REF], hence the measure arrives with some delay that increases up to τ max ≈ 6.6 seconds .

B. Acoustic data transmission over a drill string

Since the acoustic wave propagation velocity in the string material is at least three times superior to that in the mud of the borehole [START_REF] Clayer | Procède de transmission acoustique de données de forage d'un puits[END_REF], and a higher transmission rate is possible (typically 6 bps), acoustic transmission seems to be the best way to emit pulses to the surface. These acoustic waves are generated torsional contractions generated by magneto restrictive rings set inside the pipe [START_REF] Drumheller | An overview of acoustic telemetry[END_REF]. In this case τ max ≈ 2.2. It is useful to note that there exists an attenuation of around 4 dB/300m [START_REF] Drumheller | Acoustical properties of drill strings[END_REF]. However, we can neglect it because there is always a possibility of setting a repeater at any joint at each 10-15 meters of the section. We consider that this fact does not add any extra considerable delay since the repeater's amplification can occur almost instantaneously.

The telemetry system sends signals directly to the surface through the channel. Usually, there is an embedded sensor measuring φb downhole. A measurement noise S(t) is added to the data and then coded all together in such a way it can be transmitted through the acoustic channel G. At surface, a receiver will read the encoded signal with the noise N(t). Furthermore, a digital algorithm is used to decode this data and make it available for the use of the observer. Both methods can be modeled by the schematic shown in Figure 4.

C. Transmission delay range and friction hypothesis

Due to technical considerations we can assume that the transmission media is, as a first approximation, like a pure delay system with delay time τ ∈ [0, τ max ]. Moreover, the well's depth increases at a very slow rate and it stops each 10-15 meters. In this procedure, the delay can be recalculated. Hence, the delay can be defined as a constant, that is τ = 0.

On the other hand, we will consider that the friction coefficient is constant or at least slow time variant μ ≈ 0. This approximation is often assumed in the context of adaptive control. This hypothesis means that the rate of variation of the rock friction coefficient does not exhibit substantial changes during drill-operation. Even if the drilled surfaces may have different friction characteristics, the rate of penetration remains small ( ḋ ≈ 0).

IV. ADAPTIVE OBSERVER DESIGN WITH DELAYED FEEDBACK SIGNALS

A general architecture is proposed in Figure 5 where Σ is the system model, Σ is the observer, x is the observed state vector, K o is the default observer gain, K m is the observer gain when the value of φb becomes available and G is the transmission channel. Here the outputs of the system are y 1 = φr , and y 2 = φb .

A. Main Results

In this section, we will focus on a particular extension of the original observer designed in [START_REF] Canudas-De Wit | D-OSKIL: A new mechanism for controlling Stick-Slip Oscilations in oil well drillstrings[END_REF] to handle the delay presence and noise requirements. In this sense, we add two gains related to the delayed measurement as follows: ẋ(t) = A x(t) + Bv(t) + Ψ(t) μ(t)

+[K + β Γ(t)Γ T (t)C T o ][y o (t) -C o x] +K x [y m (t -τ) -C m x(t -τ)] (4) μ(t) = β Γ T (t)C T o [y o (t) -C o x] +K µ [y m (t -τ) -C m x(t -τ)] (5) Γ(t) = (A -KC o )Γ(t) + Ψ(t) (6)
where, y m (tτ) is the new delayed output vector, C m ∈ R 1×3 , K x ∈ R 3 and K µ ∈ R are the output matrix, and gain matrices for the delayed measurement. Γ(t) ∈ R 3 is a matrix generated by the system free of delay ( 6) such that it generates the states dependent of µ. Here, β is a positive scalar and v(t) (and then Ψ(t)) is persistently exciting, for example ∃δ , T > 0 such that the following inequality holds:

t+T t Γ T (ξ )C T o C o Γ(ξ )dξ > δ > 0 (7)
Theorem 1. The unknown parameter adaptive observer described by ( 4)-( 6) guarantees exponential stability for small delays τ ∈ [0, τ) if there exists some gain matrices K x and K µ such that

K x (t) = Γ(t)K µ (8)
and the following conditions hold simultaneously:

(i) (A -KC o )
is Hurwitz, and (ii) the inequalities

0 < a min -K µ C m Γ * (9) τ < τ = 2 a min -K µ C m Γ * | K µ C m Γ * | (2 | K µ C m Γ * | +1 + a 2 max )
, [START_REF] Kyllingstad | A study of stick/slip motion of the bit[END_REF] with:

a min = min t≥0 (β Γ T (t)C T o C o Γ(t)), a max = max t≥0 (β Γ T (t)C T o C o Γ(t)).
Proof 1. Following the method described in [START_REF] Zhang | Adaptive observer for MIMO linear time varying systems[END_REF] with the constraint (11), we introduce ( 5) into ( 4) and then we have the expression of the system [START_REF] Lu | Experimental study of the D-OSKIL mechanism for controlling the stick-slip oscillations in a drilling labratory testbed[END_REF].

K x (t) = Γ(t)K µ (11) ẋ(t) = A x(t) + Bv(t) + Ψ(t) μ(t) +K(y o (t) -C o x(t)) + Γ(t) μ(t) ( 12
)
Let x(t) = x(t)x(t), μ(t) = μ(t)µ with μ = 0 and since μ(t) = μ(t) -μ we obtain:

ẋ(t) = (A -KC o ) x(t) + Ψ(t) μ(t) + Γ(t) μ(t) ( 13
)
We define the following variable transformation:

η(t) = x(t) -Γ(t) μ(t) (14)
then we obtain:

η(t) =(A -KC o )(η(t) + Γ(t) μ(t)) + Ψ(t) μ(t) -Γ(t) μ(t) =(A -KC o )η(t) + [(A -KC o )Γ(t) + Ψ(t) -Γ(t)] μ(t)
Due to the system (6), we have:

η(t) = (A -KC o )η(t). (15)
The system (15) can be assumed as strictly stable e.g. the constant pair (A -KC o ) is detectable. Now, we will focus on analyzing the behavior of μ:

μ(t) = β Γ T (t)C T o (y o (t) -C o x(t)) + K µ C m x(t -τ) = -β Γ T (t)C T o C o x(t) + K µ C m x(t -τ) = -β Γ T (t)C T o C o (η(t) + Γ(t) μ(t)) + K µ C m x(t -τ) (16) Since x(t -τ) = η(t -τ) + Γ(t -τ) μ(t -τ) then μ(t) = -β Γ T (t)C T o C o (η(t) + Γ(t) μ(t)) +K µ C m [η(t -τ) + Γ(t -τ) μ(t -τ)]
Define now:

a(t) = β Γ T (t)C T o C o Γ(t) (17) b(t) = -K µ C m Γ(t -τ) ( 18
)
where 0 ≤ a(t) ≤ max t≥0 (a(t)) = a max and b(t) are scalars.

On the other hand, we can approximate

Γ * (t) = -(A -KC o ) -1 Ψ * (t) = -(A -KC o ) -1 Hu * (t)
where Γ * , Ψ * and u * denote the steady state of Γ(t), Ψ(t) and u(t) respectively. Physically, u(t) is a vertical force moving heavy rotating masses with friction. Thus, this movement can not be neither so fast nor sudden. Therefore, we can consider that u * and then Γ * vary slowly and they can be treated as constants. Since µ changes faster than the variability of the gain, we can assume b(t)

≈ -K µ C m Γ * b. Then μ(t) = -a(t) μ(t) -b μ(t -τ) + f [η(t), η(t -τ)] (19)
Summarizing, we start with the autonomous system η. We get an adaptation equation by injecting it into that of μ. With these "present" and "past" values we compute a estimation of the state vector. The development of the observer's equations are resumed in Figure 6: Since f [η(t), η(tτ)] → 0 while t → ∞ the problem becomes, at the next step, by analyzing the homogenous part of the equation, to find conditions such that system μ(t) = -a(t) μ(t)b μ(tτ) is exponentially stable, and then explicitly compute the gains. For sake of writing simplicity let us define χ(t) = μ(t). Then the system becomes:

χ(t) = -a(t)χ(t) -bχ(t -τ) (20)
We have the following result (see the appendix for the proof): Lemma 2. The time variant delay system described by ( 20) is asymptotically stable for all τ ∈ [0, τ) if the following conditions hold:

a min + b > 0 (21) τ < 2 a min + b 2b 2 + (a 2 max + 1) | b | (22)
Hence, if there exists ε > 0, τ > 0 and a max > 0, a min ≥ 0 verifying the conditions in Lemma 2, it is possible to design a gain b > 0 that makes system (20) asymptotically stable.

Let σ , ε, λ three positive scalars. We know that

0 < σ χ(t) 2 ≤ V (χ t ) V (χ t ) ≤ -ε χ(t) 2 ≤ 0,
where V represents the Lyapunov-Krasovskii functional allowing us to conclude on the exponential stability of (20):

V (χ t ) V (χ t ) ≤ - ε σ
integrating in a positive sense we get:

t 0 V (χ ρ )) V (χ ρ )) dρ ≤ t 0 - ε σ dρ (23)
Finally, we obtain:

V (χ t ) ≤ λ e -ε σ t (24)
This means that V (χ t ) decreases exponentially, that implies that the system converges to zero exponentially fast. Thus, if the conditions shown above are held, the system is globally exponentially stable. This concludes proofs of theorems 1 and 2.

V. EXPERIMENTS This section reports some preliminary experimental results carried out in the test-bed shown in Fig. 2 and Fig. 3. The details of this test bed are described in Section II-B. Further details of this experimental test-bed can also be fund in [START_REF] Lu | Experimental study of the D-OSKIL mechanism for controlling the stick-slip oscillations in a drilling labratory testbed[END_REF].

The experiments correspond to a total run of 77 sec under drilling conditions with constant WoB. These data are collected during the whole run to compute the mean and variance of the estimation error in the two compared cases: a) Observer without delayed measurements.This corresponds to the case where the gains are set to zero K x = K µ = 0 b) Observed using delayed measurements as proposed in this paper with K x , K µ = 0. Table I shown the statistic improvements of our proposed observer over the one that do not use the information.

VI. CONCLUSIONS

We have discussed the design of an observer for well drilling. We analyzed the system model, the observer's technological feasibility with respect to communication constraints, and operational conditions. Then we designed the general architecture that the new observer should have. We designed an unknown parameter adaptive observer, adding a gain related to a delayed measurement. To prove its exponential stability, we found some simple conditions between the gains of the known and unknown state equations that lead us to a problem of the stability of a linear time-variant Fig. 7. Time-evolution of the measured bit rotational velocity and the estimation using (lower) and not using (upper) the delayed measurements.

system with delay. We focused in the analysis of such a problem and we found some specific relationships between the former observer's gain and that of the delayed term with respect to the delay time through the Lyapunov's method.

VII. ACKNOWLEDGEMENTS This work has been supported by the FeedNetBack European project FP7-ICT-2007-2.

APPENDIX

Let us rewrite [START_REF] Zhang | Adaptive observer for MIMO linear time varying systems[END_REF] as

χ(t) = -a(t)χ(t) -bχ(t) + bχ(t) -bχ(t -τ) = -(a(t) + b)χ(t) + b 0 -τ χ(t + θ )dθ = -(a(t) + b)χ(t) -b 0 -τ a(t + θ )χ(t + θ )dθ -b 2 0 -τ χ(t + θ -τ)dθ (25)
Since according to [START_REF] Niculescu | Sytèmes à retard[END_REF], (25) can be viewed as an equation in the C v 1,2τ space we have: Let us define a Lyapunov-Krasovskii functional candidate V (χ t ) of the form

V (χ t ) = V 1 (χ(t)) +V 2 (χ t ) +V 3 (χ t ), (27) with:
V 1 (χ(t)) = 1 2 χ(t) 2 ≥ 0 (28) Finally, adding all Lyapunov-candidates and using the property that 0 ≤ a min ≤ a(t) ≤ a max , it follows that V (χ(t)) = V1 (χ(t)) + V2 (χ(t)) + V3 (χ(t))

V 2 (χ t ) = 1
≤ -(a min + b -b 2 τ - | b | a 2 max 2 τ - | b | 2 τ)χ(t) 2 ≤ -ε χ(t) 2 (33)
where ε > 0. Hence, according to the Lyapunov-Krasovskii theorem, we need to find a condition for (33) to be negative definite so we can have asymptotic stability for the system (20):

0 < a min + b (a 2 max + 1 + 2 | b |) | b | τ < 2(a min + b)
This concludes the proof.

Fig. 1 .

 1 Fig. 1. Drill string two-coupled masses model.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Experimental setup

Fig. 4 .

 4 Fig. 4. Block diagram of the transmission channel.

Fig. 5 .

 5 Fig. 5. General Architecture

Fig. 6 .

 6 Fig. 6. Observer's equation development Block Diagram

Figure 7

 7 Figure 7 show the time-history of the bit rotational velocity estimations. The Figures only show a zoom of a period of 1.5 sec. The upper figure compares the true measured value of this velocity with the estimated without using the delayed measurements. The lower figure show the improvements that the proposed method.

  χ(t) = -(a(t) + b)χ(t)b 0 -τ a(t + θ )χ(t + θ )dθ -b 2 -τ -2τ χ(t + θ )dθ(26)

) 2 + θ ) 2 dθ

 22 dξ dθ ≥ 0 (30)We take the derivative of V 1 (χ(t)) along the system (20):V1 (χ(t)) = -(a(t) + b)χ(t) )a(t + θ )χ(t + θ )dθ -)χ(t + θ )dθ Since we know that: + θ ) 2 dθ(31)we can rewrite V1 as:V1 (χ(t)) ≤ -(a(t) + b)χ(t) )a(t + θ )χ(t + θ )dθ +Let us find the derivative of V 2 as:V2 (χ t ) = 1 2 b 2 τ χ(t) 2 -+ θ ) 2 dθ

TABLE I STATISTICS

 I OF THE OBSERVATIONS ERROR OVER A WHOLE RUN OF 77

SEC.


J M Azaïs

S Bercu

J C Fort

A Lagnoux

P Lé

Simultaneous confidence bands in curve prediction applied to load curves

Keywords: Gaussian vectors, maximum of Gaussian sequences, simultaneous confidence bands, thresholding, load curve

In numerous contexts, one has to forecast a full curve from some explanatory variables. For that purpose, one aims at deriving simultaneous confidence bands. In this article theoretical and numerical results on the maximum of Gaussian sequences are used to construct those simultaneous confidence bands.

Introduction

In curve prediction, one is generally interested in deriving simultaneous confidence bands: regions in which the entire curve lies with some given probability. Simultaneous confidence bands have been widely studied in the literature and the topic has known an incredible expansion since the first articles published in the early fifties. The first studies have been conducted in the context of linear models by Tukey [START_REF] Tukey | The problem of multiple comparisons[END_REF] and [START_REF] Scheffé | The analysis of variance[END_REF]. In the last paper, the construction of the simultaneous confidence bands relies on the Fisher distribution and the confidence region is an ellipsoïd. The principle of Scheffé's method has been also generalized to non linear regression and usually relies on Bootstrap and asymptotic results (see e.g. the articles of Claeskens et al. [START_REF] Claeskens | Bootstrap confidence interval for regression curve and their derivatives[END_REF], Hall [START_REF] Hall | On bootstrap confidence intervals in nonparametric regression[END_REF] and the references therein). These techniques have been also extended and adapted for various problems (see e.g. estimation of distribution [START_REF] Bickel | On some global measures of the deviations of density function estimates[END_REF][START_REF] Hall | Effect of bias estimation on coverage accuracy of bootstrap confidence intervals for a probability density[END_REF], probability functions [START_REF] Csörgo | Confidence bands from censored samples[END_REF], elements of the spectral density matrix of autoregressive processes [START_REF] Newton | Simultaneous confidence bands for autoregressive spectra[END_REF], approximation of an integral by Monte Carlo method [START_REF] Kendall | Brownian confidence bands on Monte Carlo output[END_REF]...). Here we introduce a new technique based on the maximum of Gaussian processes whose distribution is estimated by a MCQMC algorithm proposed by Alan Genz ( [8,9]) on one hand and Sidak's inequality on the other hand, which seems to be new.

We implement this technique in the numerical context of load curve prediction: power producers like EDF, the electrical French group, use the information contained in the load curves of their customers to plan electricity production and to offer them an appropriate tarification. Since this 1 information is not available for all the clients, EDF needs to estimate these load curves from a learning sample. To appreciate the accuracy of these estimations we are interested in deriving precise simultaneous confidence bands. For that purpose, we construct a model (involving three main parameters, one of which is the weekly cycle) to represent precisely the load curves from explanatory variables. After a linear regression of the curves on the explanatory variables, we use the proposed technique. Obviously, this technique may apply to any context in which one need to estimate the evolution of some observed quantity during the time by giving simultaneous confidence bands. As a consequence, the construction of such bands represents a great challenge in numerous and various applied domains and particularly in industry (e.g. signal detection by Hall et al. [START_REF] Hall | Nonparametric confidence intervals for receiver operating characteristic curves[END_REF] or Macskassy et al. [START_REF] Macskassy | ROC confidence bands : An Empirical Evaluation[END_REF]) or biology (e.g. Zouhourian et al. [START_REF] Zouhourian-Saghiri | Loi d'usure mandibulaire chez Locusta Migratoria (Orthopt. Acrididae). Son utilisation pour la datation des ailés[END_REF])... The article is organized as follows. Section 2 describes the general framework and the procedure to construct the estimates and to derive the simultaneous confidence bands. Section 3 is dedicated to the computation of the maximum of the absolute value of a Gaussian series. Section 4 presents the framework of the EDF context, while numerical results are given in Section 5. The article ends with conclusion and prospects in Section 6.

Estimation of curves in high dimension using explanatory variables

In this section, we consider a family of curves R k (x) ∈ R associated to the individual k (k = 1, . . . , K) and for a generic time x, integer in X = 1, . . . , x.

In the EDF context detailed in Section 4, it can be for example

• the entire (annual) load curve (LC k (t)) t of the customer k, k = 1, . . . , K, over time t = 1, . . . , 17520. t represents any half-hour in the year.

• the centered and standardized week load curve (S k (h)) h of this customer, over time h = 1, . . . , 336. h represents any half-hour in the week.

• the mean of every week load curve (M k (w)) w of this customer, over time w = 1, . . . , 52. w represents any week in the year.

Linear regression

We assume that each curve R k (.) belongs, up to an error term, to some linear functional space i.e.

R k (x) = j∈J a j k φ j (x) + k (x) (k = 1, . . . , K) ( 1 )
where k (x) is an extra term assumed to be a Gaussian white noise and (φ j ) j∈J is an orthonormal basis. We denote J = card(J ).

We do not observe directly the coefficients a j k but

ãj k = a j k + β j k , k = 1, . . . , K, j = 1, . . . , J, (2)
where (β j k ) is an extra error term.

Moreover we suppose that to each individual are associated P explanatory variables and that for each j a linear regression model is used:

a j k = P p=1 α j p V k,p + η j k = V (k)α j + η j k , k = 1, . . . , K, (3)
where V = (V k,p ) (k=1...K,p=1...P ) is the matrix of the regression model and V (k) the k-th row of V . Thus we have J = card(J ) linear models of the form (3) that can be rewritten as

a j = V α j + η j , for j = 1, . . . , J. (4)
Eventually,

ãj k = V (k)α j k + η j k + β j k , k = 1, . . . , K, j = 1, . . . , J. (5)
Remark 1

• In most cases, the model will contain a constant: V k,1 ≡ 1.

• We assume that the errors of the J models, defined by [START_REF] Bickel | On some global measures of the deviations of density function estimates[END_REF], are independent with different variances. Since k (x) is a Gaussian white noise, the errors of the J models, defined by [START_REF] Csörgo | Confidence bands from censored samples[END_REF], are also independent with different variances σ 2 j . • Linear regression is only an example; more complicated models, as nonlinear regression, can be studied at the cost of the use of some asymptotic results.

Let αj p the least square estimator of α j p , j = 1, . . . , J and p = 1, . . . , P . The estimation of a j k is then given by:

âj k = P p=1 αj p V k,p = V (k)α j . ( 6
)
with Var((â j k ) k=1,...,K ) = σ 2 j V (V V ) -1 V and obviously an estimation or a prediction of the response R k (x) is

Rk (x) = j∈J âj k φ j (x), (7)
with

Var Rk (x) = j∈J φ 2 j (x)σ 2 j V (V V ) -1 V k,k , (8)
where

V (V V ) -1 V k,k is the (k, k)th element of V (V V ) -1 V . Note that Rk (x) is an unbiased estimator of R k (x).

Non simultaneous confidence bands

Each estimation σ2 j of the variance σ 2 j will be considered as exact because of the size of the samples. The construction of a non simultaneous confidence interval for R k (x) relies on the following fact

Z k (x) := Rk (x) -R k (x) J j=1 φ 2 j (x)σ 2 j (V (V V ) -1 V ) k,k
∼ N (0, 1), thanks to the CLT. As a consequence, if z α is the α quantile of the standard Gaussian distribution, we have for each (x, k)

|Z k (x)| < z 1-α/2 ⇔ R k (x) ∈ Rk (x) ± z 1-α/2 J j=1 φ 2 j (x)σ 2 j (V (V V ) -1 V ) k,k .
with probability 1α and we get confidence intervals for R k (x) as x varies which are not simultaneous.

Simultaneous confidence bands

In this section, our aim is to construct simultaneous confidence bands for each individual k i.e. to determine for the given k ∈ K a fractile function F k 1-α (x) such as

P ∀x ∈ X : R k (x) ∈ Rk (x) ± F k 1-α (x) = 1 -α.
Different techniques to build simultaneous confidence bands are encountered in the literature as mentioned in the introduction; e.g. Scheffé method which is not tractable in our framework since we do not require simultaneous confidence bands both in k and x. We propose in this section a new technique based on the supremum of Gaussian random variables, see Section 3.

Under the previous hypothesis and with the same notation, for all k, remind that (Z k (x)) x defines a Gaussian centered and standardized sequence with covariance given by:

Cov Z k (x), Z k (x ) = J j=1 φ j (x)φ j (x )σ 2 j J j=1 φ 2 j (x)σ 2 j J j=1 φ 2 j (x )σ 2 j 4
So if we are able to derive the value S 1-α such that

P sup x=1,...,x |Z k (x)| ≤ S 1-α = 1 -α,
we will get simultaneous confidence bands as required. The distribution of the supremum of a Gaussian process is of great interest but very few exact theoretical results can be found in the literature.

The computation of the maximum of the absolute value of a Gaussian series

In the following Z 1 , . . . , Z n consist of n observations of a centered Gaussian series which is not necessarily stationary but has a constant variance. Without loss of generality, we assume that this variance is 1. We are interested in the distribution of the variables

M * n := sup i=1,...,n |Z i | or M n := sup i=1,...,n Z i .

Classical inequalities

In this section we review classical bounds. The first one developed by Knowles in his 1987 Stanford dissertation is based on the work of Naiman [START_REF] Naiman | Conservative bands in curvilinear regression[END_REF] and Hotelling [START_REF] Hotelling | Tubes and spheres in n-spaces, and a class of statistical problems[END_REF] on the volume of tubes. Define the correlation between X j-1 and X j by ρ j and the length of the sequence (X 1 , . . . , X n ) by L := n j=2 arccos(ρ j ). The length-based bound on P{M n ≥ y} is then

P{M n ≥ y} ≤ 1 -Φ(y) + e -y 2 /2 L 2π (9)
(where Φ is the distribution function of a standard Gaussian variable) and usually beats the crude Bonferroni bound (Miller [19]) which is known to be rough:

P{M n ≥ y} ≤ n(1 -Φ(y)). (10)
The following bound (called W in the sequel) is a natural improvement of the length-based one and have been developed by Efron [START_REF] Efron | The length heuristic for simultaneous hypothesis tests[END_REF]. It still involves the correlation between two successive variables:

P{M n ≥ y} ≤ 1 -Φ(y) + φ(y) n j=1 Φ(yL j /2) -1/2 y/2 ( 11
)
where L j = arccos(ρ j ) and φ is the density function of a standard Gaussian variable.

Finally, Sidak's inequality [START_REF] Sidak | Rectangular confidence regions for the means of multivariate normal distributions[END_REF] (called Sidak in the sequel) applies in any context (in particular when the correlation structure is unknown) and shows that under weak conditions the independent case is the worst case in the sense that

P{M * n ≤ y} ≥ (2Φ(y) -1) n (12)
So using the critical value for an independent sequence will always lead to a conservative test.

For very large n one must use an efficient version of the quantile function of the normal distribution. This is in general easy and it is the case for example in Matlab with the function norminv.

n ≤ (say) 500 or 1000: almost exact calculation by MC-QMC

In such a case the MCQMC Matlab program QSIMVNV written by Genz [START_REF] Genz | Numerical computation of multivariate normal probabilities[END_REF] allows the direct calculation of Gaussian probability over hyper-rectangle for dimension up to 500-1000. It consists of transforming the integral into an integral over the hyper-cube [0, 1] n and then using quasi-Monte-Carlo (QMC) integration with lattice rule. In a final step the procedure is randomized using a Monte-Carlo quasi-Monte-Carlo (MCQMC) method. See Genz [START_REF] Genz | Numerical computation of multivariate normal probabilities[END_REF] or Azaïs and Genz [START_REF] Azaïs | Computation of the distribution of the maximum of stationary Gaussian sequences and processes[END_REF] for more details. The routine QSIMVNV provides an estimation of its numerical error.

Asymptotic band from extreme value

It is well known [START_REF] Leadbetter | Extremes and related properties of random sequences and processes[END_REF] that under very weak assumptions, for stationary sequences,

P a n (M n -b n ) ≤ y → exp(-exp(-y)) as n → +∞ (13) with a n := (2 log n) 1 2 b n := (2 log n) 1 2 -1 2 (2 log n) -1 2 (log log n + log 4π) .
Using symmetry, a bound can be deduced for M * n using the α/2 bound for M n . These results can be easily extended to non-stationary series using the tools of Azaïs and Mercadier [START_REF] Azaïs | Asymptotic Poisson character of extremes in non stationary Gaussian models[END_REF] (In fact this paper considers the more complicated case of random processes). The quality of this approximation is considered in Section 3.4.

Note that in the case where n is very large, the considered levels are very high and some small deviation from normality can heavily affects the results. Remark also that the crude approximation √ 2 log n is often used in practice.

Numerical experiment

The first application considers the comparison of the exact calculation of Genz with Sidak and W for an AR(1) series with parameter 0.9. Here the correlation between two successive variables can easily be determined. Figure 1 compares the 5% critical values as a function of n. It shows clearly that the exact calculation provides an improvement for small values and that improvement decreases with the size. Moreover W is more accurate than Sidak. Table 1 gives three examples. The interpretation of Figure 2 is not easy since the true result is not known and may depend on the particular covariances of the processes. Nevertheless we believe that this dependence is small because the extreme result (Extreme) shows it is true asymptotically. The main points are

• Sidak is better than Extreme since it controls the level and is smaller than Extreme.

• There is some coherence between Extreme, W and Sidak and we believe that for most processes the true result is close to these three. As a consequence the √ 2 log n approximation seems very crude.

• Our conclusion is that for very large n, Sidak and W are the safer choices and give very similar results. In the application section we have chosen Sidak for simplicity.

Application to the EDF load curves 4.1 General settings

As mentioned in the Introduction, we apply the technique presented in the article in the context of EDF load curves prediction. A load curve is a chart showing the amount of electrical energy that a customer uses over the course of time. Power producers like EDF, the electrical French group, use this information to plan how much electricity they need to make available at a given time. But this load curve is only available for customers with automated meter reading. For the others, it must be estimated by EDF in order to know their profitabilities and offer them an appropriate setting of prices. In this EDF context, different approaches have been already studied: a non parametrical approach developed by Misiti et al. [START_REF] Misiti | PREVISOR : un logiciel de prévision non paramétrique[END_REF], a wavelet analysis by Misiti et al. [START_REF] Misiti | Analyse de l'apport des ondelettes pour la prévision des séries chronologiques dans le contexte électrique[END_REF].

To achieve this prediction problem, we used a data set of 917 customers containing for each of them their annual load curve valued every half-hour (which corresponds to 17520 points) and some explanatory variables (17 exactly). These explanatory variables are of different nature: qualitative (activity sector, geographical localization...) and quantitative (total annual consumption, ratio of consumption summer/winter, ratio rush/slack hours...). Customers are in general firms and will be denoted as "individuals" in the sequel. We split data into a learning sample of 832 individuals and a test sample of 85 individuals. These sizes have been chosen heuristically.

For the sake of simplicity, we remove the index k of individual in the sequel.

The decomposed model

In this paragraph, we describe a model, called "decomposed" model, that takes into account the weekly cycles and lies on three main quantities:

• the centered and standardized week load curve (S(h)) h of the given customer;

• the mean of this week load curve (M (w)) w ;

• the standard deviation of this week load curve (σ(w)) w .

More precisely, the decomposed model consists first in constructing a typical week for each of the individuals: for a given individual k and a half-hour h of the week, the value of the typical week at time h is simply the median of all the values of the load curve of the individual k at the same half-hour h of the week. The obtained curve is centered and standardized and the reconstruction is then done by multiplying the typical week by a standard deviation and adding it a mean level. We explain that procedure more precisely in the rest of the section. "Time" is here measured in half an hour.

First of all, let us introduce some notation (recall that k is fixed and thus removed for clarity):

• LC(t) the entire load curve at time t = 1, . . . , T = 17520,

• w(t) = 1, . . . , 53 and h(t) = 1, . . . , 336 are respectively the index of the week and the intra week time associated to Time t. Note that in year 2002, the first and last weeks, with index 1 and 53, are incomplete and contain respectively 6 and 2 days. For these weeks h(t) varies from 1 to 288 and from 1 to 96 respectively.

• C T ot := 17520 t=1 LC(t) the total consumption which is known. Each load curve is decomposed in the following manner • To work on the rescaled profiles of the annual load curves (without dimension), it is normalized by dividing by the total annual consumption. We denote by RS 1 (t) the remaining signal.

• The mean of every week

M (w) = 1 336 w(t)=w RS 1 (t)
(with obvious modification for the first and last weeks) is subtracted

RS 2 (t) = RS 1 (t) -M (w(t))
• Each week is standardized. We compute first the empirical variance of week w σ 2 (w) = 1 336

w(t)=w (RS 2 (t)) 2
and we define

RS 3 (t) = 1 σ(w(t)) RS 2 (t)
• A typical centered and standardized week is computed

S(h) = median{RS 3 (t), h(t) = h}
In fact a small modification of this formula is done to consider bank holidays as Sundays (see [START_REF] Lagnoux | Estimation de courbes de charge et bandes de confiance simultanées[END_REF] for more details).

Eventually the reconstruction is performed in the following way

LC(t) = [S(h(t))σ(w(t)) + M (w(t))] C T ot + ν(t) ( 14
)
for t = 1, . . . , 17520.

Formula ( 14) must be considered as a decomposition of each observed curve LC and not as an estimation in the statistical sense. We assume only that the resulting noise ν, which is due to the replacement of the week by the median, is Gaussian and with constant variance and decaying covariances. In practice, this assumption is approximately satisfied since ν(t) is in general small. This decomposition has shown to be efficient in previous studies at EDF see e.g. [START_REF] Lagnoux | Estimation de courbes de charge et bandes de confiance simultanées[END_REF].

Remark 2

For the functions R k = S(h), σ 2 (w) and M (w), we use a preliminary compression on the Fourier basis. Fourier basis has been chosen because of the periodicities that are encountered in the problem. We choose the coefficients, globally for all individuals, by a hard thresholding procedure to get model [START_REF] Azaïs | Computation of the distribution of the maximum of stationary Gaussian sequences and processes[END_REF]. Since this technique is not the purpose of the current paper, the reader should refer to [START_REF] Donoho | Wavelet Shrinkage: Asymptopia?[END_REF] for details and other references therein.

Estimation and simultaneous confidence bands of the load curve

• To have a comparison scale (which is not the main point of this article), we also consider the model with no decomposition called the global model. In this model, the method is based on a compression of (LC(t)) t :

LC(t) = J j=1 a j φ j (t) + (t) =: LC(t) + (t)
First, LC(t) is estimated by LC(t) following the procedure detailed in Section 2.1. Second, we estimate the variance of (t) by

σ 2 = 1 K K k=1 1 T -1 T t=1 LC k (t) -LC k (t) 2 .
We suppose moreover this estimation exact because it is based on 832 × 17520 observations.

the fact that LC(t) -LC(t)

Var LC(t)

∼ N (0, 1),
where Var LC(t) is the expression (8).

As a consequence for any individual,

LC(t) -LC(t) = LC(t) -LC(t) + (t)
has the variance Var LC(t) + σ 2 . And we are able to construct confidence bands:

⎡ ⎣ LC(t) ± S 1-α j∈J φ 2 j (t)σ 2 j (V (V V ) -1 V ) k,k + σ 2 ⎤ ⎦
where S 1-α is simply 1.96 in the case of non simultaneous confidence intervals and is given by Sidak (Section 3.1, equation ( 12)) in the case of simultaneous confidence bands.

• For the decomposed model, we use the decomposition ( 14) and set

LC(t) = S(h(t))σ(w(t)) + M (w(t)) C T ot + (t).
Now we apply the method of Section 3.2 (MCQMC/Genz method for n = x = 336, 53 and 53 respectively) to every component, namely

S, σ, M

with level α/3, to derive estimation, lower and upper bounds for these quantities. Recall that C Tot is known. We then plug these bounds into decomposition [START_REF] Kendall | Brownian confidence bands on Monte Carlo output[END_REF]to get upper and lower bounds of LC(t) :

LC(t) ∈ LC(t) ± f 1-α (t)
with probability greater than 1α. Then we estimate the variance of

(t) = LC(t) -LC(t)
and we use Sidak to get

| (t)| = |LC(t) -LC(t)| ≤ σ B Sidak ,
where B Sidak is the bound obtained following procedure of Section 3.1.

For any individual, we construct a simultaneous confidence region of the form

| LC(t) -LC(t)| ≤ f 1-α (t) + σ B Sidak .

Prediction quality and covering ratios

To evaluate the accuracy of the predictions, we define an error criterion weighted by the energy cost w. For each individual, we compute a quadratic relative error normalized by the typical price curve of energy:

Err = t w(t) LC predicted (t) -LC real (t) 2 t w(t) LC real (t) 2
Err is the error done for an individual whose real load curve is LC real , predicted curve is LC predicted and w represents a possible weighting (e.g. by prices).

To evaluate the prediction performance on the test sample, we define two covering rates: for the simultaneous confidence bands (SCB), SCR := #{k : the entire curve is in the SCB} #{k} ;

and for the non simultaneous confidence bands (NSCB),

NSCR := #{(k, t) : LC k (t) ∈ NSCB} #{(k, t)} .
Moreover, for ρ ∈ [0, 1], we define SCR(ρ) as the ratio of load curves having at least a percentage ρ of their points in the simultaneous confidence bands. Obviously, with these definitions SCR = SCR(1).

Numerical results

Learning sample

In this section, we present briefly a few results.

Decomposed model

In the study of typical weeks S, hard thresholding suggests us to achieve 99% of the variance which corresponds to keep the 105 best frequencies (i.e. 210 best coefficients). We plot in Figure 3 the typical week of individual 11, its estimation and confidence bands (simultaneous and non simultaneous).

One can find the same type of figure for other individuals in [16].

We determine estimate and confidence bands for the other components in the same way. The construction of the estimates of the entire curve is then done in a natural way using the equation of reconstruction [START_REF] Kendall | Brownian confidence bands on Monte Carlo output[END_REF] given in Section 4.2. Concerning the confidence bands, the operation is more tricky. The three quantities being non negative, we multiply the lower and upper confidence bounds for each of the components studied following the same equation [START_REF] Kendall | Brownian confidence bands on Monte Carlo output[END_REF]. Nevertheless, in that way the level is not conserved; thus in order to have a global confidence level of 1-α, we need to construct for each of the components confidence bands with a confidence level of 1α/3. The same kind of procedure is valid for the simultaneous confidence intervals.

Remark 3

We estimate the 210+53+53=316 parameters (after compression) with 832×17=14144 data (that is about 45 data by parameter) which appears to be completely reasonnable.

Global model

First of all the load curves are normalized by the total consumption, centered and standardized. Remind that in the global model, we study the entire load curves with 17520 points; there is no decomposition.

Compression phase

We proceed as explained in Section 4.2 and the test of significance suggests us to keep 55% of the variance which corresponds to keep the 15 best frequencies (i.e. 30 best coefficients) and to a reduction factor of more than 580.

Regression phase Once the compression has been done, we proceed to the regression phase on the conserved coefficients. Finally, the estimates of the centered standardized load curves are given by:

LC k (t) = j∈J âj k φ j (t).
where J represents the set of the coefficients kept after compression. Since the dimension of the curves is too large leading to the failure of Genz algorithm, the confidence bands are constructed using the asymptotic result given in Section 3.3.

Comments and comparison of the results

We plot the different results in Figures ?? (load curve and estimations) and 5 (load curves and confidence bands) of individual 11 on [200, 600] (for the sake of clarity) for the decomposed and the global models. We remark that the estimation given by the decomposed model seems to be more accurate than the other one and fits better the load curve shape and so do the simultaneous confidence bands. More precisely, if we compare for example the medians of the errors, the decomposed model provides better estimates than the global model with 55% (model suggested by hard thresholding and significance tests techniques). On top of that, the decomposed model is more efficient, easy to compute and last but not least less consuming in terms of cost of simulation.

Test sample

We summarize the different results obtained for the test sample in Table 2 and Table 3. In terms of error, the decomposed model provides better results than the global model (corresponding to a gain about 0.27% for the median error). Concerning the covering ratios, remind that SCR(ρ) represents the ratio of individuals having 100ρ% of their curve in the SCB. For instance, 56 out of the 85 individuals of the test sample have their whole curve contained in the SCB for the decomposed model; 83 have 95% of their curve in the SCB for the global model.

We have to pay somewhere all the assumptions we have made to get our confidence bands. Table 3 shows that if we tolerate the curve to be out of the bound in 1% of the case, we get for the test sample a covering rate 82/85 which is consistent with the nominal value. In our opinion it justifies a posteriori our assumptions.

Conclusion and prospects

We have presented a method which is new in the sense that it is a combination of Sidak's inequality and numerical computation for Gaussian vectors. The method is easy to implement using few line Matlab programs. It has shown to be efficient in our problem of curve prediction for load curves. Moreover there is no doubt that the techniques developed here can be applied to many situations when one wants to predict a set of curves using a basis of functions.

Figure 1 :

 1 Figure1: For the AR(1) series with correlation 0.9, comparison of the 5% critical values computed with Sidak, W and the exact algorithm (from top to bottom) as a function of the length of the series.

Figure 3 :

 3 Figure 3: Estimation of the typical week S of individual 11 and associated confidence bands : the typical week is represented by the thick line, its estimation is the central curve and on both sides of this curve the confidence bands (the closest are the non simultaneous bands and the largest are the simultaneous bands)

Figure 4 :

 4 Figure 4: Comparison between the two models on the interval [200,600] (load curve in thick line, estimations in dashed line for the decomposed model and in solid line for the global model)

Figure 5 :

 5 Figure 5: Comparison between the two models on the interval [200,600] (load curve in thick line, simultaneous confidence bands in dashed line for the decomposed model and in solid line for the global model)

Table 1 :

 1 Numerical experiment

	value of n	10	500 900
	Exact bound 2.56 3.74 3.90
	Sidak	2.80 3.88 4.02
	difference	0.24 0.14 0.12
	W	2.64 3.79 3.93
	difference	0.08 0.05 0.03

Figure 2: For the 5% critical value, comparison of Extreme, Sidak, W (with ρ = 0.8) and 2 log(n) (from top to bottom) for large values of n

Table 2 :

 2 Errors for the test sample

	Models	Median error Mean error ± 5% conf. interval	Error range
	Decomposed model	0.0826	0.1218 ± 0.0260	[0.0077,0.6860]
	Global model	0.1047	0.1335 ± 0.0264	[0.0071,0.7185]

Table 3 :

 3 Covering rates for the test sample

	Model	SCR(1.00) SCR(0.99) SCR(0.95) SCR(0.90)
	Decomposed	56/85	82/85	84/85	85/85
	Global	49/85	73/85	83/85	84/85