Rémi Colom

Lei Xu

Loris Marini

Frédéric Bedu

Igor Ozerov

Thomas Begou

Julien Lumeau

Andrey Miroshnishenko

Dragomir Neshev

Boris Kuhlmey

Stefano Palomba
email: stefano.palomba@sydney.edu.au

Nicolas Bonod
email: nicolas.bonod@fresnel.fr

Enhanced Four-Wave Mixing in Doubly Resonant Si Nanoresonators

Frequency conversion is one of the main applications of nonlinear optical processes in which a signal is produced at a different wavelength from the excitation wavelength.

In particular, four-wave mixing (FWM) is a third order non-linear optical process that allows, for instance, the generation of visible frequencies by tuning near-infrared laser pumps. Here, in order to augment the very weak FWM conversion efficiency, we design silicon Mie resonators that exhibit two resonances of the internal electric field intensity 1 around the frequency range of the laser pumps. The linear extinction spectrum of the individual Si resonator is first measured by bright field spectroscopy and compared with numerical simulations to confirm the existence of the two resonances corresponding to electric and magnetic dipoles excitations. The FWM signal is then measured for a single Si nanoresonator when the first pump is set to the electric resonance, while tuning the frequency of the second pump across the magnetic dipolar resonance. We show that the FWM signal generated in the visible spectrum is maximum when the frequency of the tunable pump corresponds to the maximum of the internal electric field intensity.

At this position, the FWM signal is enhanced by more than two orders of magnitude compared with the FWM signal generated by the unpatterned silicon film.

Introduction

Silicon is an ubiquitous material in opto-electronics with applications in image sensors or photovoltaics. However, the spectral range of applications of silicon is limited in the visible spectrum since its band gap is located near 1.1 eV for crystalline Si. As a consequence, silicon is transparent in the near infrared, and behaves as a dielectric characterized by a high refractive index ≈ 3.4. Making silicon sensitive to near-infrared (NIR) radiation is very challenging and multiphoton absorption is a promising way to address this challenge. Four wave mixing (FWM) is a third order nonlinear effect where four photons of different frequencies are mixed together and, for instance, light in the NIR can generate visible (VIS) radiation. [1][2][3][4][5] Such nonlinear effects are generally weak since the nonlinear susceptibility χ (3) is small. However, given that the nonlinear field is proportional to the product between the nonlinear susceptibility of the nonlinear medium and the applied electric field, the nonlinear field can be enhanced near the resonances. 6 For this purpose, optical cavities can be designed to trap light inside the nonlinear material. In particular, high refractive index subwavelength-sized particles can host low order Mie resonances that are associated with strong electric field intensities inside or in the vicinity of the particle. [7][8][9][10][11][12][13][14][15][16] Silicon nanoparticles behave as optical antennas able to efficiently couple the far field radiation to the near field, with applications in fluorescence or in biosensing. [17][18][19][20][21][22][23][24] The possibility to engineer the near-field intensity via optical antennae has opened a novel route to develop nonlinear optics with nanostructured semi-conductors. 25,26 The first demonstrations were with Mie resonators designed in group IV semiconductors (silicon, germanium) yielding third harmonic generation (THG). [27][28][29][30][31][32] Alldielectric nanostructures made of crystalline silicon 31 or III-V semiconductors [33][34][35][36] can also be engineered to enhance the process of second harmonic generation (SHG) due to their non-centrosymmetric crystalline structure. III-V semiconductors feature excellent nonlinear properties and the coupling between their intrinsic nonlinearity with internal field enhancements in Mie resonators has strong potential for the development of highly efficient nonlinear all-dielectric metasurfaces. Besides SHG and THG, four-wave mixing (FWM) is a versatile nonlinear process that allows to tune the frequency of the nonlinear signal by modifying the frequency of one pump. 37,38 In the context of Mie resonators, FWM allows for spectrally tuning the pump around the frequency of a Mie resonance in order to maximize the nonlinearity enhancement by compensating for material inhomogeneities and fabrication inaccuracies.

Here we generate visible photons by degenerate FWM when two NIR pumps are exciting an individual silicon nanodisk (Fig. 1a,b). The first pump is delivered by a tunable Ti-Sapphire laser (Coherent Chameleon Ultra) generating 140 fs pulses with a central wavelength ranging between 750 nm and 880 nm. The second pump is obtained by an Optical Parametric Oscillator (OPO) (Coherent Chameleon) coupled to the main oscillator (Chameleon Ultra), as shown in Fig. 1c. The mean wavelength of the second pump can also be tuned between 1050 nm and 1500 nm. Since the two femtosecond laser pulses have to overlap both in space and in time, the arrival time of the first beam with frequency ω 1 is tuned by an optical delay line whereas the spatial overlap is controlled by positioning the point spread functions of the two pumps at the image plane of a home-built confocal microscope by a 100× and high numerical aperture (NA= 0.85) air microscope objective (Olympus LCPLN100XIR). The visible nonlinearly generated photons from individual resonators are collected by the same objective and analyzed by an optical spectrometer (Princeton Instruments, Isoplane 320) coupled to a visible CCD camera (Princeton Instruments, Pixis 1024) (Fig. 1a). A schematic representation of the FWM process used is depicted in Fig. 1b. Degenerate FWM is a third order nonlinear process in which a signal at the frequency ω s = 2ω 1 -ω 2 is generated when the sample is illuminated by the two pumps at frequencies ω 1 and ω 2 .

Assuming that the material we consider is amorphous (this assumption shall be discussed hereafter), the non-linear polarization associated with this process has the following form:

P N L,ωs (r) = 6 0 χ (3) 1122 E ω 1 (r) • E ω 2 (r) E ω 1 (r) + 3 0 χ (3) 1221 (E ω 1 (r) • E ω 1 (r)) E ω 2 (r). (1)
where χ

Results and discussion

Silicon nanodisks are designed to exhibit two resonances yielding internal field maxima near the wavelengths of the two chosen pumps ω 1 ≈ 850 nm and ω 2 ≈ 1200 nm. The nanodisks are fabricated by nanopatterning of a thin film of amorphous silicon. However, illumination of amorphous silicon by femto-second laser pulses leads to a crystallization of silicon as will be discussed later. 39,40 The actual FWM measurements will then be performed after the crystallization of silicon. Subsequently, in our simulations we use the dispersion data of crystalline silicon for calculating the optical response of the silicon nanodisks on a glass substrate.

Let us consider the average of the internal field enhancement E int /E 0 over the whole volume of the silicon disk:

|E int | = V D |E int |dr V D |E 0 |dr , (2)
where V D is the volume of the scatterer while E 0 is the amplitude of the excitation field. The internal field enhancement and its average in the particle |E int | are calculated by Finite Element Methods (FEM) (COMSOL Multiphysics). We found that a silicon cylinder with a diameter D = 340 nm and a height h = 240 nm on a glass substrate, exhibits two resonances of the internal field |E int | at ω 1 and ω 2 , for which the spectrum is plotted in Fig. 2a.

The two maxima in the spectral field distribution are located in the NIR. The nature of these two resonances can be determined by computing a multipolar decomposition of the optical response of the nanodisk. In Fig. 2c, the extinction cross-section of the nanodisk (D = 340 nm and a height h = 240 nm) and its multipolar expansion are plotted as a function of the wavelength.

The spectral peak at λ ≈ 1205 nm is associated with a magnetic-dipolar resonance.

This multipolar analysis is confirmed by observing that the calculated field distribution at λ = 1205 nm in Fig. 2d is similar to a magnetic dipole field distribution. A counterpart to the second peak of the internal field average can also be found in the extinction spectrum located at λ ≈ 840 nm. The multipolar analysis of this second peak reveals that it is mostly due to an electric dipolar response combined with a weaker magnetic quadrupolar contribution.

Two peaks of the ED can be observed in the spectrum near λ ≈ 840 nm and λ ≈ 1050 nm that are associated with very different electric field distributions inside the nanoresonator.

Different from the resonant peak at λ ≈ 1050 nm where the electric field is mostly enhanced around the surface of the nano-resonators, the electric field for the peak at λ ≈ 840 nm, penetrates inside the nano-resonator, resulting in a strong internal field enhancement which can further facilitate light-matter interactions and augment the nonlinear process.

A careful analysis of Eq. 1 shows that the nonlinear polarization depends on two terms

E ω 1 (r) • E ω 2 (r)E ω 1 (r) and (E ω 1 (r) • E ω 1 (r)) E ω 2 (r).
We therefore plot in Fig. 2b the overlap integrals:

I N L 1 = V D E ω 1 (r) • E ω 2 (r) E ω 1 (r)d 3 r , (3)
I N L 2 = V D (E ω 1 (r) • E ω 1 (r)) E ω 2 (r)d 3 r . ( 4
)
The color maps in Fig. 2b provide information on the FWM generation with respect to the ω 1 and ω 2 frequencies showing four areas in which the FWM signal is enhanced. The peak that occurs in the area corresponding to the spectral range of the two pumps, i.e.

λ 1 = [800 -900] nm and λ 2 = [1150 -1300] nm, is more pronounced I N L 2 . As the FWM signal occurs at the wavelength λ s = λ 1 λ 2 2λ 2 -λ 1 , the peak observed on the top left of the I N L1 and I N L2 color maps, i.e. for λ 1 ≈ 840 nm and λ 2 ≈ [1200 -1300] nm, will yield a FWM signal in the visible at λ s ≈ [620 -646] nm for λ 1 ≈ 840 nm. The femtosecond pumps induce a crystallization of silicon after illumination of the nanodisk. The FWM signal in monocrystalline silicon should be rigorously modeled by three terms 42 while Eq. (1) that includes two terms is valid in amorphous silicon. However, the femtosecond pumps induce a polycrystallization of the silicon nanodisk with a random orientation of the crystalline meshes. The non-linear tensor used when modeling the FWM signal in our experiment can Before investigating the nonlinear optical properties of the fabricated silicon Mie resonators, we verified that the geometrical and optical properties of the samples are in agreement with the expected morphology and optical behaviour in the linear regime. The former is first investigated by a scanning electron beam microscope and the diameter of the disk is found in reasonable agreement with the targeted value D = 340 nm (see Fig. 3a). Second, we verified that the linear response of isolated nanodisks before the crystallization of silicon, i.e. before illuminating the silicon nanodisk with the femtosecond pumps, matches the numerical prediction of the extinction cross-section calculated with amorphous silicon (See SI).

Since the generation of a FWM signal requires an excitation of the sample by two highpower femtosecond laser pumps, a modification of the extinction spectrum is observed after the first illumination by the pumps. This change corresponds to a blue-shift of the resonance wavelengths in the extinction spectrum. This blue shift of the spectra is assumed to result from the crystallization of the initially amorphous silicon nanodisks due to the pulsed laser illumination 39,43 and is confirmed by numerical calculations of the extinction spectra performed with the refractive index of crystalline silicon (See SI). We verified that after the crystallization induced by the first illumination from the two laser pumps, the resonances do not exhibit any further spectral shift. Therefore, we show in Fig. 3b,c the stable spectra obtained after several illuminations (spectra acquired in the visible and NIR are obtained with two different cameras which explains why they are plotted in two different figures). The noise observed in Fig. 3b,c comes from the fact that extinction spectra are measured at the single resonator scale. Hence, the FWM measurements are performed only after illuminating the nanoresonators and verifying that crystallization has occurred which stabilizes the spectral position of the resonances.

The intensity of the FWM signal is measured with the camera while sweeping the wavelength λ 2 of the second pump across the NIR resonance of the nanoresonator after crystallization has occurred. The first pump is set at wavelength λ 1 = 810 nm, matching the high frequency electric dipolar resonance of the Si disk.

In order to assure that the pump excitation has not further modified the properties of the silicon resonator, the extinction of the silicon nanodisk is measured before and after each excitation of the FWM signal. The FWM signal enhancement yielded by the individual silicon resonator is compared with the FWM signal measured from an unpatterned silicon thin film. This measurement of the FWM signal on the pristine silicon film is used as a reference. In addition to this reference, all measurements are normalized by the intensity of the two pumps (normalization of the FWM signal is detailed in Methods). The FWM spectra obtained from the pristine thin film and from the silicon Mie resonator are displayed in Fig. 4 for six different wavelengths λ 2 of the pump 2. The FMW signal generated by individual nanoresonators is more than two-orders of magnitude higher than the one measured from the unpatterned Si film (Fig. 4a,b).

The spectral dependence of the FWM signal generated by the silicon film may be due to several factors. First, the third-order susceptibility is certainly dispersive. This makes the nonlinear response of the thin film dependent both on the wavelengths of the pumps and on the wavelength of the signal. Second, the dependence of the FWM signal may also be caused by the excitation of Fabry-Perot resonances inside the silicon slab. Regarding the FWM response of the silicon nanodisk, it has a maximum at λ 2 = 1200 nm and decreases for shorter and larger values of λ 2 . The significance of the Mie resonances for enhancing the FWM signal is clear when observing the two spectra in Fig. 4a,b since the FWM signal generated by the silicon particles can be two-orders of magnitude higher than that measured on the unpatterned Si film.

In order to better understand the enhancement of the FWM intensity, we plot in Fig. 5 the FWM signal intensity as a function of the wavelength λ 2 of the second pump. This signal intensity is compared to the square of the internal field enhancement defined in Eq. ( 2) calculated with the dielectric permittivity of crystalline silicon in Fig. 5. [START_REF] Green | Solar Energy Materials and Solar Cells[END_REF] A good agreement is found between the internal field enhancement and the FWM signal. This is in agreement with the fact that a FWM signal is proportional to film.

∝ |(E ω 1 (r) • E ω 1 (r)) E ω 2 (r)| 2 . This
In Fig. 5(b), we plot the FWM signal and the extinction measurement on the same figure in order to compare the peak of the FWM signal to the resonance of the nanodisk extinction.

There appears indeed to be a red-shift of the FWM peak with respect to the resonance of the extinction of the crystallized nanodisk. This red-shift is expected since the resonance of the internal field is classically red-shifted compared to the resonance of the extinction cross section. 14 The stronger extinction cross-section observed in Fig. 5(b) at shortest wavelengths, around 1100 nm, may be due to a larger scattering of light caused by the resonant excitation of an electric dipole inside the Si particle. This electric dipole is associated with a strong enhancement of light intensity around the surface of the particle, as displayed in Fig. 2d, which is expected to increase the sensitivity to surface scattering.

Conclusion

In conclusion, we combined linear and non-linear spectroscopy to report the generation of visible light by exciting a single silicon particle in the NIR spectrum with two pumps.

Linear characterization

For the linear characterization of the isolated silicon nanodisks, the extinction spectrum is measured with a confocal microscope. The sample is illuminated by a collimated white source slightly focused in order to combine intensity and homogeneity of the beam on the sample.

The light from the sample is subsequently collected by a 100× magnifying objective with a numerical aperture equal to 0.85 (Olympus LCPLN100XIR) and then focused on the slit of the photospectrometer (Princeton Instrument, Isoplane 320) that acts as a spatial filter.

In fact, only the light coming from the image of the slit onto the sample is collected by the photospectrometer. Light scattered by an isolated Si nanodisk is then collected by adequately closing the slit to spatially filter the light signal in one direction and by using the binning option of the camera to spatially filter in the other direction. Extinction measurements are then performed by measuring the light spectrum transmitted through (i) an isolated disk I d and (ii) the substrate I s and the background noise I b . Extinction spectrum is defined by

Ext(λ) = -I d (λ)-Is(λ) Is(λ)-I b (λ) .

Nonlinear characterization

Four-wave mixing measurements are performed using the experimental setup displayed in The reflected FWM signal is then collected by the same objective and subsequently sent to a spectrometer.

We measure the FWM signal from the Si nanoresonators and comparing it to the FWM signal obtained from a Si thin-film. The row data from the spectrometer are normalized in order to take into account the change of average intensity on the sample when the wavelength of the pumps is changed. For that purpose, we measured the average powers of the two pumps and normalized them by the surface of the focused spot. This area was estimated by using the radius of a diffraction limited system: r(λ) = 0.61 λ N A where N A is the numerical aperture of the objective that in our case is equal to N A = 0.85. The normalization factor for the FWM measurements is consequently equal to N norm = P λ 1 πr(λ 1 ) 2 2 P λ 2 πr(λ 2 ) 2 . Finally the signal from the silicon nanodisks and the FWM signal from the layer are both normalized by the maximum FWM intensity from the silicon layer: 3.5 * 10 22 .

Numerical calculation method

For the numerical simulation of linear scattering, we use a finite-element-method solver in the frequency domain from COMSOL Multiphysics to model the linear response from our Si resonator. As the presence of substrate makes it difficult to perform the multipolar analysis using the scattered fields directly, the multipolar decomposition analysis here is performed using the polarization currents induced inside the Si resonator. 44 By considering the fields radiated into air, we calculate and predict the contributions associated with each multipole from our Si resonator.

Figure 1 :

 1 Figure 1: (a) Schematic representation of the degenerate four-wave mixing (FWM) experiment: two laser pumps of frequency ω 1 and ω 2 are mixed in a silicon nanodisk to generate a FWM signal at frequency ω s = 2ω 1 -ω 2 . (b) FWM principle. (c) Schematic representation of the experimental set-up: the laser pump beam ω 1 is in red and ω 2 is in brown. The FWM beam is in yellow. D. M.: Dichroic Mirror.

  elements of the third-order susceptibility tensor of amorphous silicon and the nonlinear signal has a frequency ω s = 2ω 1 -ω 2 .

Figure 2 :

 2 Figure 2: (a) Spectrum in the near-infrared of the simulated internal field enhancement averaged in the whole volume of the crystalline silicon disk, 41 |E int | , chosen to feature two resonances at λ 1 and λ 2 : diameter D = 340 nm and height h = 240 nm. (b) Overlap integrals I N L1 (left) and I N L2 (right) calculated with Eqs. 3-4 and plotted with respect to λ 1 and λ 2 . (c) Scattering spectrum and multipolar decomposition calculated for the chosen Si particle (D = 340 nm, h = 240 nm) on a substrate of refractive index 1.5. Full calculation: black curve; colored curves: multipolar decomposition. Blue curves: electric (full line) and magnetic (dotted line) dipolar contribution. Red curves: electric (full line) and magnetic (dotted line) quadrupolar contribution. (d) Electric near-field distributions for wavelength position A and B as depicted in (c) corresponding to the excitation of 2 nd ED and 1 stMD resonances, respectively. All calculations are carried out with the permittivity data of crystalline silicon.[START_REF] Green | Solar Energy Materials and Solar Cells[END_REF]

Figure 3 :

 3 Figure 3: (a) Image of a fabricated Si resonator obtained by Scanning Electron Microscopy (SEM). The measurement of the diameter of the Si nanodisk gives a diameter equal to 342 nm. (b,c) Comparison between the measured and the calculated extinction spectra in the visible (b) and near infrared (c) spectrum. Full black line, left scale: experimental measurements of the extinction spectrum; blue full line, right scale: numerical calculations of the extinction cross-section of a silicon cylinder of diameter D = 340 nm and height h = 240 nm with crystalline silicon. 41

Figure 4 :

 4 Figure 4: Comparison of the normalized FWM signal measurement on (a) an unpatterned silicon layer and (b) on a silicon nanodisk of diameter D = 340 nm and height h = 240 nm when the wavelength of the second pump λ 2 is changed. The values in the legend are the different values of λ 2 .

  explains why the two fold enhancement of the internal field intensity |E int | 2 induces a nearly two fold enhancement of the FWM signal. Let us stress that we tuned only λ 2 while the first pump matches the electric dipole resonance. This explains why the gain efficiency brought by the magnetic resonance is limited by 2, but that combining the two resonances leads to a FWM signal enhancement of two orders of magnitude compared with the thin

Figure 5 :

 5 Figure 5: FWM signal normalized by the pump intensities plotted in red dots with respect to the wavelength λ 2 of pump 2 and compared with (a) the spectrum of the internal field enhancement ( |E int | ) 2 calculated with Eq. (2) and with (b) the measured extinction crosssection calculated with the permittivity of crystalline silicon.

  This particle was designed to exhibit two different resonances, namely an electric dipolar resonance associated a magnetic quadrupolar resonance excited by the first pump and a magnetic dipolar resonance FWM experiments with silicon. Higher order resonances could also certainly yield strong FWM signal enhancements, but at the cost of a much larger fingerprint of the photonic nanostructure. The frequency of one pump was swept over the magnetic resonance to verify that the two fold enhancement of the internal field in the volume of the Si particle yielded an almost two fold enhancement of the FWM intensity.This experiment was performed while keeping the frequency of the other pump at the electric dipole frequency. The excitation of the two different resonances leads to an enhancement of the FWM intensity in the visible spectrum of two orders of magnitude compared with a thin Si film of same thickness. Our study proves the importance of internal electric field enhancement yielded by Mie resonances to enhance the frequency conversion efficiency from the near-infrared to the visible spectrum. All-dielectric nanophotonics therefore appears as a promising route for upconverting near infrared signals to higher frequencies and for detecting infrared photons with indirect bandgap semiconductors.MethodsNanodisk fabricationSilicon nanodisks are fabricated by e-beam lithography (EBL) along with reactive-ion etching. A 240 nm thick layer of amorphous silicon is deposited onto a 150 µm thick glass substrate by plasma assisted reactive magnetron sputtering (Buhler, HELIOS). A 60 nm thick layer of a PMMA (poly(methyl-methacrylate)) positive e-beam resist is spin-coated on top of the silicon film. Arrays with a 10 µm pitch of nanocircles with diameter ranging from 250 nm to 360 nm are then patterned by e-beam lithography (Pioneer, Raith). The sample is then developed in a commercial MIBK-isopropanol solution. A layer of nickel, 16 nm in thickness, is then evaporated on top of the sample and acts as a hard mask during the reactive-ion etching step. After the resist development, the remaining PMMA and the nickel layer on top of it are removed in a lift-off process by putting the samples into a ultrasonic bath of ethyl lactate. After this step, the nickel mask remains only on the nanodisks patterned during the e-beam lithography step. The zones of the silicon layer that are not covered by the nickel mask are then etched by reactive-ion etching (MG-200, Plassys) by a gas mixture containing SF 6 , O 2 and CHF 3 and with a pure O 2 plasma. Finally, the remaining nickel mask is removed in an acid solution of HCl and FeCl 3 .

Fig. 1 .

 1 Fig.1. The first pump is a femtosecond pulse centered at the wavelength λ 1 generated by a Ti-sapphire laser (Coherent Chameleon Ultra). The second pump is a femtosecond pulse centered at the wavelength λ 2 , produced by an optical parametric oscillator (OPO) pumped

Acknowledgement

Research conducted within the context of the International Associated Laboratory ALPhFA: Associated Laboratory for Photonics between France and Australia. This work has been carried out thanks to the support of the A*MIDEX project (no. ANR-11-IDEX-0001-02) funded by the Investissements d'Avenir French Government program, managed by the French National Research Agency (ANR). We acknowledge support by the Australian Research Council. Nanofabrication processes were performed in PLANETE CT PACA cleanroom facility (CINaM, Marseille). The authors thank Prof. Yuri Kivshar and Dr Brian Stout for fruitful discussions and valuable comments.


Dipayan Chaudhuri

David Granton

Dominic Xiang Wang

Sharon Einav

Yigal Helviz

Tommaso Mauri

Jean-Damien Ricard

Jordi Mancebo

Jean-Pierre Frat

Sameer Jog

Moderate Certainty Evidence Suggests the Use of High-Flow Nasal Cannula Does Not Decrease Hypoxia When Compared With Conventional Oxygen Therapy in the Peri-Intubation Period

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


Andrea G B Tettamanzi
email: andrea.tettamanzi@univ-cotedazur.fr

David Emsellem
email: david.emsellem@kcitylabs.fr

Célia Da Costa

Alessandro Pereira

Giovanni Venerandi

Fusco

Célia Da Costa Pereira
email: celia.da-costa-pereira@univ-cotedazur.fr

Alessandro Venerandi
email: alessandro.venerandi@univ-cotedazur.fr

Giovanni Fusco
email: giovanni.fusco@univ-cotedazur.fr

Possibilistic

Possibilistic Estimation of Distributions to Leverage Sparse Data in Machine Learning

Keywords: Possibility Theory, Machine Learning, Weakly Supervised Learning

come

Introduction

Supervised learning is the machine learning task of learning a function that maps an input to an output based on example input-output pairs [START_REF] Russell | Artificial Intelligence: A Modern Approach[END_REF]. Each example consists of an input record, which collects the values of a number of input variables, and the associated value of the output variable (also called the supervisory signal). The learnt function can then be used to "predict" the value of the output variable for new unlabeled input records, whose output value is not known.

In many real-world problems, obtaining a fully labeled dataset is expensive, difficult, or outright impossible. An entire subfield of machine learning, called weakly (or semi-) supervised learning has thus emerged, which studies how datasets where the supervisory signal is not available or completely known for all the records can be used for learning. According to a useful taxonomy of classification problems that can arise in that field [START_REF] Hernández-González | Weak supervision and other nonstandard classification problems: A taxonomy[END_REF], four broad classes of problems can be identified:

single-instance, single-label (SISL), which corresponds to the standard setting where all the examples consist of a single instance, to which a single (i.e., certain) class label is assigned; single-instance, multiple-label (SIML), where some examples, consisting of single instances, are assigned a (disjunctive) set of possible class labels, including, to an extreme, the set of all the class labels, which corresponds to the case of a missing supervisory signal; multiple instance, single-label (MISL), when examples may consist of sets of instances, being assigned a single label as a whole; multiple instance, multiple-label (MIML), when, in addition, some examples are assigned a (disjunctive) set of possible values.

This taxonomy of course assumes that the output variable of the underlying objects, which one seeks to predict, can only have a single true value. It should be mentioned that other problems exist, called multi-label classification [START_REF] Tsoumakas | Multi-label classification: An overview[END_REF], where, for a given underlying object, described by an input record, multiple values can be active at the same time, which would then be described by a conjunctive set of output values. In that case, the learnt function is set-valued.

The above taxonomy can be extended to regression or "predictive modeling" problems, where the "label" is a number, ranging on a discrete or continuous interval.

In the framework of an interdisciplinary research project applying machine learning and urban morphology theory to the investigation of the influence of the urban environment on the value of residential real estate [START_REF] Venerandi | A machine learning approach to study the relationship between features of the urban environment and street value[END_REF], we faced the problem of incorporating into a predictive model uncertain information associated with prices, addressing issues of data sparsity, a problem that falls within the SIML category of the above taxonomy. This prompted us to propose an original method to deal with output variable uncertainty in predictive modeling, based on possibility theory, which we present below.

As it has been argued for example by Bouveyron et al. [START_REF] Bouveyron | Robust supervised classification with mixture models: Learning from data with uncertain labels[END_REF], while the problem of noise in data has been widely studied in the literature on supervised learning, the problem of label noise remains an important and unsolved problem in supervised classification. Nigam et al. [START_REF] Nigam | Impact of noisy labels in learning techniques: A survey[END_REF] proposed Robust Mixture Discriminant Analysis (RMDA), a supervised classification whose aim is to detect inconsistent labels by comparing the labels given for labeled data set with the ones obtained through an unsupervised modeling based on the Gaussian mixture model. To solve the problem of automatic building extraction from aerial or satellite images with noise labels, Zhang et al. [START_REF] Zhang | Gis-supervised building extraction with label noise-adaptive fully convolutional neural network[END_REF] propose to capture the relationship between the true label and the noisy label, a general label noise-adaptive (NA) neural network framework consisting of a combination of a base network with an additional probability transition module (PTM) introduced to capture the relationship between the true label and the noisy label. Other researchers prefer to focus on constructing a loss function that is robust to noise [START_REF] Ghosh | Making risk minimization tolerant to label noise[END_REF].

However, the uncertainty brought about by data sparsity in our problem is hard to characterize as a probability distribution, making the application of one of the probabilistic approaches to weak supervised learning proposed in the literature unattractive.

Vannoorenberghe et al., instead, proposed a different approach [START_REF] Vannoorenberghe | Handling uncertain labels in multiclass problems using belief decision trees[END_REF] in which the induction of decision trees is based on the theory of belief functions. In their framework, it is supposed that the training examples have uncertain or imprecise labels. In the same spirit, Quost et al. [START_REF] Quost | Learning from data with uncertain labels by boosting credal classifiers[END_REF] also proposed a belief-function-based framework to be used for supervised learning, in which the training data are associated with uncertain labels. They supposed that each example in the training data set is associated to a belief assignment that represents the actual knowledge of the actual class of the example and used a boosting method to solve the classification problem. Denoeux et al. [START_REF] Denoeux | Handling possibilistic labels in pattern classification using evidential reasoning[END_REF] introduce a category of learning problems in which the labels associated to the examples in the training data set are assessed by an expert and encoded in the form of a possibility distribution. Although this work is very relevant to what we are proposing here, the authors obtain their possibility distributions from human experts, which can be expensive and difficult, whereas the method we propose automatically computes those distributions from data.

Traditionally, housing market valuations are modeled, in a linear fashion, through a combination of intrinsic and extrinsic features evaluated for each dwelling. Although such linear models provide easy-to-read results, they are severely limited as they assume linearity and independence among variables. However, this might not be the case. For example, a specific variable might change its behaviour for different subsets of observations. More recently, researchers have applied Machine Learning (ML) techniques to study the same phenomenon [START_REF] Park | Using machine learning algorithms for housing price prediction: The case of Fairfax County, Virginia housing data[END_REF][START_REF] Chiarazzo | A neural network based model for real estate price estimation considering environmental quality of property location[END_REF][START_REF] Gerek | House selling price assessment using two different adaptive neuro-fuzzy techniques[END_REF]. However, their aim was mainly predictive. Thus, although linearity and independence among variables were tackled through the use of such algorithms, their results lacked interpretability. Finally, the intrinsic/extrinsic dichotomy does not hold when the goal of the analysis is the valuation of urban subspaces (like neighborhoods or street segments) instead of individual dwellings. To tackle these issues, we have devised an approach rooted in Urban Morphology to explain housing values at a fine level of spatial granularity, that of street segments and we designed a sequence of appropriate ML techniques that output interpretable results. To be more specific, the proposed approach, firstly, computes street-based measures of housing values, urban form, functions, and landscape and then models the relationship between them through an ensemble method comprising of Gradient Boosting (GB) [START_REF] Friedman | Greedy function approximation: a gradient boosting machine[END_REF], topological Moran's test [START_REF] Moran | Notes on continuous stochastic phenomena[END_REF], and SHAP [START_REF] Lundberg | A unified approach to interpreting model predictions[END_REF], a recently developed technique to interpret outputs of ML algorithms. The approach has been used to explain the median valuation of street segments in the French Riviera, using housing transactions from the period 2008-2017, through more than 100 metrics of urban form, functions, and landscape.

One difficulty this approach runs into is that transaction data, which are the only source of observations of the output variables (the measures of housing values), are rather sparse at the scale of the street segment. One possible way to overcome this difficulty would be to limit our study to those street segments for which enough observations are available, e.g., at least ten transactions within the observation period; however, this drastically reduces the number of street segments that can be studied and introduces a bias toward neighborhoods having a relatively high turnover. An alternative and more attractive solution, which is the subject of this paper, is to use all available observations, while taking into account the uncertainty brought about by data sparsity.

Essentially, the solution we propose is to model the uncertainty relevant to the output variable within the framework of possibility theory and modify the loss function of the ML technique, used to model the phenomenon, so that it can weight the error based on the uncertainty of the output variable. This approach is very much in the same spirit as the fuzzy loss function proposed by Hüllermeier [START_REF] Hüllermeier | Learning from imprecise and fuzzy observations: Data disambiguation through generalized loss minimization[END_REF][START_REF] Hüllermeier | Learning from imprecise data: Adjustments of optimistic and pessimistic variants[END_REF]. An important advantage of such solution is that it is readily transferable to any supervised or semi-supervised ML technique using a loss function (which is the case for the vast majority of such techniques).

The rest of the paper is organized as follows: Section 2 provides some background on possibility theory, which is required in order to understand the proposed approach; Section 3 states the main question we address, as it emerges from the real-estate price study that motivated our proposal. The proposed solution itself is presented in Section 4, while Section 5 discusses its empirical validation. Section 6 draws some conclusions and proposes some ideas for further research.

Background on Possibility Theory

Fuzzy sets [START_REF] Zadeh | Fuzzy sets[END_REF] are sets whose elements have degrees of membership in [0, 1]. Possibility theory [START_REF] Dubois | Possibility Theory-An Approach to Computerized Processing of Uncertainty[END_REF] is a mathematical theory of uncertainty that relies upon fuzzy set theory, in that the (fuzzy) set of possible values for a variable of interest is used to describe the uncertainty as to its precise value. At the semantic level, the membership function of such set, π, is called a possibility distribution and its range is [0, 1]. A possibility distribution can represent the available knowledge of an agent. π(I) represents the degree of compatibility of the interpretation I with the available knowledge about the real world if we are representing uncertain pieces of knowledge. By convention, π(I) = 1 means that it is totally possible for I to be the real world, 1 > π(I) > 0 means that I is only somehow possible, while π(I) = 0 means that I is certainly not the real world.

A possibility distribution π is said to be normalized if there exists at least one interpretation I 0 s.t. π(I 0 ) = 1, i.e., there exists at least one possible situation which is consistent with the available knowledge. Definition 1. (Possibility and Necessity Measures) A possibility distribution π induces a possibility measure and its dual necessity measure, denoted by Π and N respectively. Both measures apply to a classical set S ⊆ Ω and are defined as follows:

Π(S) = max I∈S π(I); (1)
N (S) = 1 -Π( S) = min I∈ S {1 -π(I)}. (2)
In words, Π(S) expresses to what extent S is consistent with the available knowledge. Conversely, N (S) expresses to what extent S is entailed by the available knowledge. It is equivalent to the impossibility of its complement S-the more S is impossible, the more S is certain. A few properties of Π and N induced by a normalized possibility distribution on a finite universe of discourse Ω are the following. For all subsets A, B ⊆ Ω:

1. Π(A ∪ B) = max{Π(A), Π(B)}; 2. Π(A ∩ B) ≤ min{Π(A), Π(B)}; 3. Π(∅) = N (∅) = 0; Π(Ω) = N (Ω) = 1; 4. N (A ∩ B) = min{N (A), N (B)}; 5. N (A ∪ B) ≥ max{N (A), N (B)}; 6. Π(A) = 1 -N ( Ā) (duality); 7. N (A) > 0 ⇒ Π(A) = 1; Π(A) < 1 ⇒ N (A) = 0; A consequence of these properties is that max{Π(A), Π( Ā)} = 1. In case of complete ignorance on A, Π(A) = Π( Ā) = 1.

Problem Statement

To make this paper self-contained, we briefly recall here some elements, relevant to the solution we are going to describe in Section 4, of the problem that motivated our proposal. The interested reader can refer to [START_REF] Venerandi | A machine learning approach to study the relationship between features of the urban environment and street value[END_REF] for a more detailed explanation.

Pre-processing

Housing transactions of different years and housing typologies cannot be directly compared due to yearly inflation, housing market cycles (e.g., economic recession, upturn), and specific market behaviours affecting different housing types. For example, bigger properties tend to be sold less frequently as they are more expensive and subject to long term investments, while smaller properties, due to their relative lower valuations, tend to be exchanged more easily and tend to be subject of shorter-term investments. The average price per square meter tends also to be structurally higher for small flats for technical reasons (even the smallest flat needs sanitary and cooking equipment, which proportionally weigh more on the average price per surface unit compared to a larger property). The very notion of average price per square meter can thus be challenged when applied to such diverse housing markets. To address these issues, instead of the conventional price per square meter, our method requires to separate the transactions by year and housing type and compute ventiles of prices for each subset year of transaction -housing type. We consider such statistics as appropriate normalized values, which account for different market segments and years, thus making transactions comparable among them.

Computation of the Median of Values

Having classified each transaction in a ventile of value, the next step requires to first assign to each data point the street segment to which they belong and, second, aggregate the information on value at the street level through the computation of a measure of central tendency (i.e., median). Such measure provides information on the central point of the distribution of the value of the housing market, for each street. We perform such computation for the ensemble of each street segment and its immediate neighbouring streets. This for two reasons: firstly, transactions located in streets directly connected to one another tend to have similar valuations (due to the influence of the same location factors, presence of properties at the intersection of several streets segments, etc.); secondly, data on house prices tend to be quite sparse, even for several years, and thus a local interpolation allows us to increase the data coverage. Nevertheless, most street segments end up having less than ten transactions per housing type, which introduces uncertainty into the computation of the median statistics.

Street-Based Metrics of the Urban Environment and Landscape

To characterize the context of each street segment in the most comprehensive way possible, we compute a set of descriptors that quantify aspects of the urban fabric, street-network configuration, functions, housing stock, and landscape. Their definition is out of the scope of this paper.

Proposed Solution

In abstract terms, we can describe the problem as follows. We are given a sample of observed variates x 1 , x 2 , . . . , x n and a number of probability distributions (the hypothesis space). We want to assess the possibility degree of these probability distributions based on the given sample, i.e., the degree to which they are compatible with the observations.

Possibility of a Price Distribution

We can limit ourselves to parametric families of distributions. In this case, this problem can be described as a sort of possibilistic parameter estimation. In the Incertimmo project, we consider distributions described by three deciles: (d 1 , d 5 = median, d 9 ). We recall that the ith decile d i is the smallest number

x satisfying Pr[X ≤ x] ≥ i 10 . By definition, d 1 ≤ d 5 ≤ d 9 .
Furthermore, for a probability distribution defined in the [a, b] interval, we know that 10% of the probability mass is in the [a, d 1 ] interval, 40% in the (d 1 , d 5 ] interval, 40% in the (d 5 , d 9 ] interval, and 10% in the (d 9 , b] interval. We make the additional simplifying assumption that the probability mass is uniformly distributed within each of the above intervals. While this might look like a very restrictive assumption, on the one hand it is motivated by the type of qualitative descriptions of price distribution that are of interest to the human geographers (i.e., by the application at hand) and, on the other hand, could easily be relaxed by selecting other parametric families of distributions without serious consequences on the proposed approach.

This yields a parametric family of probability distributions on the [a, b] interval whose density (in the continuous case) or probability (in the discrete case) function is

f (x; d 1 , d 5 , d 9 ) =          1 10(d1-a) , a ≤ x ≤ d 1 ; 2 5(d5-d1) , d 1 < x ≤ d 5 ; 2 5(d9-d5) , d 5 < x ≤ d 9 ; 1 10(b-d9) , d 9 < x ≤ b. (3)
The degenerate cases that arise when d 1 = d 5 , d 5 = d 9 , or d 9 = b are treated by adding to the result the contributions of the lines whose condition becomes empty.

In the specific application described in Section 3, as we have seen, the observed variates are ventiles of the general distribution of housing prices, taking up values in the discrete set {1, 2, . . . , 20}.

Given the sample x 1 , x 2 , . . . , x n , it is easy to compute the probability that it is produced by a given distribution of the parametric family, having parameters (d 1 , d 5 , d 9 ). This is

Pr[x 1 , x 2 , . . . , x n | (d 1 , d 5 , d 9 )] = n i=1 f (x i ; d 1 , d 5 , d 9 ). ( 4
)
This probability is in fact a likelihood function over the distribution of price ventiles, which we will denote by L(d 1 , d 5 , d 9 ). The link between likelihoods and possibility theory has been explored in [START_REF] Dubois | A semantics for possibility theory based on likelihoods[END_REF]. The main result of that study was that possibility measures can be interpreted as the supremum of a family of likelihood functions. It should be stressed that this is an exact interpretation, not just an approximation. Based on this result, we transform the likelihood function L into a possibility distribution over the set of parametric probability distributions of the form (d 1 , d 5 , d 9 ).

Notice that the parameters d 1 , d 5 , d 9 of the probability distributions are the elementary events of this possibility space. A possibility distribution over that space is obtained by letting

π(d 1 , d 5 , d 9 ) = L(d 1 , d 5 , d 9 )/L max , (5)
where

L max = max 1≤x≤y≤z≤20 L(x, y, z),
with x, y, z ∈ {0, 1, . . . , 20}, so that all the maximum-likelihood probability distributions have a possibility degree of 1, thus yielding a normalized possibility distribution. Alternatively, the logarithm of the likelihood could be used instead.

We restrict the values of the three parameters to the set {0, 1, . . . , 20}, because prices are relative and expressed in ventiles in the application at hand. Of course, for other applications, different ranges of values should be considered. Now, specific values or intervals of one parameter will constitute complex events, i.e., sets of elementary events, and their possibility and necessity measures can be computed as usual, as the maximum of the possibilities of the distributions that fit the specification and 1 -max of all the others, respectively. For instance, the possibility measure over the median price ventile (d 5 ) of a given street segment will be given, for all d 5 ∈ {1, . . . , 20}, by

Π(d 5 ) = max 1≤x≤d5≤y≤20 π(x, d 5 , y), (6)
where π(•, •, •) is defined as in Equation 5.

Loss Function under Possibilistic Uncertainty

The error made by the model predicting that the median price for a street segment is in ventile ŷ when all we know is the possibility distribution π over the probability distributions of the prices for that segment, can be defined as

L(ŷ, π) = 1 0 min Π(y)≥α (ŷ -y) 2 dα, (7)
where Π(y) is the possibility measure of the distributions having y as their median. Equation 7 is based on an underlying square error function e(y) = (ŷ -y) 2 , but it could be easily generalized to use an arbitrary error function.

In practice, if Λ = (0 = λ 1 , λ 2 , . . . , 1) is the list of possibility levels of π, such that ∀i > 1, ∃(z, y, z) : π(z, y, z) = λ i , Equation 7 can be rewritten as

L(ŷ, π) = Λ i=2 (λ i -λ i-1 ) min Π(y)≥λi (ŷ -y) 2 . (8)
This loss function has been coded in Python in such a way that it could be provided as a custom evaluation function to an arbitrary machine learning method offering this possibility. Most of the loss computation requires iterating through all price distributions of the parametric family (with the three parameters 0 ≤ d 1 ≤ d 5 ≤ d 9 ≤ 20, there are 1,540 of them). To optimize performance, we pre-computed the loss function into a lookup table. Since we are using Gradient Boosting Gradient descent (Newton version in XGBoost), we have provided also functions that return its gradient and Hessian.

Experiments and Results

In this section we report the experiments we carried out to validate our method.

We use real-world data consisting of all the housing transactions made on the French Riviera over the period 2008-2017. 1 Each record contains detailed intrinsic features of the dwelling that was sold/bought, including its address and the price paid. From these data, we compiled a dataset whose records correspond to street segments, described by more than 100 metrics of urban form, functions, and landscape, and a distribution of price ventiles for each type of dwelling.

Our goal is to compare the performance of a predictive model trained on this dataset, where the labels are uncertain, due to sparseness of transaction data, to the performance of a predictive model trained on a dataset where the labels are certain (in our case, estimated based on a sufficient number of transactions). If the model trained under uncertainty is able to obtain results similar to those of the model trained without uncertainty, this can be taken as evidence that our method is successful at compensating for the loss of transaction data.

Experimental Protocol

We proceeded as follows. From our dataset, we extracted the set of street segments having at least 10 recorded transactions. These are the street segments for which we consider that the distribution of prices can be estimated in a reliable way. Let us call this dataset D. We then constructed a second dataset D by randomly subsampling the transactions from the street segments of D, in such a way as to obtain a similar distribution of the number of transactions per street segment as in the full dataset (i.e., the dataset including also segments having fewer than 10 recorded transactions). For example, if in the full dataset 15% of the street segments has more than 10 recorded transactions (by the way, this subset of the full dataset is what we have called D), then also 15% of the street segments of D has more than 10 transactions. In general, if the percentage of street segments in the full dataset Table 1. Distribution of available transactions per street segment for different subsampling rates. Column "≥ 10" gives the number of street segments with at least 10 transactions, " [START_REF] Chiarazzo | A neural network based model for real estate price estimation considering environmental quality of property location[END_REF][START_REF] Ghosh | Making risk minimization tolerant to label noise[END_REF]" the number of street segments with 3 to 9 transactions, and so on. that have n transactions is x, then D is samples so that the percentage of street segments in D that have n transactions be as close as possible to x. This way, we can say that D is to D as D is to the set of all street segments in the study area and, as a consequence, any observation about the predictive power on D of models trained on D can provide, by proportional analogy, an indication of the predictive power on the full dataset of models trained on D. Figure 1 graphically illustrates this sampling mechanism and Table 1 provides some statistics about the dataset D that we obtain depending on the chosen sampling rate. Notice that the two datasets D and D have the same size (in our case, 14,317) and consist of exactly the same street segments; what differs between them is the number of transactions available to estimate the distribution of prices for each street segment. In dataset D, the distribution is known precisely, and it has the form (d 1 , d 5 , d 9 ). In dataset D , instead, all is known is a possibility distribution π, as explained in Section 4.

Rate

Validation of the Possibility Distribution

To show that the possibility distribution π defined as per Equation 5, as well as its associated possibility and necessity measures, does indeed qualitatively describe the actual price distribution, we studied the possibility (computed according to Equation 6) of the observed median price ventile m of street segments. Ideally, Π(m) should be 1 for every street segment, if π perfectly described how the prices are distributed.

Figure 2 shows the probability distribution of Π(m) for different sampling rates of the set of transactions. We can observe that when the possibility distributions π of transaction prices for each street segment is constructed using all the recorded transactions (which are at least 10 for any one of the 14,317 street segments considered for this study), the median is assigned a possibility of 1 for most street segments, with some exceptions, which, upon inspection, turned out to be street segments whose price distribution is not unimodal. Since the parametric family of distribution used to fit the actual price distributions is unimodal, the most likely values for their parameters are those that make d 5 correspond to one of the modes. This is an intrinsic limitation introduced by the particular choice of a unimodal family of distribution, which was made to simplify the geographical interpretation of the result; however, despite this limitation, the results of the study seem to confirm the validity of the method used to construct the possibility distributions. Unsurprisingly, as the sampling rate decreases, the number of street segments for which the median is assigned a high possibility decreases.

Empirical Test of the Method

To conduct our tests, we selected, among all possible predictive modeling methods, XGBoost [START_REF] Chen | XGBoost: A scalable tree boosting system[END_REF], which is the one that gave the best results when applied to dataset D based on a critical comparison and benchmarking of the most popular methods available. The rationale of this choice is that we wanted our solution to "prove its mettle" on a very challenging task, namely to prevent the degradation of the accuracy of the strongest available method when the supervisory signal becomes uncertain.

We trained a model to predict the median of ventile prices (i.e., d 5 ) on dataset D by using XGBoost regression with the standard loss function and we trained another model on dataset D by using XGBoost regression instrumented with the proposed possibilistic loss function. We compare the results given by these two models when applied to a test set consisting of street segments not used to train the models. We treat the model trained on dataset D as the ground truth and we measure the deviation of the model trained on D from such a target.

As a measure of prediction error, we compute the RMSD of the median (d 5 ) predicted by the model trained on D for each segment, with respect to the median of that segment in D. We used a sampling rate of 20% to generate D from D, i.e., D contains only 20% of the transactions available in D. We split D into a training set containing 80% of the street segments and a test set containing the remaining 20% of the street segments. After training the model on the training set, we obtain an RMSD of 2.778556 on the test set. Figure 3 shows a plot of predicted vs. actual price ventiles for the test set. For comparison, XGBoost using the standard loss function based on MSE trained on 80% of D (therefore, with full information), gives an RMSD of 2.283798 when tested on the remaining 20%. In other words, the possibilistic loss function allows the prediction error to increase by less than 22%, even though 80% of the transactions were removed from the dataset!

Conclusion

We have proposed a method based on possibility theory to leverage sparse data, which can be combined with any machine learning method to approach weakly supervised learning problems. The solution we propose constructs a possibilistic loss function, which can then be plugged into a machine learning method of choice, to account for an uncertain supervisory signal.

Our solution is much in the same spirit as the fuzzification of learning algorithms based on the generalization of loss function proposed in [START_REF] Hüllermeier | Learning from imprecise and fuzzy observations: Data disambiguation through generalized loss minimization[END_REF], in that it pursues model identification at the same time as data disambiguation, except that in our case ambiguity (i.e., uncertainty) affects the output variable only, which is only partially observable in the available data, while all the input variables are perfectly known and, thus, non-ambiguous. Furthermore, as in [START_REF] Hüllermeier | Learning from imprecise and fuzzy observations: Data disambiguation through generalized loss minimization[END_REF], the predictive model is evaluated by looking at how well its prediction fits the most favorable instantiation of the uncertain labels of the training data (this is the sense of the minimum operator in the possibilistic loss function definition).

The development of the method we presented has been motivated and driven by a very specific application, namely by the need to leverage sparse data in a human geography setting. However, its working principle is quite general and could be extended to suit other scenarios. Indeed, distilling a completely general method is the main direction for future work.

Fig. 1 .

 1 Fig. 1. A graphical illustration of the sampling mechanism used to create subsamples of the original dataset.

Fig. 2 .

 2 Fig. 2. Distribution of Π(m) for different subsampling rates.

Fig. 3 .

 3 Fig. 3. Results obtained on the test set by applying XGBoost regression to D with the possibilistic loss function.

Extracted from the PERVAL database, https://www.perval.fr/.

Acknowledgments

This research was funded by the IDEX "UCA JEDI", within the AAP Partenariat 2016 (action 6.5).

Andrea Tettamanzi has been supported by the French government, through the 3IA Cte dAzur "Investments in the Future" project managed by the National Research Agency (ANR) with the reference number ANR-19-P3IA-0002.

The authors would like to thank Denis Overal, director of the R&D at Kinaxia, for his support and insightful suggestions.


Tudor-Bogdan Airimitoaie

Patrick Lanusse

Asma Achnib

Robust Anticipative Controller Design with Application to A High Dynamic Engine Testbed

Keywords: robust control, discrete-time systems, generalized predictive control, feedforward anticipative filter, preview systems

This paper presents a robust controller design method for reference tracking in multi-input multi-output preview systems. In the context of preview systems, it is supposed that future values of the reference signal are available a number of time steps ahead. The objective of the controller being to minimize a quadratic error between the reference and the system's output, the optimal solution needs to take into account the known future values of the reference. Furthermore, it is desired to maintain the control signal at an acceptable level. The feedforward preview is obtained by solving a mixed L 2 /L ∞ optimisation, where the L ∞ constraint is used to reduce the level of the control. The proposed solution combines a robust feedback controller with a feedforward preview filter. The feedback controller's purpose is to assure robustness of the closed-loop system to model uncertainties and is not detailed here. The focus of this paper is on the design of the feedforward preview filter taking into account that a feedback controller is present. The proposed solution is validated in simulation on a high dynamic engine testbed.

INTRODUCTION

In a number of automatic control applications (known also as preview systems), future values of the reference signal are available in advance. Examples can be found in aeronautical systems (trajectory tracking of unmaned aerial vehicles in [START_REF] Park | A new nonlinear guidance logic for trajectory tracking[END_REF]), autonomous vehicles (explicit path tracking in [START_REF] Shin | Explicit path tracking by autonomous vehicles[END_REF]) or mobile robots [START_REF] Freitas | Terrain model-based anticipative control for articulated vehicles with low bandwidth actuators[END_REF]. In other applications, disturbance signals can also be known beforehand, such as in water level control (for irrigation and drainage ducts in [START_REF] Schuurmans | Simple water level controller for irrigation and drainage canals[END_REF]; [START_REF] Li | Offtake feedforward compensation for irrigation channels with distributed control[END_REF]). Therefore, it is of interest to propose appropriate control algorithms that take into account this information.

The present paper is focused on the development of discretetime robust control methods for preview systems. Throughout this paper, it is considered that future values of the reference signal are known a number of time steps ahead. The objective of preview is that of designing control algorithms capable of minimizing a quadratic criterion on the error between the plant output and the known future reference. As such, the controller anticipates the future actions that need to be taken with respect to changes of the reference signal and the plant output reacts to these changes before they occur.

The concept of anticipative control (known also as preview control) was proposed initially in [START_REF] Sheridan | Three models of preview control[END_REF] and further developed in [START_REF] Tomizuka | Optimal continuous finite preview problem[END_REF]. A complete review of anticipative control is presented in [START_REF] Liao | Robust preview tracking control for a class of uncertain discrete-time systems[END_REF].

Nevertheless, only a few number of research articles have attempted to deal with the problem of robust anticipative control.

In [START_REF] Makarov | Modeling and preview H ∞ control design for motion control of elastic-joint robots with uncertainties[END_REF] an anticipative control problem is considered in the context of motion control of robots. The proposed solution is based on the H ∞ design methodology. Weighting norms are used to assure robustness of the feedback control law.

A free-weighting matrices technique is used with Lyapunov stability theory to derive robust anticipative controllers in [START_REF] Liao | Robust preview tracking control for a class of uncertain discrete-time systems[END_REF]; [START_REF] Li | H ∞ preview control of a class of uncertain discrete-time systems[END_REF]. Both reference anticipative tracking and disturbance rejection are treated; however, an initial design step involves a quadratic criterion which depends on weighting matrices that can be chosen arbitrarily. The proposed method is dependant on the choice of these weighting matrices, which makes it difficult to interpret the obtained results.

The contribution of this paper is to extend the results presented in Achnib et al. (2018a,b) to multi-input multi-output (MIMO) systems and to present a complete theoretical development of the proposed approach. As in Achnib et al. (2018a,b), the robustness considerations take into account both parametric uncertainties and control effort limitations.

A two degrees of freedom feedforward-feedback control schema represents the basis for the proposed controller. A feedforward filter with anticipative behaviour is introduced to take into account future values of the reference signal. It is supposed that the time window of known future reference values is known. The preview window's length of the anticipative feedforward filter can be adjusted accordingly. The feedforward part is obtained by solving a frequency domain optimisation problem with mixed L 2 performance and L ∞ constraints for robustness. The interesting aspect of the proposed approach is that the anticipative feedforward action can be added on top of any robust feedback controller.

Simulation results on a 2-inputs 2-outputs high dynamic engine testbed from [START_REF] Lanusse | Development of a fractional order based mimo controller for high dynamic engine testbeds[END_REF] validate the proposed approach. Testbeds are used in the automobile industry to tune and assess automotive engines. As more High Dynamic (HD) scenarios need to be simulated, the testbeds require MIMO and robust control systems. Furthermore, to improve the accuracy of the results, it is desired to follow the test scenarios as closely as possible. Since the configuration of each scenario is known in advance, it is possible to use this information in a preview controller in order to minimize the tracking error.

The paper is organised as follows. The problem formulation that is addressed in the present paper is described in Section 2. Section 3 introduces the factorisation of the preview feedforward filter that will be used in this paper. Then the optimisation of the unknown part of the filter is presented in Section 4. A solution to reduce the size of this filter is described in Section 5. Simulation results on the high dynamic engine testbed are given in Section 6. Section 7 concludes this paper.

PROBLEM FORMULATION

In this paper we use square discrete-time transfer function matrices that we denote with upper-case letters. Signals are denoted with small-case letters in the time-domain and uppercase letters in the frequency-domain. Z-transform variable z is used in the frequency-domain, while in the time-domain we use the unit advance and delay operators q and q -1 . RH m×m Figure 1 gives a representation of the problem addressed in this paper. The plant, represented by G(z) ∈ RH m×m ∞ belonging to a model-set M G , is uncertain (or with slowly time-varying parameters) with delays on inputs and/or outputs. The sampling time is T s and the sampling frequency is f s . K(z) denotes the robust feedback controller designed to assure the stability of the feedback loop in the presence of parametric uncertainties and disturbances

d u (t) ∈ R m and d y (t) ∈ R m . G n (z) is the nominal plant model chosen from M G such that G n (z) ∈ RH m×m ∞ has only input delays. If the true plant G(z) = G n (z)
and in the absence of disturbances d u (t) and d y (t), it can be seen from Fig. 1 that ε(t) = 0 and u f b (t) = 0. F(z) is the preview feedforward filter that needs to be designed such that the output y(t) ∈ R m of the plant G(z) follows the reference y r (t) and minimises the ∞-norm of the error

ε y (t) = y r (t) -y(t).
(1)

It is assumed that at any time a + d samples of the reference signal y r (t) ∈ R m (see also Fig. 1) are known beforehand, where d is the maximum input delay in G n (z); however, no other assumption is made about the dynamics that generate y r (t).

The term "preview" is used to indicate that the filter has terms depending on future values of its input. The interpretation of "preview" is better understood in the time-domain, where F(q) has terms in q -1 and also in q which depend on future values of this filter's input. This is possible in applications for which future values of the reference or of the disturbances are known beforehand (see discussion in the introduction section).

For robustness, the variability of the true plant model and the limits of the control actuator have to be taken into account in the design of F(z). From the schema in Fig. 1, it can be seen also that the control u(t) is given by the sum of the feedback and feedforward control actions:

u(t) = u f b (t) + u f f (t).
Thus, the design of F(z) has to take into account also the effects of the feedback loop.

FACTORISATION OF THE FEEDFORWARD PREVIEW FILTER

In this section, we introduce a factorisation of F(z) that is useful in the development of the algorithm. We begin by recalling the inner-outer factorisation.

Lemma 1. (Inner-Outer factorisation [START_REF] Vidyasagar | Control System Synthesis: A Factorization Approach[END_REF]). Let

G(z) ∈ RH p×m ∞ have full rank for all z = e jω . Then G(z) has inner-outer factorisation G(z) = G i (z)G o (z), where G i (z) ∈ RH p×n ∞ , n = min(p, m), is said to be inner if G T i (z -1 )G i (z) = I and G o (z) ∈ RH n×m ∞ is outer if it has a right inverse G -1 o (z) which is analytic outside the unit circle. If p ≥ m, then G o (z) is square, whereas if p ≤ m then G i (z) is square.
The proof is given in [START_REF] Vidyasagar | Control System Synthesis: A Factorization Approach[END_REF]. Using this lemma, the nominal model G n (z) with input delays has an inner-outer factorisation

G n (z) = G i (z)G o (z)G d (z), (2)
where G i (z) is the inner factor, G o (z) is the outer factor and G d (z) is a m × m diagonal matrix regrouping the input delay terms.

This suggests the following factorisation of the preview feedforward filter F(z)

F(z) = F o (z)T F (z), with F o (z) = G d (z -1 )G -1 o (z). (3) In the previous equation, G d (z -1
) is a preview filter depending on future values of its input.

Using the factorisations of G n (z) and F(z), the schema in Fig. 1 can be simplified as in Fig. 2. In (3), T F (z) is a square diagonal preview filter. Each diagonal term of T F (z) can be written:

T F (z)(k) = t F k,-s z -s + . . . + t F k,-1 z -1 + t F k,0 + t F k,1 z + . . . + t F k,a z a , ( 4
)
where k indicates the diagonal element (1 ≤ k ≤ p). In the previous equation, s and a give, respectively, the size of the storage and the anticipation horizon of filter T F (z). The rest of this paper is focused on the design of the preview filter T F (z) in order to minimize the ∞-norm of the error ε y (t) under input u FF (t) constraints and uncertainties of the real plant model G(z).

DESIGN OF THE PREVIEW FEEDFORWARD FILTER

This section describes the method used to find the parameters of the preview filter T F (z). Let denote

θ k = [t F k,-s , . . . ,t F k,-1 ,t F k,0 ,t F k,1 , . . . ,t F k,a ] T (5)
the adjustable parameters of the diagonal element in (4) and

θ = [θ T 1 , . . . , θ T p ] T (
6) the vector of all adjustable parameters. Two closed-loop transfer function matrices in Fig. 2 are of interest. Let us define, in the frequency-domain, the transfer matrix from reference to control input as

H uy r (z) = = (I m + K(z)G(z)) -1 G i (z -1 )G -1 o (z) + K(z) T F (z), (7
) and the transfer matrix from reference to system output as

H yy r (z) = G(z)H uy r (z), (8)
where I m is the identity matrix of size m.

Using ( 1) and ( 8), the following objective function is defined

J F (θ ) = sup t ε y (t) (9)
= sup t (I p -H yy r (q)) y r (t) , ∀G(q) ∈ M G .

For finite energy reference signals y r (t), a least upper bound of the previous objective function can be obtained by using the 2-norm 1 (see [START_REF] Doyle | Feedback control theory[END_REF]):

J F (θ ) ≤ I p -H yy r (e jωT s ) 2 y r (t) 2 . ( 11
)
As stated in section 2, the dynamics that generate y r (t) are unknown, as such the previous objective should be satisfied for any y r (t), thus only the first term in the right-hand side of the inequality (11) can be optimised. The following minimisation problem is proposed for the unknown parameters θ : θ = arg min

θ I p -H yy r (e jωT s ) 2 , ∀G(z) ∈ M G . ( 12
)
Proposition 2. For diagonal preview filter T F (z), the minimisation problem given in ( 12) can be decomposed in p independent minimisations, one for each diagonal element of T F (z) as in θk = arg min

θ k I p (:, k) -H yy r (e jωT s )(:, k) 2 , ( 13
)
where the notation (:, k) is used to indicate the k-th column of the matrix.

Proof. Introducing the p-inputs and p-outputs transfer matrix

H yy f (z) = G(z) (I m + K(z)G(z)) -1 G i (z -1 )G -1 o (z) + K(z) , in (7), one obtains H yy r (z) = H yy f (z)T F (z), (14)
1 For discrete-time operators, the frequency response is obtained by replacing q through e jωTs = e j2π f / fs .

where

T F (z) = diag(T F (z)(i)) =    T F (z)(1) . . . T F (z)(p)    . (15)
After multiplication on the right-hand side of ( 14), one obtains H yy r (z) = H yy f (z)(:, 1)T F (z)(1), . . . ,

H yy f (z)(:, p)T F (z)(p) . ( 16
)
It is straightforward to see that the L 2 -norm in (12) simplifies to the sum of the L 2 -norms for each column element as in ( 13) and that the optimum θ is given by the collection of the θk .

Rewritting (4) as T F (z)(k) = φ T (z)θ k , where θ k has been defined in (5) and φ (z) = [z -s , . . . , z -1 , 1, z, . . . , z a ] T , column k of H yy r (z) is given by

H yy r (z)(:, k) = H yy f (z)(:, k)φ T (z)θ k , (17)
which shows that it is linear in the parameters vector θ k .

Using proposition 2, we define p constrained optimisation problems (one for each θ k ) as θk = arg min

θ k I p (:, k) -H yy r (e jωT s )(:, k) 2 , (18a) s.t. a ∑ l=-s θ k (l) = 1 (18b)
and

W uy r (e jωT s )(i, k)H uy r (e jωT s )(i, k) ∞ ≤ 1, (18c) ∀ 1 ≤ i ≤ m, ω ∈ [0, π f s ] and G(z) ∈ M G .
Two constraints are added to the optimisation problem in ( 18). The first one, (18b), is introduced to ensure that the steady state gain of T F (z) has unit value.

In the second one, (18c), the weighting function W uy r (z) is used to introduce a frequency constraint on the control input u(t). In practice, this is necessary in order to limit the control action at frequencies where the plant model has very low gain.

REDUCING THE NUMBER OF PARAMETERS OF THE PREVIEW FEEDFORWARD FILTER

Often in applications, the dynamics of the preview feedforward filter need not be as fast as the other blocs in Fig. 2. It is therefore recommended to choose a slower sampling period for the preview feedforward filter T F (z) in order to reduce the size of vector θ i .

Let r denote a sampling time multiplier for the preview feedforward filter. The Z-transform variable and the unit time advance operator for the sampling period T s r = r • T s become, respectively, z r = z r and q r = q r . Then T F (z r ) denotes the anticipative part of the feedforward filter that is sampled at r • T s seconds with diagonal element:

T F (z r )(k) = t F k,-s z -s r + . . . + t F k,-1 z -1 r + t F k,0 + t F k,1 z r + . . . + t F k,a z a r , (19
) Taking into account the sampling period multiplier r, the sampling frequency for the T F (z r ) filter in ( 19) is f s r = f s r . As such, z r becomes in the frequency domain e jω r (rT s ) , where ω r ∈ 0, π f s r .

To reduce aliasing effects due to T F (z) which is sampled at a slower frequency, an anti-aliasing filter H f (z) is introduced at the sampling period T s as in Fig. 3. The cut-off frequency of H f (z) is 0.8 f s 2r . Fig. 3 shows also the separation of sampling times in the control schema. The vector valued signal y a (t) is obtained from y a r (t) by upsampling at r times the initial sampling rate and introducing zeros at the new sampling times. Then a filtering through

H up (z) = r-1 ∑ n=0 z -n , ( 20
)
allows to obtain, in the time-domain, the same result as the initial signal:

y a ((n • r + v) • T s ) = y a r (n • T s r ), ∀ 0 ≤ v ≤ r -1.
The optimisation problem in ( 18) involves the frequency responses of filters H yy r and H uy r which now include both terms sampled at T s and T s r .

Let's denote H uy a (z) the transfer from Y a (z) to U(z). The following corollary gives the equation defining this transfer. Lemma 3. The frequency response of the transfer from Y r (z r )

to U(z), ∀ω ∈ 0, π f s r , is given by H uy r (e jωT s ) = H uy a (e jωT s ) 1 r r-1 ∑ n=0 e -jωT s n T F (e jωrT s ). ( 21)

Proof. Using the definition of H uy a (z), one obtains the frequency response of this transfer function matrix as H uy a (e jωT s ), ∀ω ∈ [0, π f s ]. One can restrict its evaluation over the frequencies given by ω r ∈ 0, π f s r :

U(e jω r T s ) = H uy a (e jω r T s )Y a (e jω r T s ), ∀ω r ∈ 0, π f s r .

At the same time, the frequency response from Y r (z r ) to Y a r (z r ) is given by

Y a r (e jω r (rT s ) ) = T F (e jω r (rT s ) )Y r (e jω r (rT s ) ), ∀ω r ∈ 0, π f s r (23)
Using the Discrete Fourier Transform (DFT) definition for N samples of the vector valued signal y a r (t):

Y a r (k) = 1 N N-1 ∑ n=0 y a r (n)e -j2π kn N , ∀ 0 ≤ k ≤ N -1. (24)
Let's denote by y a r up (t) the result of the up-sampling of y a r (t):

y a r up (n • r + v) = y a r (n), if v = 0 0, otherwise. (25)
One can compute the DFT of rN samples of y a r up (t):

Y a r up (k) = 1 rN rN-1 ∑ n=0 y a r up (n)e -j2π kn rN , ∀ 0 ≤ k ≤ rN -1 (26)
Introducing the change of notation n = w • r + v, with 0 ≤ w ≤ N -1 and 0 ≤ v ≤ r -1, and using (25) for 0 ≤ k ≤ N -1:

Y a r up (k) = 1 rN N-1 ∑ w=0 r-1 ∑ v=0 y a r up (w • r + v)e -j2π k(w•r+v) rN (27) = 1 rN N-1 ∑ w=0 y a r (w)e -j2π kwr rN (28) = 1 r Y a r (k), ∀ 0 ≤ k ≤ N -1. (29)
From ( 20), ( 22), ( 23), and (29), one obtains U(e jω r T s )/Y r (e jω r (rT s ) ) = U(e jω r T s )/Y a (e jωT s )• Y a (e jωT s )/Y a r (e jω r (rT s ) ) •Y a r (e jω r (rT s ) )/Y r (e jω r (rT s ) )

= H uy a (e jω r T s ) 1 r r-1

∑ n=0 e -jω r T s n T F (e jω r (rT s ) ), ∀ω r ∈ 0, π f s r , which proves Lemma 3. Corollary 4. The frequency response of the transfer from Y r (z r )

to Y (z) for ω ∈ 0, π f s r is given by

H yy r (e jωT s ) = = G(e jωT s )H uy a (e jωT s ) 1 r r-1 ∑ n=0 e -jωT s n T F (e jωrT s ). (30)
Using ( 21) and ( 30), the optimisation problem ( 18) can be solved for ω ∈ 0, π f s r .

APPLICATION TO A HIGH DYNAMIC ENGINE TESTBED

In this section, we present simulation results on a high dynamic engine (HDE) testbed obtained using the preview feedforward approach described in this paper. First, section 6.1 describes the testbed used for simulations. Finally, section 6.2 shows results obtained using the preview feedforward filter. The testbed and the feedback controller have been previously presented in [START_REF] Lanusse | Development of a fractional order based mimo controller for high dynamic engine testbeds[END_REF].

Description of the High Dynamic Engine Testbed

Dynamic combustion engine testbeds are used for many purposes, such as calibration, driving cycles, control development and research. The simulation system used in this paper is based on the HDE testbed used in [START_REF] Lanusse | Development of a fractional order based mimo controller for high dynamic engine testbeds[END_REF], which is available at the PRISME laboratory in Orléans, France. It can be used to test combustion, hybrid or electric vehicle engines.

A HDE testbed consists of a vehicle engine to be tested and a dynamic mechanical load provided by an electrical machine.

Here, the test engine is a Peugeot S.A. EB2, 1.2 liter 3-cylinder Spark Ignition (SI) engine with a maximum power of 60 kW, a maximum torque of 116 Nm and a maximum speed of 4500 rpm. This engine is controlled by a fully open engine control unit (ECU) with T s = 5 ms sampling time. The load is provided by an asynchronous MDA 250 KW machine driven by ABB technology (see also schema in Fig. 4).

The outputs to be controlled are the engine speed and torque (see Fig. 4). The control inputs are the desired current that feeds the electrical machine and the throttle (accelerator pedal) position. Thus, the HDE testbed is a 2-inputs 2-outputs system: 20 operating points have been taken into account for system identification and control. These points have been chosen such that the engine speed varies from 2150 to 3800 rpm (rotations per minute) and the torque from 15% to 75% of the highest torque that the engine can provide for each speed (more details are given in [START_REF] Lanusse | Development of a fractional order based mimo controller for high dynamic engine testbeds[END_REF]). An operating point at 2590 rpm and 52 Nm has been arbitrarily chosen to define a nominal model. A frequency domain identification of the 21 models has been done and Fig. 5 shows the magnitude and phase plots of all models. One can remark significant variation in both magnitude and phase. The nominal model transfer function matrix is given by (32).

Y 1 (z) Y 2 (z) = G 11 (z) G 12 (z) G 21 (z) G 22 (z) U 1 (z) U 2 (z) (31)

Simulation results using the proposed anticipative approach

The preview feedforward filter approach is evaluated in simulation using the HDE testbed and CRONE feedback controller developed in [START_REF] Lanusse | Development of a fractional order based mimo controller for high dynamic engine testbeds[END_REF].

A sampling period multiplier r = 20 has been chosen. The anticipation order is a = 5 while the storage coefficient of the feedforward filter is s = 4. The optimisation has been done using fmincon in Matlab. All 20 models plus the nominal model have been used. The T F (z r ) filter that satisfies the optimisation problem is given in (33). T F 11 (z r ) =0.0183z 5 r + 0.0464z 4 r + 0.0955z 3 r + 0.202z 2 r + 0.375z 1 r + 0.0852 + 0.151z -1 r + 0.00637z -2 r + 0.0495z -3 r -0.0289z -4 r , (33a)

T F 12 (z r ) =0, T F 21 (z) = 0, (33b)
T F 22 (z r ) =0.178z 5 r + 0.0997z 4 r + 0.126z 3 r + 0.0495z 2 r + 0.0899z 1 r + 0.0135 + 0.0672z -1 r + 0.109z -2 r + 0.0953z -3 r + 0.172z -4 r . (33c)
Figures 6, 7, 8, and 9 show simulation results using the preview feedforward + robust feedback control presented in this paper.

The reference signals y r (t) are the same as in the previous subsection. The simulations show the variation of control inputs and system outputs around the equilibrium of each linearised model. It can be observed that the two outputs anticipate the changes of the references before they actually occur. Furthermore, the feedforward reduces the chattering in the outputs and input signals with respect to the feedback controller alone.

CONCLUSIONS

In this paper, we have presented an algorithm for the design of preview feedforward filters for MIMO systems. The proposed approach considers robustness to plant model uncertainties and control level constraints. Simulation results on a high dynamic engine testbed validate the proposed control. As shown in [START_REF] Lanusse | Development of a fractional order based mimo controller for high dynamic engine testbeds[END_REF], a judicious choice of a simplified nominal model allows to synthesise a feedback CRONE controller with improved performances for the engine speed control and an improved decoupling from the torque control. The resulting control-system will be assessed using a cycle based on the dynamic part of the non-road transient cycle (NRTC) test.

G 11 (z) =z -3 -0.644 + 1.98z -1 -2.12z -2 + 0.669z -3 + 0.294z -4 -0.173z -5 1 -4.29z -1 + 7.56z -2 -6.82z -3 + 3.15z -4 -0.591z -5 (32a)

G 12 (z) =z -15 0.148 -0.508z -1 + 0.665z -2 -0.383z -3 + 0.064z -4 + 0.0135z -5 1 -4.59z -1 + 8.64z -2 -8.35z -3 + 4.16z -4 -0.851z -5 (32b)

G 21 (z) =z -2 0.0919 + 0.00408z -1 -0.293z -2 + 0.223z -3 -0.0255z -4 -0.00105z -5 1 -4.01z -1 + 6.56z -2 -5.46z -3 + 2.27z -4 -0.371z -5 (32c)

G 22 (z) =z -20 1.07e -06 + 0.149z -1 -0.524z -2 + 0.702z -3 -0.427z -4 + 0.0995z -5 1 -4.34z -1 + 7.74z -2 -7.08z -3 + 3.32z -4 -0.635z -5 (32d)

∞Fig. 1 .

 1 Fig. 1. Feedforward-feedback control schema used for robust anticipative control.

Fig. 2 .

 2 Fig. 2. Feedforward-feedback control schema using factorisation of nominal model and preview feedforward filter.

Fig. 3 .

 3 Fig. 3. Feedforward-feedback control schema used for robust anticipative control with sampling rate change and lowpass filter.

Fig. 4 .

 4 Fig. 4. High dynamic engine testbed schema.

Fig. 5 .

 5 Fig. 5. Bode diagram for 20 operating points (solid black line) and for the nominal model (dashed blue line).

Fig. 6 .

 6 Fig. 6. Engine speed control using feedback and preview feedforward.

Fig. 8 .

 8 Fig. 8. Current input using feedback and preview feedforward.

Fig. 9 .

 9 Fig. 9. Throttle input using feedback and preview feedforward.


Olivier Devillers
email: olivier.devillers@inria.fr

Charles Duménil
email: charles.dumenil@inria.fr

A Poisson sample of a smooth surface is a good sample

. Using this result, we prove that when points are Poisson distributed on a surface under the same hypothesis, whose expected number of vertices is λ, the expected size is O(λ log 2 λ).

Introduction

While the complexity of the Delaunay triangulation of n points is strictly controlled in two dimensions to be between n and 2n triangles (depending on the size of the convex hull) the gap between the lower and upper bound ranges from linear to quadratic in dimension 3. The worst case is obtained using points on the moment curve 1 and the best case by using the center of spheres defining a packing. 2To get a more precise result on the size of the 3D Delaunay triangulation, it is possible to make different kinds of hypotheses on the point set. A first possibility is to assume a random distribution in 3D and if the points are evenly distributed in a sphere [START_REF] Dwyer | The expected number of k-faces of a Voronoi diagram[END_REF], (resp. in a cube [START_REF] Bienkowski | Average case complexity of voronoi diagrams of n sites from the unit cube[END_REF]), Dwyer (resp. Bienkoswski et al.) proved that the expected size is Θ(n). But this hypothesis of random distribution is not relevant for all applications, for example when dealing with 3D reconstruction the Delaunay triangulation is an essential tool and it is much more natural to assume that the points are not distributed in space but on a surface [START_REF] Cazals | Delaunay triangulation based surface reconstruction[END_REF]. If the points are evenly distributed on the boundary of a polyhedron, the expected size was proved to be Θ(n) in the convex case [START_REF] Golin | On the average complexity of 3d-Voronoi diagrams of random points on convex polytopes[END_REF] and between Ω(n) and Õ(n) in the non convex case by Golin and Na [START_REF] Golin | The probabilistic complexity of the Voronoi diagram of points on a polyhedron[END_REF].

Instead of using probabilistic hypotheses one can assume that the points are a good sampling of the surface, namely an ( , η)-sample where any ball of radius centered on the surface contains at least one and at most η points of the point-set. Under such hypothesis Attali and Boissonnat proved that the complexity of the Delaunay triangulation of a polyhedron is linear [START_REF] Attali | A linear bound on the complexity of the Delaunay triangulation of points on polyhedral surfaces[END_REF]. Attali, Boissonnat, and Lieutier extend this result to smooth surfaces verifying some genericity hypotheses with an upper bound of O(n log n) [START_REF] Attali | Complexity of the Delaunay triangulation of points on surfaces: The smooth case[END_REF]. The genericity hypothesis is crucial since Erickson proved that there exists good sample of a cylinder with a triangulation of size Ω(n √ n) [START_REF] Erickson | Dense point sets have sparse Delaunay triangulations or ". . . but not too nasty[END_REF]. In the example by Erickson the point set is placed in a very special position on an helix, nevertheless, even with an unstructured point set it is possible to reach a supra-linear triangulation since Erickson, Devillers, and Goaoc proved that the triangulation of points evenly distributed on a cylinder has expected size Θ(n log n) [START_REF] Devillers | Empty-ellipse graphs[END_REF].

40:2

Poisson sample is good Contribution

In this paper we prove that a Poisson sample of parameter λ on a smooth surface of finite area is an (ε, η)-sample for ε = 3 » log λ λ and η = 1000 log λ with high probability. Using the result of Attali, Boissonnat, and Lieutier, it yields that the complexity of the Delaunay triangulation of a Poisson sample of a generic surface is O(λ log 2 λ) losing an extra logarithmic factor with respect to the case of good sampling (see Section 3).

Notation, definitions, previous results

We consider a surface Σ embedded in R 3 , compact, smooth, oriented and without boundary. At a point p ∈ Σ, for a given orientation, we denote by κ 1 (p) and κ 2 (p) the principal curvatures at p with κ 1 (p) > κ 2 (p). We assume that the curvature is bounded and define

κ sup = sup p∈Σ max(|κ 1 (p)|, |κ 2 (p)|).
We denote by σ(p, R) the sphere of center p and radius R. We denote by B(σ) the closed ball whose boundary is the sphere σ, by E the interior of a set E and, for p ∈ Σ, by D(p, R) the intersection between Σ and the B(σ(p, R)). Abusively we call D(p, R) a disk. For a discrete set X, we denote (X) the cardinality of X. If X is a set of points, Del (X) denotes the Delaunay triangulation of X. In the 3D case, (Del (X)) is the sum of the number of tetrahedra, triangles, edges and vertices belonging to the Delaunay triangulation.

Without loss of generality, we assume that Area(Σ) = 1 and consider that the set of points X is a Poisson point process with parameter λ > 0 over Σ.

We recall classical properties of a Poisson sample:

Observation 2.1. For two regions R and R of Σ,

P [ (X ∩ R) = k] = (λ Area(R)) k k! e -λ Area(R) , E [ (X ∩ R)] = λ Area (R), R ∩ R = ∅ ⇒ (X ∩ R
) and (X ∩ R ) are independent random variables.

In particular, we have P [ (X ∩ R) = 0] = e -λ Area(R) and E [ (X)] = λ. We consider the same definition of genericity as Attali, Boissonnat and Lieutier, roughly: the set of points where one of the principal curvatures is locally maximal is a finite set of curves whose total length is bounded and, the number of contacts of any medial ball with the surface is finite.

Then we define what is a good-sampling of a surface and precise the result by Attali, Boissonnat and Lieutier.

Definition 2.2 (Good sample).

A point-set on a surface is an (ε, η)-sample if any ball of radius ε centered on the surface contains at least one and at most η points of the sample.

Theorem 2.3 ([2]). The 3D Delaunay triangulation of an (ε, η)-sample of a generic smooth surface has complexity

O Ä η 2 ε 2 log 1 ε ä .
While the result of Attali et al. provides a bound O(N ln N ) on complexity of the Delaunay triangulation of an (ε, η)-sample of N points and a constant η, by looking more carefully at the result [2, Eq.( 14)], we notice that the actual complexity can be expressed by C( η ε ) 2 log(ε -1 ) for C being a constant of the surface.

3

Is a random sample a good sample?

In a Poisson sampling of parameter λ on the surface, a disk of radius ε = 1 √ λ is expected to contain π points, but with constant probability it can be empty or contains more than η points. Thus with high probability there will be such disks even if their number is limited. Thus such a sample is likely not to be a good sample with ε 2 = 1 λ and η constant. Nevertheless, it is possible to not consider η as a constant, namely, we take η = Θ(log(λ)). In a first Lemma, we bound the area of D(p, R), for any p ∈ Σ and R > 0 sufficiently small. Lemma 3.1. Let Σ be a smooth surface of curvature bounded by κ sup , and consider p ∈ Σ and R > 0 smaller than 1 κsup . The area of D(p, R) is greater than 3 4 πR 2 . Proof. The bound is obtained by considering the fact that the surface must stay in between the two tangent spheres of curvature κ sup tangent to the surface at p. The tangent disk at p of radius

√ 3 2 R > √ 3 2 1
κsup is included in the projection of D(p, r) on the tangent plane and thus has a smaller area than D(p, R). Lemma 3.2. Let Σ be a C 3 surface of curvature bounded by κ sup . For R small enough, Area(D(p, R)) < 5 4 πR 2 . Proof. Let z = f (x, y) := 1 2 κ 1 x 2 + 1 2 κ 2 y 2 + O(x 3 + y 3 ) be the Monge of Σ patch [START_REF] Peter W Hallinan | Twoand three-dimensional patterns of the face[END_REF] at a point p. We denote by dσ an element of surface and by A(p, R) the projection of D(p, R) on the xy-plane. Since on D(p, R) the slope of the normal to Σ is bounded, we have:

Area (D(p, R)) = D(p,R) dσ = A(p,R) 1 + ( ∂f ∂x (x, y)) 2 + ( ∂f ∂y (x, y)) 2 dxdy
That is smaller than

Area(D(p, R)) ≤ 2π θ=0 R r=0 r » 1 + 2(rκ sup ) 2 drdθ = π 3
(2(Rκ sup ) 2 + 1) Noticing that (a + 1)

3 2 -1 = a a+ √ a+1+2 √ a+1+1 ≤ 15
8 a for a < 1, we can conclude that for any R small enough,

Area(D(p, R)) ≤ π 3 15 4 (Rκ sup ) 2 κ 2 sup = 5 4 πR 2 . Lemma 3.3. Let Σ be a C 3 surface with Area(Σ) = 1. Let M R be a maximal set of k R disjoint disks D(p i , R) on Σ. If R is small enough then k R ≤ 4 3πR 2 .
Proof. By Lemma 3.1, for R small enough, we have D(p, R) ≥ 3 4 πR 2 . Thus: does not contain any point of X is O(λ -1 ).

k R • 3 4 πR 2 ≤ i=k R i=1 Area (D(p i , R)) ≤ Area (Σ) =
Proof. We prove that a Poisson sample has no empty disk of radius 3 » log λ λ with probability O(λ -1 ). In a first part we use a packing argument. On the one hand, for any ε > 0 small enough and given a maximal set M ε/3 and any point p ∈ Σ, the disk D(p, ε) contains entirely one of the disks D(p i , ε

3 ) belonging to M ε/3 . Indeed, by maximality of M ε/3 , the disk D(p, ε/3) intersects a disk of M ε/3 whose diameter is 2ε/3 so D(p, ε) contains it entirely. On the other hand, remember from Lemma 3.1 that if ε is small enough then Area(D(p, ε))≥ 3 4 πε 2 . Then we can bound the probability of existence of an empty disk for ε small enough:

Figure 1

 1 Figure 1 Illustration of the proof of Lemma 3.1 for the 2D case.

3 Figure 2

 32 Figure 2 A disk of radius ε always contains a disk of a maximal set of disks of radius ε 3 ,

  Let X be a Poisson sample of parameter λ distributed on a C 3 smooth closed surface Σ of area 1. If λ is large enough, the probability that there exists p ∈ Σ such that D p, 3 »

	1,
	and we can deduce the following bound: k R ≤ 4 3πR 2 .
	Lemma 3.4. log λ
	λ

The moment curve is parameterized by (t, t

, t

). When computing the Delaunay triangulation of points on this curve, any pair of points define a Delaunay edge.[START_REF] Attali | Complexity of the Delaunay triangulation of points on surfaces: The smooth case[END_REF] The kissing number in 3D is 12, thus in such a point set, the number of edges is almost 6n.

= 12 πε 2 e -λ 1πε 2 12 .

By taking ε = 3 » log λ λ we get:

We have proved that when a Poisson sample is distributed on a surface, the points sufficiently cover the surface, i.e., there is no large empty disk on the surface with high probability. Now we have to verify the other property of a good sample, namely, a Poisson sample does not create large concentration of points in a small area. Lemma 3.5. Let X be a Poisson sample of parameter λ distributed on a C 3 closed surface of area 1. If λ is large enough, the probability that there exists p ∈ Σ such that D(p, 3

Proof. Consider an M ε maximal set, we can notice that for any p ∈ Σ, the disk

Then we can bound the probability of existence of a disk containing more than η points:

We use a Chernoff inequality [START_REF] Mitzenmacher | Probability and computing: Randomized algorithms and probabilistic analysis[END_REF]

From Lemmas 3.1 and 3.2, we have that: 27 4 πε 2 ≤ Area(D(c, 3ε)) ≤ 45 4 πε 2 for ε small enough. Consequently we can say that the expected number of points v 0 in D(c, 3ε) verifies