Julien Fatome

Esben Bertrand Kibler

Ravn Andresen

Herve Rigneault

Christophe Finot
email: christophe.finot@u-bourgogne.fr

Esben R Andresen

Hervé Rigneault

All-fiber spectral compression of picosecond pulses at telecommunication wavelength enhanced by amplitude shaping

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Optical fibers are an ideal medium to observe the consequences of Kerr nonlinearity. Indeed, thanks to the guided propagation and light confinement, the low third-order nonlinear coefficient of silica is counterbalanced by a large available interaction distance. Therefore, optical fibers have been widely used to observe large self-phase modulation (SPM) of transform-limited pulses, which usually leads to a spectral expansion of the pulse [START_REF] Stolen | Self-phase modulation in silica optical fibers[END_REF]. If this effect has ordinarily a negative impact on long-haul telecommunications or high-power fiber amplifiers, various specific applications may on the contrary benefit from this spectral enhancement. For example, highly-coherent multi-wavelength continuum sources have been proposed spanning the whole band of conventional optical telecommunications [START_REF] Yousoff | A 36channel x 10-GHz spectrally sliced pulse source based on supercontinuum generation in normally dispersive highly nonlinear holey fiber[END_REF]. Another widely used application is nonlinear temporal pulse compression and relies on the linear compensation of the nonlinear chirp induced by SPM. Thanks to this technique, femtosecond few-cycles pulses have been generated [START_REF] Tomlinson | Compression of optical pulses chirped by self-phase modulation in fibers[END_REF]. However, the quality of the results depends directly on the temporal intensity profile of the initial pulses. Therefore, it has been shown that initial pulse shaping may significantly improve the performance of these devices [START_REF] Parmigiani | Ultra-flat SPM-broadened spectra in a highly nonlinear fiber using parabolic pulses formed in a fiber Bragg grating[END_REF][START_REF] Clarke | Reconfigurable Optical Pulse Generator Employing a Fourier-Domain Programmable Optical Processor[END_REF][START_REF] Schimpf | Controllong the influence of SPM in fiber-based chirped-pulse amplification systems by using an actively shaped parabolic spectrum[END_REF] by taking advantage of the parabolic pulse features [START_REF] Finot | Optical parabolic pulse generation and applications[END_REF].

A less-known potential application of SPM (even if it was described as soon as 1978 [START_REF] Stolen | Self-phase modulation in silica optical fibers[END_REF]) relies on the configuration where the initial pulse is negatively linearly chirped corresponding to an anomalous dispersive medium. In such a configuration, the SPM-induced chirp cancels the initial frequency modulation and long and short wavelengths redistribute toward the center wavelength, leading to a spectral narrowing of the pulse [START_REF] Markaryan | Spectral compression of ultrashort laser pulses[END_REF][START_REF] Planas | Spectral narrowing in the propagation of chirped pulses in single-mode fibers[END_REF][START_REF] Oberthaler | Spectral narrowing of ultrashort laser pulses by selfphase modulation in optical fibers[END_REF][START_REF] Andresen | Spectral compression of femtosecond pulses in photonic crystal fibers[END_REF]. Concentrating the energy in the central part of the spectrum has for example been found of practical interest for the generation of high power picosecond pulses from femtosecond oscillators [START_REF] Limpert | SPM-induced spectral compression of picosecond pulses in a single-mode Yb-doped fiber amplifier[END_REF] or for biophotonics applications with an enhanced spectral resolution [START_REF] Andresen | Tunable light source for coherent anti-Stockes Raman scattering microspectroscopy based on the soliton self-frequency shift[END_REF]. However, once again, the performances of the spectral compression process crucially depends on the initial pulse shape and routinely used hyperbolic secant or Gaussian temporal intensity profiles lead to the presence of large spectral pedestals. To overcome this issue, the optimum pulse shape has been found to be theoretically parabolic [START_REF] Finot | Spectral compression of optical parabolic similaritons[END_REF][START_REF] Finot | Parabolic pulse evolution in normally dispersive fiber amplifiers preceding the similariton formation regime[END_REF], which has recently been experimentally validated by using femtosecond pulses evolving in a passive photonic crystal fiber around 1 µm [START_REF] Andresen | Transform-limited spectral compression by self-phase modulation of amplitude shaped pulses with negative chirp[END_REF].

In this contribution, we provide a new illustration of the benefits brought by an initial shaping of the temporal intensity profile. Our experiments explore a completely different range of input parameters (initial wavelength, duration and peak power as well as typical fiber lengths) and a high-quality spectral compression by a factor above 10 is demonstrated in a passive fiber.

We first recall the basis of spectral narrowing and quantify the process by means of numerical simulations. We then describe our all-fiber experimental setup and discuss the experimental results that have been obtained.

Principle of operation and numerical simulations

Spectral narrowing occurs when a linearly chirped pulse (with a chirp provided by an anomalous dispersion medium such as a pair of gratings [START_REF] Markaryan | Spectral compression of ultrashort laser pulses[END_REF][START_REF] Limpert | SPM-induced spectral compression of picosecond pulses in a single-mode Yb-doped fiber amplifier[END_REF], a prism [START_REF] Oberthaler | Spectral narrowing of ultrashort laser pulses by selfphase modulation in optical fibers[END_REF][START_REF] Andresen | Spectral compression of femtosecond pulses in photonic crystal fibers[END_REF] or an hollow-core photonic bandgap fiber [START_REF] Rusu | All-fiber picosecond laser source based on nonlinear spectral compression[END_REF]) is affected by self-phase modulation in a highly nonlinear Kerr medium.

Basically, if the initial pulse has the temporal complex field envelope ψ given by: 2 ( , 0) ( , 0) exp 2

T C T P T i T ψ   = -     (1)
with P being the temporal intensity profile and C T the chirp parameter, then the SPM induced after a propagation length L is given by the following nonlinear term that adds to the initial phase, γ is respectively the nonlinear coefficient of the fiber under test :

( , ) ( , 0) NL T L P T L ϕ γ = (2)
This results in the following pulse (here the impact of other effects such as dispersion or attenuation has been neglected) :

2 ( , ) ( , 0) exp ( , 0) 2 T C T L P T i T P T L ψ γ     = - -         (3)
From this simple equation, it is clear that starting from a linearly chirped pulse, one can obtain a transform-limited pulse of the same duration (therefore with a narrower spectrum) if this condition is fulfilled:

2 ( , 0) 2 T C T P T L K γ - = (4)
K being a constant. If this requirement is not fulfilled, the input chirp compensation will only be partial with a time-dependent term, resulting in the large spectral pedestals that are usually observed [START_REF] Oberthaler | Spectral narrowing of ultrashort laser pulses by selfphase modulation in optical fibers[END_REF][START_REF] Limpert | SPM-induced spectral compression of picosecond pulses in a single-mode Yb-doped fiber amplifier[END_REF][START_REF] Rusu | All-fiber picosecond laser source based on nonlinear spectral compression[END_REF].

Following Eq. ( 4), it is important to note that this last condition is perfectly fulfilled for an input temporal intensity profile having a parabolic shape : 2 ( , 0) 2

T C P T T K L γ = - (5)
that can be written using the usual notation for a parabolic bright pulse as (for |T| < T P ) : 2 2 ( , 0) 1

P P T P T P T   = -     (6)
with P P = -K being the peak power of the pulse and T P being the half-total width of the parabolic pulse such that 2 2

P P T L P T C γ = - (7)
We can make out that K should be negative and that a negative input chirp coefficient C T is needed. From ( 7) it also appears that the peak power of the pulses is directly proportional to the chirp coefficient and inversely proportional to the length of the fiber and to the nonlinear coefficient. Note that using dips of light on a continuous background could also be a possible solution and dark parabolic profiles [START_REF] Hirooka | Bright and dark 40 GHz parabolic pulse generation using a picosecond optical pulse train and an arrayed waveguide grating[END_REF] can also fulfill Eq. ( 5) to spectrally compress positively chirped pulse.

In order to get more precise insights into the nonlinear dynamics of the spectral compression, we numerically integrate the standard nonlinear Schrödinger equation that takes into account also the second order dispersion of the fiber β 2 :

2 2 2 2 2 T i z β ψ ψ γ ψ ψ ∂ ∂ ∂ = - ∂ (8)
In the following, we base our discussion on the experimental parameters that will be further detailed in section 3. We will compare hyperbolic secant and parabolic pulses having a full width at half-maximum (FWHM) temporal duration of 70 ps and a linear chirp coefficient C T of -1.93 GHz/ps. Both pulses have a similar FWHM spectral width of 130 GHz. The spectral compression will take place in an optical fiber having the following properties: β 2 = 9.6 10 -4 ps 2 /m, γ = 10 /W/km and a fixed length L of 1 km.

Figure 1 describes the evolution of the output spectral profiles as a function of the input pulse peak-power for the two pulse shapes under investigation (hyperbolic secant and parabolic pulses). This figure outlines the significant influence of the initial shape on the compression behavior. If the spectral width and the power required for optimum compression seem comparable in both cases, the presence of sidelobes and oscillations is clearly accentuated in the case of a sech 2 profile. One can also make out that the spectral brightness is enhanced using a parabolic input pulse. The spectral expansion following the point of optimum compression also differs significantly. Figure 2 describes more quantitatively the differences between the two initial pulse shapes. More precisely, we have plotted the evolution of two quantities as a function of the input peak-power: the FWHM of the optical spectrum (Fig. 2a) and the Strehl ratio (Fig. 2b) that we use in order to assess the quality of the spectral compression. This Strehl ratio, often used to assess the spatial performance of imaging system, can also be used to quantify the quality of a temporal compression [START_REF] Zaouter | Third-order spectral phase compensation in parabolic pulse compression[END_REF]. In our case, we define the Strehl ratio as the ratio between the actual spectral peak power and the spectral peak power obtained assuming a flat temporal phase of the pulse. Therefore, the Strehl ratio is comprised between 0 and 1, 1 defining an ideal spectral compression. From the results of these two figures, we can make out several important conclusions. First, we can note that the optimum compression (both in terms of spectral width and Strehl ratio) is obtained for an input peak power of 1.55 W for sech and parabolic pulse shapes, respectively. This value corresponds to input pulse energies of 105 and 125 pJ respectively and is in agreement with the prediction of Eq. ( 8). In both cases, a compression factor close to 13 is achieved and one can note that between 1.2 W and 1.8 W of input peakpower, the FWHM spectral width does not exhibit major differences. Therefore, the main parameter to definitively assess the best configuration and power will be the Strehl ratio. A significantly higher Strehl ratio is indeed obtained in the case of parabolic pulses, with value as high as 0.99, synonym of transform-limited pulses. It has to be compared to the performances achieved with a picosecond sech pulse with a poor resulting Strehl ratio of 0.24.

Interestingly, in the case of a parabolic pulse, the root-mean-square (rms) width of the pulse spectrum according to the input peak power follows a simple parabolic law. This result fully agrees with the evolution of the temporal rms width in a dispersive medium : in [START_REF] Anderson | Analytic study of pulse broadening in dispersive optical fibers[END_REF], it has indeed been shown that the temporal rms of any pulse propagating in a dispersive medium (i.e. a medium that imposes a spectral parabolic phase) scales following a non-trivial parabolic law. In the case of a parabolic pulse, we can clearly draw a parallel between the impact of the temporal parabolic SPM-induced phase and the action of second order dispersion acting as a parabolic phase in the spectral domain [START_REF] Finot | Parabolic pulse evolution in normally dispersive fiber amplifiers preceding the similariton formation regime[END_REF]. Details of the output temporal and spectral intensity and chirp profiles generated at the optimum compression point are shown in Fig. 3. The temporal features of the spectrally compressed pulse confirm that for an initial parabolic shape, a null temporal chirp is obtained, synonym of the expected flat phase. On the contrary, a non-negligible residual chirp is visible for the hyperbolic secant pulse. Comparing the initial and compressed pulses, we can note that there is no change in the temporal intensity profile, thus confirming that the impact of chromatic dispersion is here negligible. Regarding the details of the compressed spectral profile (Fig. 3c), as already described in [START_REF] Finot | Parabolic pulse evolution in normally dispersive fiber amplifiers preceding the similariton formation regime[END_REF], the parabolic shape leads to significantly lower sidelobes, with a very low ratio of energy contained in the wings of the spectrum. The recompressed spectral profile can be theoretically described by a Bessel function of the first kind [START_REF] Finot | Parabolic pulse evolution in normally dispersive fiber amplifiers preceding the similariton formation regime[END_REF]. In the case of the hyperbolic pulse, some noticeable sidelobes are visible containing a non-negligible part of the pulse energy. Consequently, the spectral peak intensity is twofold under the one achieved with a parabolic pulse.

Experiments

Experimental setup

The experimental setup is sketched in Fig. 4. An initial picosecond pulse train is delivered by a modelocked erbium doped fiber laser running at a repetition rate of 10 GHz at a central wavelength of 1552.5 nm. The initial pulses are transform-limited with an hyperbolic secant intensity profile typical of soliton pulses and with a FWHM temporal duration of 2.4 ps. The repetition rate is then decreased down to 1.25 GHz by means of an intensity modulator driven by a 1:7 duty cycle pulse sequence provided by a 10-Gbit/s electrical pattern generator. Decreasing the repetition-rate of the pulse train enables us to have access, for the same average power, to a higher peak-power. Note that this pulse picking technique has also the additional advantage to avoid any spectrum comb structure of 10-GHz spaced lines. Indeed, with a 1.25 GHz repetition rate and given our optical spectrum analyzer resolution, the resulting spectrum appears as continuous. A phase modulator driven by a sinusoidal beating at a frequency around 100 MHz is also used to avoid the deleterious effects of Brillouin backscattering. The amplitude shaping is achieved by a liquid crystal on silicon (LCoS) device [START_REF] Clarke | Reconfigurable Optical Pulse Generator Employing a Fourier-Domain Programmable Optical Processor[END_REF] that enables to tailor both the amplitude and phase spectral profiles of the incoming pulse train. Note however that in our case, we will only benefit from the amplitude shaping capabilities of the devices. Indeed, the linear spectral chirp that is targeted is too high to be directly synthesized by the LCoS device. Consequently, in order to imprint the linear chirp, we use a piece of anomalous standard single mode fiber with a length L SMF of 3900 m and a second-order dispersion β 2SMF = -21.9 ps 2 /km. The fiber rigorously leads to a parabolic spectral phase β 2SMF L SMF ω 2 / 2, ω being the radial frequency. Due to the high level of parabolic spectral phase, this leads in the temporal domain to the expected nearly quadratic temporal phase. The initial 2.4-ps pulses are thus stretched to a temporal duration of 70 ps. The resulting pulses are characterized in intensity by a fast oscilloscope with a 63 GHz bandwidth, high enough to obtain the shape of the intensity profile. The low incoming power in the SMF and its low nonlinear coefficient ensure that no unwanted nonlinear effects affect the pulse propagation in this segment. Then the pulses are amplified to the required average power (up to 30 dBm) by a high-power erbium doped fiber amplifier (EDFA) and the output power is varied by an optical variable attenuator in order to avoid power-dependent nonlinear effects in the amplifier. The optical attenuator is controlled by a computer (PC) and synchronized to the optical spectrum analyzer, which enables us to obtain systematic mapping of the output spectral shape. Finally, the spectral narrowing takes place in a highly nonlinear fiber (commercially available HNLF distributed by OFS) with a length of 1 km, an attenuation of 1 dB/km and a nonlinear coefficient of 10/W/km. The chromatic dispersion is chosen to be normal (β 2 = 9.6 10 -4 ps 2 /m is used) in order to prevent from any modulation instability process [START_REF] Andresen | Spectral compression of femtosecond pulses in photonic crystal fibers[END_REF] . An optical spectrum analyzer (resolution of 5 GHz) has been used to record the resulting spectra.

Experimental results

The experimental temporal intensity profile recorded after pulse shaping and stretching is plotted in Fig. 5. We can clearly see that the initial intensity spectral shaping enables to efficiently shape the hyperbolic secant pulses into a parabolic profile. The impact of the power launched into the fiber is demonstrated in Fig. 6 where we have plotted similarly to Fig. 2, the output spectrum according to the input peak power of the incoming pulses. This map is made of 200 spectra.

Those experimental results agree well with the qualitative trends that have been described in the previous section: the spectral narrowing is progressive and becomes optimum for input peakpowers close to 1.5 W. The significant difference in behavior previously discussed between the two types of pulse shape is also confirmed.

Fig. 5 : Initial pulse profiles after stretching in the 3.9-km long SMF segment. In absence of amplitude spectral shaping, the recorded intensity profile (grey line) is compared to a sech 2 intensity fit (grey circles). Using an initial spectral shaping, the experimental results (black line) are compared to a parabolic intensity profile (black circles). Fig. 6 : Experimental results: Evolution of the spectral profile as a function of the input peak power for an initial sech 2 intensity profile (a) and for a parabolic intensity profile (b). Spectra are normalized so that the energy of the spectra equals 1.

More details on the spectral shape of the recompressed pulses are given in Fig. 7a and stress that the parabolic shaped pulses lead to an enhanced spectral compression providing lower sidelobes and increased peak spectral power. Compared to input pulses, the power density has then been increased by a factor 9. Moreover, in agreement with numerical predictions, an improvement of spectral brightness by a factor 2.1 is achieved compared to the sech 2 pulses. A compression factor ratio of 12.2 is also measured, consistent with the factor expected from the numerical simulations. Note however that this value is slightly impacted by the resolution of our OSA (we believe the real value of the compressed spectral width to be slightly lower). We have checked that the temporal intensity profiles of the output pulses were not affected by the compression process. Therefore, the resulting time-bandwidth product is found to be 0.745, to be compared to the value of transform-limited parabolic pulse of 0.66. We are thus only 12% higher than the bandwidth limited case. On the other side, for hyperbolic secant pulses, the resulting time-bandwidth product is more than twofold above the transform limit condition (0.315 for a transform-limited hyperbolic secant pulse).

Regarding the evolution of the pulse quality, the Strehl ratio demonstrates once again the superiority of the parabolic-shaped pulses with a Strehl ratio 3 times higher than in the sech pulse case. Finally, the evolution of the rms spectral width measured for the parabolic pulse shape as a function of input peak power confirms the parabolic law described in the previous section.

Discussion

In our experiment, we have demonstrated a x12 compression ratio. To the best of our knowledge, this result represents a record for generation of near transform limited pulse: the previous works dealing with the same technique in a passive fiber have reported a compression factor below 10 [START_REF] Planas | Spectral narrowing in the propagation of chirped pulses in single-mode fibers[END_REF][START_REF] Oberthaler | Spectral narrowing of ultrashort laser pulses by selfphase modulation in optical fibers[END_REF][START_REF] Limpert | SPM-induced spectral compression of picosecond pulses in a single-mode Yb-doped fiber amplifier[END_REF][START_REF] Andresen | Tunable light source for coherent anti-Stockes Raman scattering microspectroscopy based on the soliton self-frequency shift[END_REF][START_REF] Andresen | Transform-limited spectral compression by self-phase modulation of amplitude shaped pulses with negative chirp[END_REF][START_REF] Rusu | All-fiber picosecond laser source based on nonlinear spectral compression[END_REF], except in [START_REF] Andresen | Spectral compression of femtosecond pulses in photonic crystal fibers[END_REF] but the high level of compression ratio was achieved at the expense of high sidelobes and pulses were far from being close to transform-limit. Regarding many aspects, the superiority of the shaped pulse has been definitively confirmed. However, the results have been below the theoretical predictions obtained for a parabolic pulse shape. In particular, higher substructures have been recorded, thus explaining the lower Strehl ratio equal to 0.57 instead of 1. We attribute this mainly to the imperfect parabolic pulse shape that has been initially synthesized. Indeed the process is intrinsically sensitive to the gradient of the intensity profile and as can be seen in Fig. 5, some deviations from a perfect inverted parabola are observed. This temporal defect could have been suppressed using a feedback loop in the shaping process to iteratively optimize the input shape. This point is however beyond the scope of the present paper, our target being to confirm the benefits of parabolic pulses.

Higher compression factors could have been obtained using a longer stretching fiber.

Given, the high ratio of β 2SMF /β 3SMF in SMF fiber, the influence of third-order dispersion β 3SMF can be neglected as a first approximation (the typical dispersion length linked to β 3SMF is more than one order of magnitude higher than the fiber length under use) : for the picosecond pulse under investigation, the cubic phase induced by β 3SMF is not as deleterious as what has been outlined when dealing with femtosecond pulses [START_REF] Oberthaler | Spectral narrowing of ultrashort laser pulses by selfphase modulation in optical fibers[END_REF][START_REF] Kutuzyan | Dispersive regime of spectral compression[END_REF]. If really needed, they can be compensated for using the phase tailoring ability of the LCoS device. Second-and third-order dispersion of the dispersion-flattened HNLF have also not been found to be here a limiting factor during the spectral compression process since dispersive effects have low influence on the narrow spectrum.

However, some practical issues have to be taken into account. First, the strong temporal broadening of the input pulses should not lead to pulse overlapping in the wings of the pulses.

Therefore, around half of the bit time seems to be an upper limit to the stretching that should be achieved. Then, increasing the temporal pulse duration according to (7) requires a higher peak power and therefore a higher average power. A first limitation is of course the power capacity of the EDFA. But another stringent effect is the Brillouin backscattering [START_REF] Smith | Optical power handling capacity of low loss optical fibers as determined by stimulated Raman and Brillouin scattering[END_REF]. Indeed, if Brillouin backscattering is often fully negligible for optical spectra of several tens of GHz, in our configuration, due to the compression process that strongly increases the power spectral density, Brillouin backscattering becomes a stringent issue that we have had to overcome by means of an additional phase modulation.

Conclusion

We have demonstrated the spectral compression of picosecond pulses in an all-fiber configuration at telecommunication wavelengths. Using widely available standard single mode fiber as a stretcher to impose the initial linear chirp, a spectral narrowing by a factor above 12 has been achieved for output pulses close to the Fourier transform limit. This compression has been based on a kilometer long highly nonlinear fiber with W-level peak powers and has involved pulse energies well below the nJ level and is compatible with GHz operation. An initial parabolic pulse shaping has also been shown to significantly improve the results. The quality of the compressed spectra has indeed been enhanced with much lower substructures and an improved Strehl ratio.

Fig. 1 :

 1 Fig.1: Simulation: Evolution of the output spectral profile according to the input peak power for an initial sech 2 intensity profile (a) and for a parabolic intensity profile (b). Spectra are normalized so that the energy of the spectra equals 1.

Fig. 2 :

 2 Fig. 2 : Simulation: Evolution of the output spectral parameters according to the initial peak-power for the parabolic pulse shape (solid black line) and sech pulse ( grey line) (a) FWHM spectral width. (b) Strehl ratio. (c) Evolution of the rms spectral width as a function of input peak power for a parabolic pulse : numerical results (solid black line) are compared to a parabolic fit (dots).

Fig. 3 :

 3 Fig. 3 : Simulation: (a) Temporal intensity and chirp profiles obtained for the optimum point of compression for hyperbolic secant and parabolic profiles (grey and black line respectively). The input temporal intensity profiles are plotted with full circles. (b) Spectral intensity profiles at the point of optimum compression for hyperbolic secant and parabolic profiles (grey and black solid lines respectively). The results are compared to the initial spectral plotted in dashed lines.

Fig. 4 :

 4 Fig. 4 : Experimental setup. MLL : modelocked laser ; IM : intensity modulator ; PM : phase modulator ; SMF : single mode fiber ; EDFA : erbium-doped fiber amplifier ; OVA : optical variable attenuator ; HNLF : highly nonlinear fiber ; OSA : optical spectrum analyzer.

Fig. 7 :

 7 Fig. 7 : (a) Optical Spectra obtained at the optimum compression point (solid line) for a sech pulse (grey line) and parabolic shaped pulse (black line) and compared to the input spectra (dashed line). (b) Evolution of the Strehl ratio according to the initial peak-power for the parabolic pulse (solid black line) and sech pulse (solid grey line) (c) Evolution of the rms spectral width according to the input peak power for a parabolic pulse : experimental results (solid black line) are compared to a parabolic law.

Acknowledgements

We acknowledge the support of the Agence Nationale de la Recherche (ANR SOFICARS, grant ANR-07-RIB-013-03) and Conseil Regional de Bourgogne (PARI Photcom). We also thank Philippe Morin, Kamal Hammani, Guy Millot, Charles-Henri Hage and Benoit Frisquet for their help and for fruitful discussions.


S Castric

L Denis-Vidal
email: lilianne.denis-vidal@math.univ-lille1.fr

Z Cherfi
email: zohra.cherfi@utc.fr

G Joly Blanchard

N Boudaoud
email: nassim.boudaoud@utc.fr

Modeling Pollutant Emissions of Diesel Engine based on Kriging Models: a Comparison between Geostatistic and Gaussian Process Approach

Keywords: Modelling, Optimization problems, modelling errors, Automotive emissions, Diesel engines

In order to optimize the performance of a diesel engine subject to legislative constraints on pollutant emissions, it is necessary to improve their design, and to identify how design parameters affect engine behaviours. One specificity of this work is that it does not exist a physical model of engine behaviour under all possible operational conditions. A powerful strategy for engine modeling is to build a fast emulator based on carefully chosen observations, made according to an experimental design. In this paper, two Kriging models are considered. One is based on a geostatistical approach and the other corresponds to a Gaussian process metamodel approach. Our aim is to show that the use of each of these methods does not lead to the same results, particularly when "atypical" points are present in our database. In a more precise way, the statistical approach allows us to obtain a good quality modeling even if atypical data are present, while this situation leads to a bad quality of the modeling by the geostatistical approach. This behaviour takes a fundamental importance for the problem of the pollutant emissions, because the analysis of these atypical data, which are rarely erroneous data, can supply precious information for the engine tuning in the design stage.

INTRODUCTION

The automotive industry faces the competing goals of producing better performing vehicles and keeping development time with low costs. It is crucial for the manufacturers to be able to produce fuel-economic vehicles, which respect pollutant emissions standards, and which meet the customersexpectations. Accordingly, the complexity of the engines responses to be optimized and the number of the parameters to be controlled during the design stage, have increased rapidly, in the last years. In order to deliver vehicles, which respond to these requirements, in a reasonable time scale, companies use design of experiments (DOE) in one side, and modeling, in the other side. DOE [START_REF] Box | Hunter Statistics for Experimenters: Design, Innovation, and Discovery 2nd Edition[END_REF]] is a powerful tool, but the cost of the experiments and their duration, particularly in the field of pollutant emissions, can be a limit to their use in automotive industry. The engine developers use two main approaches to model engine behaviour. The first one is based on chemical and physical models, via differential systems, see Borg [START_REF] Borg Jonathan | Alex On the application of Wiebe functions to simulate normal and knoking spark-ignition combustion[END_REF] et [START_REF] Brahmi | Talon Identifiabilité et estimation des paramètres du modèle de combustion d'un moteur diesel CD[END_REF]. This approach is not the subject of this article, because in our case, we don't have such models. Further-more, even when these models are available, generally, they are time-consuming, impractical for multi-objective optimisation routines, and fail to capture all the trends in the engine system described by measured data. Statistical modeling based on carefully chosen measured data of engine performance, according to an experimental design, is an important alternative technique. Strategies based on Lolimot [START_REF] Castric | Schimmerling Diesel engine combustion model for tuning process and a calibration method[END_REF]] (Local Linear Model Tree) and Zeldovich mechanisms [START_REF] Heywood | Internal combustion engine fundamentals[END_REF] have been developed in order to predict emissions of NOx. In the first case, the corresponding model can lead to singular points, which reduces the precision of the results. In the second case, the results are not satisfactory enough. Therefore, several statistical modeling methods offer interesting possibilities, see [START_REF] Edwards | The optimization of common rail FIE equipped engines through the use of statistical experimental design, mathematical modelling and genetic algorithms[END_REF]. The methods based on statistical training such as the neural networks [START_REF] Fang | Sudjanto Design and modeling for computer experiments Chapman and Hall/CRC[END_REF], knew a striking success. At the same time, the use of the computer modeling [START_REF] Santner | The design and analysis of computer experiments[END_REF] and [START_REF] Bayarri | [END_REF]] is strongly developed. The methods based on the surface responses methodologies are studied and compared, particularly by [START_REF] Jin | Simpson Comparative studies of metamodeling techniques under multiple modeling criteria[END_REF], taking into account various levels of complexity of the problem in term of nonlinearity and of the number of variables. However, most of these approaches cannot be used in our case, because we are drastically limited by the small number of experiments which our industrial partner is able to do. The advantages offered by the kriging methods [Matheron (1963), Zhu et al. (2010) andJack (2009)] brought us to choose them for our study. Two approaches of Kriging model are considered in this work. The first one is based on a geostatistical approach. Softwares such as R and Matlab contain a toolbox for the use of this method, but unfortunately, they are restricted to less than 3 dimensions. The method has been adapted to higher dimensions, and more generally an innovative approach for functional data has been proposed. The performance of this model is due to its ability to take into account the spatial dependence of data [Cressie (1991)], and has minimal variance estimators without bias. But this method may be difficult to couple with optimization process of calibration. The second one is a Gaussian process metamodel approach [START_REF] Baillargeon | Fitzback interpolation statistique multivariable de donnes de prcipitations dans un cadre de modlisation hydrologique Colloque Géomatique 2004: un choix stratégique[END_REF]. In this case, the system response in consideration is modeled by the sum of two quantities f (x) and δ(x), where f (x) corresponds to a deterministic function which is a linear combination of known functions (in our case a one-degree polynomial model) and δ(x) is the stochastic part and corresponds to a Gaussian stationary process of order 2, see [START_REF] Iooss | Contribution au traitement des incertitudes en modélisation numérique: propagation d'ondes en milieu aléatoire et analyse statistique d'expériences simulées HDR[END_REF]. The aim of this paper is to bring response to the problem of the polluting emissions. At first, we wish to know if the two kriging models are adapted to the modeling of the three pollutants studied here: NOx, smokes and CO. Then, the question which arises is to know if, according to the values taken by the coefficient used to estimate the quality of the modeling, conclusions can be established as for the presence or not of atypical points. The detection of these atypical points and their interpretation are very important, because these points are often revealing specific information for engine tuning. This paper is organized as follows: In the second section, the Kriging techniques, geostatistical approach and Gaussian process metamodel, are recalled. In the third section, the engine behaviour and the importance of controlling pollutant emissions are described. In the last section, numerical results and some elements of discussion are given.

KRIGING TECHNIQUES

Kriging methods are used frequently for spatial interpolation of soil properties, see [START_REF] Krige | A statistical approach to some basic mine valuation problems on the Witwatersrand[END_REF].Kriging is a linear least squares estimation algorithm, which is used for interpolation. The aim is to estimate the value of an unknown real function y at point x * , given the values of the function y at some other points x (i) ∈ IR for each i ∈ {1, ..., n}. The Kriging estimator is defined by

ŷ(x * ) = n i=1 λ i y(x (i) )
(1)

In this paper, we will note y(x

(i) ) = y (i)
Where n is the number of surrounding observations x (i) is a vecteur composed of the d values (x

(i) 1 , ..., x (i)
d ) of the factors at point i λ i is the weight corresponding to the observation y(x (i) ). The weights are estimated in order to make the estimator unbiased with minimal variance.

The system response is treated as a realization of a random function y(x). This model can be written as:

y(x) = f (x) + δ(x) .
The deterministic function f (x) provides the mean of y(x) and δ(x) is the stochastic part which verifies some assumptions depending on the choice of Kriging model.

Geostatistical approach: Ordinary Kriging

In this approach the weights are determined such that the following Kriging variance var(ŷ(x * ) -y(x * )) is minimal under the unbiased constraint given by

n i=1 λ i = 1 .
(2)

It leads to a classical optimization problem with equality constraint. The Lagrange multiplier theory is used in order to work out this problem. This leads to a linear system to be solved, see [START_REF] Davis | Statistics and Data Analysis in Geology 2nd[END_REF]. This system produces, under certain assumptions specified below, a function called variogram described in the following paragraph.

Variogram The variogram is a function representing the spatial dependency. It is obtained from the stationarity definition. Indeed, this stationarity hypothesis is an indispensable condition for the use of the Kriging method. In the case of ordinary Kriging, the expression of the variogram is obtained from the following definition of intrinsic stationarity:

E(y(x (i) + h) -y(x (i) )) = 0 , 1 ≤ i ≤ n var(y(x (i) + h) -y(x (i) )) = 2γ(h) , 1 ≤ i ≤ n .
More precisely, the expression of the theoretical variogram is deduced from the second condition of intrinsic stationarity. With the assumption of isotropy, the variation of a data set will be dependent only on distance ||h|| between two locations. To infer the variogram from observed data, we will then use the common formula for the experimental variogram, see Cressie (1991):

γ(r) = 1 |N (r)| N (r) y(x (i) ) -y(x (j) ) 2 (3)
where

N (r) = (i, j), ||x (i) -x (j) || = r .
The quantity |N (r)| is the pair number of N (r) and the function γ(r) is the experimental variogram.

Variogram Model

The experimental variogram presented in equation ( 3) estimates the theoretical variogram, for only a finite number of distances. Moreover, it does not necessarily form a valid variogram. This means that maybe, it does not concern a negative conditionally function. Indeed, this condition is necessary to ensure the positivity of the variance of a sum of random variables, see [START_REF] Christakos | On the problem of permissible covariance and variogram models[END_REF]. The experimental variogram is then modelled by a function of negative conditional type and is defined for all distances. It is named the variagrophic model. This model must be selected among the various forms of the variogram models, which exist in the literature and adjusted to the experimental variogram [START_REF] Arnaud | Estimation et interpolation spatiale Hermes Science Publications[END_REF]. It means that the parameters of the model must be estimated. This adjustment is done with an estimation method such as the weighted least squares or maximum likelihood method. Once the variographic model is chosen, and its parameters estimated, the weights λ i which appear in (1) are computed by solving the following system:

AΛ = B with A =      γ(r 11 ) γ(r 12 ) . . . γ(r 1n ) 1 γ(r 21 ) γ(r 22 ) . . . γ(r 2n ) 1 . . . . . . . . . . . . . . . γ(r n1 ) γ(r n2 ) . . . γ(r nn ) 1 1 1 1 1 0      Λ = [λ 1 λ 2 . . . λ n λ] t B = [γ(r 01 ) γ(r 02 ) . . . γ(r 0n ) 1] t and r ij = ||x (i) -x (j) || , r 0j = ||x * -x (j)
|| where λ is the Lagrange multiplier. The function γ(r ij ) is the variogram model used for adjusting the experimental variogram.

Kriging Emulator Validation

The true test of the quality of the fitted emulator model is its ability to predict the response at untried factor values. In order to maximally exploit the data to aid model fitting, the emulators are validated using leave-one-out cross validation. This process involves taking the fitted model and re-fitting it to a subset of non used experimental data. More precisely, for an experiment with d design factors, the set of n experimental design points and corresponding responses, contain the information used to build the Kriging model. A cross validation involves predicting at each design point in turn when that point is left out of the predictor equations. Let ŷ(x (i) ) be the estimate of the y(x (i) ) based on all the design points except x (i) . The prediction error is then calculated by a coefficient of determination R 2 used in multiple regression:

R 2 = 1 - n i=1 |y(x (i) ) -ŷ(x (i) )| 2 n i=1 |y(x (i) ) -y| 2 (4)

Gaussian Process Metamodel

In this approach, the deterministic function f (x) is a linear combination of known functions In this study, a one-degree polynomial model is used and f(x) can be rewritten as follows:

f β (x) = β 0 + d 1 β i x i
where β = (β 0 , . . . , β d ) t is the regression parameter vector. The stochastic part δ(x) is a Gaussian, stationary of order two process, then:

E(δ(x)) = 0 E(δ(x), δ(z)) = σ 2 δ R(x -z) where σ 2
δ denotes the variance of Y and R is the correlation function. Among correlation models given in the literature, see Abr (1997), we choose to use the generalized exponential correlation function:

R θ,q (x -z) = exp(- d k=1 θ k |x k -z k | q k ) or R θ,q (x -z) = d k=1 exp(-θ k |x k -z k | q k )
where θ = (θ 1 , . . . , θ d ) t and q = (q 1 , . . . , q d ) t are the correlation parameters (hyperparameters) with θ k ≥ 0 and 0 ≤ q k ≤ 2, ∀k ∈ {1, . . . , d}. This choice is motivated by wide spectrum of shapes that such a function offers. If a new point x * = (x * 1 , . . . , x * d ) is considered, the following unbiased linear predictor is obtained:

ŷ = f β (x * ) + r(x * )Γ -1 (Y -F )
where

• Y is the vector of observations Y = (y (1) , . . . , y (n) ) t , • F = (f β (x (1) ), . . . , f β (x (n) )) t , • Γ = σ 2 δ R θ,q (x (i) -x (j) ) (i,j)∈{1,...,n} 2 is the covari- ance matrix of δ n = (δ(x (1) ), . . . , δ(x (n) )) t , • r(x * ) = σ 2 δ (R θ,q (x (1) -x * ), . . . , R θ,q (x (n) -x * )), since r(x * ) = E[δ(x * )δ n ].
In the following, unknown regression and correlation parameters β, σ δ , θ and q are estimated by maximizing a likelihood function, see [START_REF] Fang | Sudjanto Design and modeling for computer experiments Chapman and Hall/CRC[END_REF]:

β = (F t R -1 θ,q F ) -1 F t R -1 θ,q Y σ2 δ = 1 n (Y -F ) t R -1 θ,q (Y -F ) ( θ, q) = arg min θ∈I R +d q∈[1,2] d σ2 δ |R θ,q | 1/n

POLLUTANTS EMISSION TESTING

The problem which we handle consists in a contribution to the modeling of the pollutants emissions, with the aim of being able to make their rate predictable, most upstream to the phase of engine development. Obviously, this problem is not new but remains critical, see [START_REF] Christakos | On the problem of permissible covariance and variogram models[END_REF]. Besides, the good choice of the modeling method and detection of outliers are important as in any domain where data exploitation is used to make prdictions, see Zhu et al. (2010).

As it has been said previously, our main objective is to establish a comparison between the two Kriging methods, especially for the industrial problem in interest. Are these two approaches equivalent? If they are not, which is the best of them regarding the specificity of our problem? To achieve this goal, we will apply those techniques over data received from a car manufacturer. The company wants to improve its engine responses owing to model based methods. In our case, the goal is to create models representing the pollutant emissions of a diesel engine. In order to be compatible with the engine tuning process, the inputs and the outputs of the models have been imposed by the cars manufacturer. The inputs of our models

x (i) = (x (i) 1 , , x (i)
d ) are the following ones:

• the opening percentage of air flux, x (i) 1

• the rail pressure, x (i) 2

• the mass fuel quantities,x (i) 3

• the instants of fuel injection, x (i) 4

• the engine speed, x (i) 5

• rhe EGR quantity, x (i) 6

• the manifold pressure, x (i) 7

The outputs of the models are the following ones:

• The NOx quantity in ppm (parts per million)

• The Particule quantity in FSN (Filter Smoke Number) • The CO quantity in ppm (parts per million)

To solve this problem, we dispose of two databases. The first one contains 317 points. To test our models, we decided to use 267 points to create the models and 50 points to validate them. Consequently, we dispose of a learning database of 267 points and a prediction/validation database of 50 points. The second data base is composed of 90 points that are divided into a learning database of 75 points and a prediction/validation database of 15 points. the second database is smaller than the first one because the engine speed parameter is not taken into account We choose these two bases because they correspond to the two functioning modes possible in the engine tuning. The adjuster builds at first a representative global model of the behaviour of engine on the entire domain which interests him. Then, he builds several models corresponding to particular and more representative functioning points.

COMPARAISON OF THE TWO KRIGING APPROACHES

At first time, we want to be sure that both types of kriging, previously presented, can give a response to the modeling of pollutants emitted by diesel engine. To do it, we began with the modeling of NOxand particles. We do not consider CO at first, because, as we shall see it afterward, this measure is rather delicate and the databases can be erroneous. Furthermore, NOx are the pollutants that are modeled the most easily because their formation depends only on the temperature of combustion. It is why they are often used to test the good accuracy of the models which we want to elaborate.

In the Geostatistic Kriging context, we used a variogram whose the formula is the following one:

y = b + c * (r) a .
Considering the first database, figure 1, 2 and 3 show the results obtained. We can see that whatever the Kriging method used, the R 2 coefficient is satisfactory and is greater than 0.85. However, in the figures 2, 3, for the NOx, the R 2 coefficient obtained by the geostatistical approach is about 0.94 and is better than the R 2 for Gaussian Kriging model (R 2 = 0.87). Moreover, we can notice that the fit of the variogram is correct with this function.

The results obtained for the particles are quite similar. We reproduce the same reasoning with the second database which contains fewer test points than the first one. In the same way, we consider a power variogram model for the geostatistical approach: For the NOx, which are not represented here, the results are quite satisfactory and completly equivalent with R 2 = 0.97 for the geostatistical kriging, and R 2 = 0.95 for the Gaussian Kriging. However, for this pollutant, the variogram is less adjusted. This can be explained by the small size of the database which we have and which limits the adjustment of the experimental variogram. As for particles, the results are just satisfactory with a value of R 2 around 0.85. This is explained partially by the insufficient size of the database but is especially due to the complexity of the physical phenomenon We can conclude that the Kriging methods offer a good response to the problem of modelling pollutants emissions such as NOx and particles. Nevertheless, it often arrives that the measures contain errors. This is due to the context in which these measures are made and to the necessity of keeping instruments of post-treatment on the engine in spite of the modifications which it leads on the conditions of the combustion.

y = b + c * (r) a . 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3

CO MODEL AND ATYPICAL POINTS

In what follows, we are interested in the pollutant CO by using the data contained in the first database. We notice that the model obtained by the geostatistical kriging does not allow a correct prediction of the CO, in spite of the fact that the variogram fits correctly the experimental data. We notice that the gaussian kriging gives a very good value of R 2 equal to 0.95 (figure 9). Indeed, the results are linked to the presence of an atypical point in the learning database. We could consider this point as an outlier. However, nothing can allow considering it as erroneous datum of a physical viewpoint. It is possible that it corresponds to a particular mode of engine functioning or to the fact that the instruments of post-treatment, need a cleaning as for example for the Anti-particles Filter (FAP) which requires regularly a regeneration operation.

To validate this supposition, we eliminate this point of our sample and we build the model again. The following results are obtained. We notice that R 2 has a highest value 0.94 (figure 11) and is comparable to that one obtained with the Gaussian Kriging. Then, two hypothesis are possible. The first one is that the point is an outlier. Thus we have to wonder why the Gaussian model did not react to its presence while the geostatistical model is very sensitive. In the second hypothesis, this value is an atypical one, revealing an important phenomenon. In that case, the capacity of the geostatistical Kriging to have a different behaviour following the presence or no of atypical data can be very useful to inform us about particular behaviour of the engine. These hypothesis are very important specially if the tests must be exploited very quickly (even in real time) to detect points of particular functioning during engine tuning and calibration. The hypothesis was tested over the second data base, still on CO modelling: The figures show that over this database it seems that there is no atypical point and in this case the two approaches produce quite the same results. The last test done is the introduction in this database of an outlier, which is illustrated by the two figures 15 and 16. In this case, the R 2 of the two models are degraded even if the statistical approach is still better than the geostatistical one. Thus, the comparison between the two approaches can provide an indication about the nature of a point: if it is an outlier or if it is an atypical point. If the statistical krigeage gives acceptable results whereas the geostatistical one not, it can be assumed that there is an atypical point that must be studied. The results we have described below, allow us to assert that in the field of diesel pollutants emissions, the kriging approaches can be mobilized. However, as outliers or atypical data are present or not, both methods have different performances. The geostatistical approach is recommended especially in our case, where often an atypical point corresponds to a particular behavior of the engine. These first results encourage us to continue the study, first to demonstrate its validity for modelling of other pollutants, secondly by studying other models of kriging, such as cokriging, for example.

Fig. 3 .

 3 Fig. 1. adjusted variogram model

Fig. 6 .

 6 Fig. 4. adjusted variogram model

Fig. 9 .

 9 Fig. 7. adjusted variogram model

Fig

  Fig. 10. adjusted variogram model

Fig

  Fig. 12. adjusted variogram model

Fig. 14 .

 14 Fig. 14. statistic kriging model -CO measures vs estimations

  Fig. 15. geostatistical kriging model -CO measures vs estimations


D Moulay

N Verdière
email: verdiern@univ-lehavre.frviden@orange.fr

L Denis-Vidal

Performance, interopérabilité et sécurité pour le développement durable Identifiablility of parameters in an epidemiologic model modeling the transmission of the Chikungunya

Keywords: Identifiability, Nonlinear models, Epidemiologic model, Chikungunya virus

In the last years, several epidemics have been reported in particular the chikungunya epidemic on the Réunion Island. For predicting its possible evolution, new models describing the transmission of the chikungunya to the human population have been proposed and studied in the literature. In such models, some parameters are not directely accessible from experiments and for estimating them, iterative algorithms can be used. However, before searching for their values, it is essential to verify the identifiability of models parameters to assess wether the set of unknown parameters can be uniquely determined from the data. Thus, identifiability is particularly important in modeling. Indeed, if the model is not identifiable, numerical procedures can fail and in that case, some supplementary data have to be added or the set of admissible data has to be reduced. Thus, this paper proposes to study the identifiability of the proposed models by (Moulay, Aziz-Alaoui & Cadivel 2011).

Introduction

The chikungunya virus is a vector-borne disease transmitted by mosquitoes of Aedes genus. Several epidemics of this tropical disease have been reported this last 50 years. Recently an unprecedented epidemic has been observed in the Réunion island (a French island in the Indian Ocean) in [2005][2006] where one third of the total population has been infected. A pic of 40 000 infected per week has been reached in february 2006. An other chikungunya epidemic has been reported in Italy in 2007. It was the first time that such disease is observed in a non tropical region. The responsible vector of these two epidemics is identified: the Aedes Albopictus mosquito [START_REF] Reiter | Aedes albopictus as an epidemic vector of chikungunya virus: another emerging problem?[END_REF]. Contrary to Aedes Aegipty, the main vector of Dengue, which also transmits the chikungunya virus, Aedes Albopictus has developed capabilities to adapt to non tropical region. Chikungunya is now a major health problem. European health authorities are now strongly engaged in the control of this disease. Since there is no vaccine nor specific treatment, efforts are mostly directed towards prevention measures and the control of mosquito proliferation. Since these events, several works and models are proposed to try to understand their emergence or re-emergence. Various fields of research are concerned, such as epidemiology, biology, medicine or mathematics. For instance, Dengue, a vector borne disease mainly transmitted by Aedes Aegipy mosquitoes was the subject of several studies [START_REF] Esteva | A model for dengue disease with variable human population[END_REF][START_REF] Esteva | Analysis of a dengue disease transmission model[END_REF]. Models for the chikungunya virus have been recently proposed [START_REF] Dumont | On a temporal model for the chikungunya disease: Modeling, theory and numerics[END_REF], [START_REF] Moulay | The chikungunya disease: Modeling, vector and transmission global dynamics[END_REF].... Since the models are recent, the not well-known parameters have not yet been studied. In this paper, we propose to take again the models proposed by [START_REF] Moulay | The chikungunya disease: Modeling, vector and transmission global dynamics[END_REF] and to do an identifiability study. In their paper, the models are uncontrolled and can be described in a general state-space form:

Γ θ = ẋ(t, θ) = f (x(t, θ)), y(t, θ) = h(x(t, θ), θ). (1)
Here x(t, θ) ∈ R n and y(t, θ) ∈ R m denote the state variables and the measured outputs, respectively and θ ∈ U p the unknown parameters vector (U p is an open subset in R p ). The functions f (., θ) and h(., θ) are real, rational and analytic for every θ ∈ U p on M (a connected open subset of R n such that x(t, θ) ∈ M for every θ ∈ U p and every t ∈ [0, T ]).

The identifiability definition of the uncontrolled model Γ θ is the following: The identifiability of models has been extensively studied [START_REF] Ljung | On global identifiability for arbitrary model parametrizations[END_REF], [START_REF] Vajda | Similarity transformation approach to structural identifiability of nonlinear models[END_REF], [START_REF] Verdière | Identifiability and estimation of pharmacokinetic parameters for the ligands of the macrophage mannose receptor[END_REF] and different approaches have been proposed for studying the global identifiability of nonlinear systems. We can mention for example, the Taylor Series approach of [START_REF] Pohjanpalo | System identifiability based on the power series expansion of the solution[END_REF].

He proposed a method based on the analysis of a power series expansion of the output which gives rise to an algebraic system constituted of an infinite number of equations. A second method is based on the local state isomorphism theorem ( [START_REF] Walter | Global approaches to identifiability testing for linear and nonlinear state space models[END_REF], [START_REF] Chappell | Structural identifiability of the parameters of a nonlinear batch reactor model[END_REF], (Denis-Vidal, Joly-Blanchard & Noiret 2001), [START_REF] Chapman | Structural identifiability of nonlinear systems using linear/nonlinear splitting[END_REF]). It leads to study the solution of a specific set of differential partial equations. A third one is a method based on differential algebra that was introduced by (M. Fliess 1993), [START_REF] Ljung | On global identifiability for arbitrary model parametrizations[END_REF] and [START_REF] Ollivier | Identifiabilité des systèmes[END_REF]). It allows one to obtain relations linking the observations, the inputs and the unknown parameters of the system. These relations can be used to obtain a first estimation of the unknown parameters without a priori any knowledge of them [START_REF] Verdière | Identifiability and estimation of pharmacokinetic parameters for the ligands of the macrophage mannose receptor[END_REF]. It is the latter method which will be used in this paper for studying the models identifiability.

The paper is organized as follows. In the second section, models describing the transmission of the chikungunya virus to human population are presented. Some results obtained in (Moulay, Aziz-Alaoui & Cadivel 2011) will be recalled since they will give us first, the framework of our study then, the steps to study the identifiability of the not wellknown parameters. In the third section, the identifiability results are given.

Presentation of the models

In [START_REF] Bacaër | Approximation of the basic reproduction number r0 for vector-borne diseases with a periodic vector population[END_REF]) the author formulate several methods to compute the basic reproduction number for epidemiological models. One of the first models describing the chikungunya transmission virus using SI-SIR type models is proposed. Moreover, some biological parameter values are given. An other approach is describe in [START_REF] Dumont | On a temporal model for the chikungunya disease: Modeling, theory and numerics[END_REF], where a global aquatic stage for the mosquito dynamics supplements a classical transmission model. In [START_REF] Dumont | Vector control for the chikungunya disease[END_REF], authors formulate an ordinary differential equation system to study control of chikunugunya virus using mechanical and chemical tools. In [START_REF] Moulay | Optimal control of chikungunya disease: Larvae reduction, treatment and prevention[END_REF], control efforts are taken into account through the formulation of an optimal control problem, where the objective is to control the mosquito proliferation and limit the number of human and mosquito infection. This papers deal with the Réunion Island epidemic.

Our model given in [START_REF] Moulay | The chikungunya disease: Modeling, vector and transmission global dynamics[END_REF][START_REF] Moulay | Optimal control of chikungunya disease: Larvae reduction, treatment and prevention[END_REF] takes into account the mosquito biological life cycle and describe the transmission virus to human population. For the reader convenience, we briefly recall the modeling steps. The mosquito biological life cycle consists in four stages: eggs, larvae, pupae and adults. We use a stage structured model to describe the following stages: eggs (E), larvae and pupae (L, two stages biologically close) and female adults (A, only females can transmit the virus) stages. The density variation of each stage is describe by the following scheme:

density variation = entering -(leaving + death)
The egg density variation is then described by the number of eggs laid by females b, by eggs becoming larvae with a transfer rate s and by eggs death with a natural mortality rate d. We assume that the number of eggs is proportional to the number of females b(t)A(t), and regulated by a carrying capacity K E since mosquitoes are able to detect the best breeding site ensuring the egg development, then b

(t) = bA(t)(1 -E(t)/K E ).
Other stages, are described in the same way. The input in the larvae stage, given with a transfer s is also assumed to be regulated by a carrying capacity K L which characterizes the availability of nutrients and space. The number of new larvae entering the L stage is then given by s(t) = sE(t)(1 -L(t)/K L ). These larvae become adult females with a transfer rate s L . Natural deaths occur with a rate d L , d m for larvae and adults respectively.

This model is then included in a classical SI-SIR epidemiological model to describe the virus transmission to human population. To this aim, the adult stage A is divided into two epidemiological states: susceptible S m and infective I m , since mosquitoes carry the infection along their life. All previous assumptions are summed up in Fig. 1. Based on our model description (see Fig. 1) and assumptions, we establish the following equations:

Adult

                                                         dE dt (t) = bA(t) 1 E(t) K E -(s + d)E(t) dL dt (t) = sE(t) 1 - L(t) K L -(s L + d L )L(t) dA dt (t) = s L L(t) -d m A(t) d Sm dt (t) = s L L(t) -d m Sm (t) -β m ĪH (t) N H (t) Sm (t) d Īm dt (t) = β m ĪH (t) N H (t) Sm (t) -d m Īm (t) d SH dt (t) = -β H Īm (t) A(t) SH (t) -d H SH (t) +b H ( SH (t) + ĪH (t) + RH (t)) d ĪH dt (t) = β H Īm (t) A(t) SH (t) -γ ĪH (t) -d H ĪH (t) d RH dt (t) = γ ĪH (t) -d H RH (t) (2
                                        E (t) = bA(t) 1 - E(t) K E -(s + d)E(t) L (t) = sE(t) 1 - L(t) K L -(s L + d L )L(t) A (t) = s L L(t) -d m A(t) (a)        S H (t) = -(b H + β H I m (t)) S H (t) + b H I H (t) = β H I m (t)S H (t) -(γ + b H )I H (t) I m (t) = -s L L(t) A(t) + β m I H (t) I m (t) + β m I H (t) (b) (3)
and it is defined on ∆ × Ω where

∆ =      (E, L, A) | 0 ≤ E ≤ K E 0 ≤ L ≤ K L 0 ≤ A ≤ s L d m K L      (4)
and

Ω = (S H , I H , I m ) ∈ (R + ) 3 | 0 ≤ S H + I H ≤ 1 0 ≤ I m ≤ 1 .
(5)

The stability analysis of the model is detailed in [START_REF] Moulay | The chikungunya disease: Modeling, vector and transmission global dynamics[END_REF]. We briefly recall some results about this model. The study was conducted in two steps and they will be taken again for the identifiability study. First, we analyze the mosquito dynamics in the absence of virus, which correspond to the subsystem (3a). The mosquito dynamic is governed by the following threshold:

r = b s + d s s L + d L s L d m (6)
obtained from computation of the equilibrium.

Theorem 2.1.

• System (3a) always has the mosquito-free equilibrium (0, 0, 0), which is globally asymptotically stable (GAS) if r ≤ 1 and unstable otherwise

• If r > 1 system (3a) has an endemic equilibrium (E * , L * , A * ) wich is GAS, where

  E * L * A *   = 1 - 1 r        K E γ E K L γ L s L d m K L γ L        γ E = 1 + (s+d)dmK E bs L K L and γ L = 1 + (s L +d L )K L sK E
MOSIM'12 -06 au 08 juin 2012 -Bordeaux -France

In both cases, the global stability is obtained by using Lyapunov function theory. Now we assume r > 1, the biological interesting case, in order to ensure the persistence of mosquito population and we consider the subsystem (3b).

The stability of equilibrium of the transmission dynamics model is described thanks to the basic reproduction number (van den [START_REF] Van Den Driessche | Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission[END_REF][START_REF] Diekmann | Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis and Interpretation[END_REF], computed in the case r > 1 which is the biologically interesting case:

R 0 = β m β H d m (γ + b H ) (7)
We show the following result Theorem 2.2. Assume r > 1 and let us denote (E * , L * , A * ) the endemic equilibrium of (3a).

• System (3b) always has the disease-free equilibrium (1, 0, 0), which is GAS if R 0 ≤ 1 and unstable otherwise.

• If R 0 > 1 system (3a) has an endemic equilibrium (S * H , I * H , S * m ) which is GAS and where

  S * H I * H I * m   =        b H β H + b H + β H (β H + b H )R 0 d m b H β m (β H + b H ) (R 0 -1) b H β H + b H R 0 (R 0 -1)       
The first part of the theorem is obtained using Lypunov function theory. The case of the endemic equilibrium needs more study. The idea here is that the mosquito dynamic system drove the transmission dynamics. It may be assimilated to master-slave system.

The coupling term is s L L(t) A(t) .

In order to study the equilibrium stability we consider the limit system associated to (3b) and use the result of [START_REF] Vidyasagar | Decomposition techniques for large-scale systems with nonadditive interactions stability and stabilizability[END_REF]) :

Theorem 2.3. Consider the following

C 1 system      dx dt = f (x) dy dt = g(x, y), (8)
with (x, y) ∈ R n × R m . Let (x * , y * ) be an equilibrium point. If x * is GAS in R n for the system dx dt = f (x) and if y * is GAS R n for the system dy dt = g(x * , y)

, then (x * , y * ) is (locally) asymptotically stable for system (8). Moreover, if all trajectories of (8) are forward bounded, then (x * , y * ) is GAS for (8).

The GAS of the endemic equilibrium (S * H , I * H , S * m ) of system (3b) where s L L(t) A(t) is replaced by s L L * A * is then shown using the theory of competitive systems [START_REF] Hirsch | Differential equations, dynamical systems, and linear algebra[END_REF], [START_REF] Hirsch | System of differential equations that are competitive or cooperative. iv: structural stability in three-dimensional systems[END_REF]) [START_REF] Smith | Monotone Dynamical Systems : An introduction to the theory of competitive and cooperative systems, American Mathematical Society[END_REF] and the Poincaré-Bendixson property [START_REF] Thieme | Convergence results and a poincare -bendixson trichotomy for asymptotically autonomous differential equations[END_REF].

Identifiability Analysis

Recall that the identifiability analysis of models parameters consists in assessing wether the set of unknown parameters can be uniquely determined from the data. Thus, it is essential to determine the state variables that can be considered as observable. In the case of the chikungunya Réunion Island epidemic, authorities have registered the average number of eggs in a number if sites. Thus, (E) can be considered as an observable variable. Furthermore, they estimate the number of new infection week by week. More generally, it seems to be realistic to assume that data about human population may be obtained. For instance, we know that the entire Réunion island before the epidemic was susceptible. Data indicating week per week new cases of the disease may be provided by the INVS (French Institute for Health Care). We know that the epidemic was declared over by April 2006. In the end the INVS counted 265,733 cases of chikungunya from March 2005 to April 2006. This represents more than 35% of the total population of the Island. That is why it seems reasonable to assume that susceptible (S H ) and infected human (I H ) are observable. The parameters whose values are not directly accessible from the fiel are: s, s l , K e , K l for the system (3a) and d l , d m for the system (3b). Let us recall main results in differential algebra for proving the parameters identifiability.

Differential Algebra

This method consists in eliminating unobservable state variables in order to get relations between outputs and parameters. Let us recall the methodology. The system Γ θ is rewritten as a differential polynomial system completed with θi = 0, i = 1, . . . , p, thus the following system composed of polynomial equations and inequalities is obtained:

       p( ẋ, x, θ) = 0, q(x, y, θ) = 0, r(x, y, θ) = 0, θi = 0, i = 1, . . . , p. (9)
Let us introduce some notations:

• I is the radical of the differential ideal generated by ( 9). I, endowed with the following ranking which eliminates the states variables:

[θ] ≺ [y, u] ≺ [x] (10)
MOSIM'12 -06 au 08 juin 2012 -Bordeaux -France is assumed to admit a characteristic presentation C (i.e., a canonical representant of the ideal) which has the following form: θ1 , . . . θp , P 1 (y, u, θ), . . . , P m (y, u, θ), Q 1 (y, u, θ, x), . . . , Q n (y, u, θ, x), }

• I θ is the radical of the differential ideal generated by ( 9) for the particular value of parameter θ and C θ is the characteristic presentation associated with the ranking [y, u] ≺ [x].

• Finally, I i0 θ is the ideal obtained after eliminating state variables, the set

C i0 θ = C θ ∩ Q(θ){U, Y } is a characteristic presentation of this ideal.
The authors in (Denis-Vidal, Joly-Blanchard, Noiret & Petitot 2001) have given some technical conditions for having the equality C θ = C(θ). Under these assumptions, the characteristic presentation C θ , that is, C i0 θ of I i0 θ is proved to contain the differential polynomials P 1 (y, u, θ), . . . , P m (y, u, θ) which can be expressed as

P i (y, u, θ) = γ i 0 (y, u) + ni k=1 γ i k (θ)m k,i (y, u) (12)
where (γ i k ) 1≤k≤l are rational in θ, γ i u = γ i v (u = v), (m k,i ) 1≤k≤l are differential polynomials with respect to y and u and γ i 0 = 0.

The list {γ i 1 (θ), . . . , γ i ni (θ)} is called the exhaustive summary of P i . The size of the system is the number of observations. The identifiability analysis is based on the following proposition (Denis-Vidal, Joly-Blanchard, Noiret & Petitot 2001).

Proposition 3.1. If for i = 1, . . . , m, P i (y, u, θ) = det(m k,i (y, u), k = 1, . . . , n i ) is not in the ideal I i0 θ , then Γ θ is globally identifiable at θ if and only if for every θ ∈ U θ ( θ = p), the characteristic presentations C i0 θ and C i0 θ are distinct.

For studying the identifiability of the parameters s, s l , K e , K l , d l , d m in (3a) and (3b), the two coupled systems can be considered as a unique system in which E, S H and I H are supposed to be observed. However, we will take again the procedure done in [START_REF] Moulay | Optimal control of chikungunya disease: Larvae reduction, treatment and prevention[END_REF] and presented in section 2, that is, decompose the identifiability analysis in two steps. Indeed, for studying the parameters of the second system it is essential to know those of the first one. Besides, our aim is to propose an identifiability study which can be used for a numerical procedure. Indeed, as it was done in [START_REF] Verdière | Identifiability and estimation of pharmacokinetic parameters for the ligands of the macrophage mannose receptor[END_REF], the use of differential algebra gives output polynomials usable for estimating the unknown parameters.

Application to the Vector population

Since E is supposed to be observed, the equation y = E is added to the sytem (3a). In using the elimination order [y] ≺ [E, L, A], the package diffalg of Maple gives the caracteristic presentation constituted of the three following polynomials (13):

P 1 = (-bK e + by)A + ẏK e + K e sy + K e dy P 2 = (bK 2 e s l -2K e bys l + by

2 s l )L -K 2 e d ẏ -K 2 e s ẏ -K 2 e ÿ + K e ÿy -K 2 e d m ẏ -K 2 e d m sy -K 2 e d m dy -K e ẏ2 + d m K e ẏy + d m K e sy 2 + d m K e dy 2 P 3 = (K 3 e s l K l d m d + K 3 e K l d l d m s + K 3 e K l d l d m d +K 3 e s l K l d m s -bK 3 e s l sK l )y + (K 3 e K l d l s + K 3 e K l d m d +K 3 e K l d l d m + K 3 e s l K l d + K 3 e s l K l s + K 3 e s l K l d m +K 3 e K l d l d + K 3 e K l d m s) ẏ + (K 3 e K l d l + K 3 e K l d +K 3 e K l d m + K 3 e K l s + K 3 e s l K l )ÿ +K l K 3 e ... y + (3K 2 e bs l sK l -2K 2 e s l K l d m d + K 3 e sd m d -2K 2 e K l d l d m s -2K 2 e K l d l d m d + K 3 e s 2 d m -2K 2 e s l K l d m s)y 2 +(-2K 2 e K l d l d m + K 3 e sd m -K 2 e K l d l s -K 2 e K l d l d + K 3 e sd -K l d m K 2 e d -K 2 e s l K l s -K 2 e s l K l d -2K 2 e s l K l d m -K l d m K 2 e s +K 3 e s 2 ) ẏy + (K 2 e K l d l + K l d m K 2 e + 2K 2 e K l s + 2K 2 e K l d+ K 2 e s l K l ) ẏ2 + (-2K 2 e K l d l -2K l d m K 2 e -K 2 e K l s -2K 2 e s l K l -K 2 e K l d + K 3 e s)ÿy + 3K l K 2 e ÿ ẏ -2K l K 2 e ... y y + (s l K l d m K e s -2K 2 e s 2 d m -3K e bs l sK l +K l d l d m K e s + s l K l d m K e d -2K 2 e sd m d + K l d l d m K e d)y 3 +(-K 2 e s 2 -2K 2 e sd m -K 2 e sd + s l K l d m K e + K l d l d m K e ) ẏ y 2 + (-K e K l d l + K 2 e s -K e s l K l -K l d m K e ) ẏ2 y +2K l K e ẏ3 + (-2K 2 e s + K l d m K e + K e K l d l + K e s l K l )ÿy 2 -3K l K e ÿ ẏy + K l K e
... y y 2 + (s 2 d m K e + sd m K e d +bs l sK l )y 4 + sd m K e ẏy 3 -sK e ẏ2 y 2 + sK e ÿy 3 .

(13)

The polynomials P 1 and P 2 used to express L and A as a function of y and the parameters of the model. The third one, P 3 , links the output with the parameters: it is the output polynomial. With the function belong to, we verifie that the functional determinant (P 3 ) is not in the ideal I i0 θ . The exhaustive summary is constituted of 21 expressions. In using the Rosenfeld-Groebner algorithm, we obtain the identifiability of the parameters K e , K l , s, s l . Thus, from the observation E, the unknown parameters can be estimated.

Application to the population Model

The third equation of (3b) links the human population to the vector population with the term L(t)/A(t). According to the previous section, they can be explicitely determined from E(t) thus the ratio L(t)/A(t) can be considered as a known input u. As previously, in adding y 1 = I H , y 2 = S H to (3b) and in considering the elimination order [y 1 , y 2 , u] ≺ [I h , S H , I m ], one gets for the two following output (14)

Only the polynomial P 4 contains the parameters β H and β m and is used for studying their identifiability. The functional determinant P 4 is proved not to be in the ideal I i0 θ . In studying the exhaustive summary of P 4 , we conclude that the parameters β H and β m are identifable.

Conclusion

In this paper, the identifiability of models describing the transmission of the chikungunya virus to human population has been studied. According to the reemergence of this virus, the chikungunya becomes a major health problem especially since the main vector has developed capabilities to adapt to non tropical region. The identifiability is an important step in the modeling. Indeed, the identifiability study enables one to know if a model is well-posed and if the unknown parameters can be assess from some observations done in the field.

Fig. 1 :

 1 Fig. 1: Compartmental model for the dynamics of Aedes albopictus mosquitos and the virus transmission to human population.

)

  Using the following variable changes S m = Sm /A, I m = Īm /A S H = SH /N H , I H = ĪH /N H and R H = RH /N H and the fact that then S m = 1 -I m et R H = 1 -S H -I H , system (2) reads as:

  MOSIM'12 -06 au 08 juin 2012 -Bordeaux -France polynomials:P 4 = y 2 ÿ2 -ẏ2 2 + (β H β m + β m b H )y 1 y 2 2 + s l uy 2 ẏ2 +s l b H uy 2 2 -s l b H uy 2 + b H ẏ2 + β m y 2 y 1 ẏ2 -β m b H y 2 y 1 P 5 = ẏ2 + ẏ1 -b H + b H y 2 + (b H + γ)y 1 .