Julien Vincent
email: julien.vincent@onera.fr

Pierre Borderies

Vincent Gobin

Jean-René Poirier

Modelling the Effects of Realistic Environments in the Nearfield and Farfield of Low Frequency Antennas with 3D FDTD Method

Keywords: FDTD, low frequency antennas, natural environments, Sommerfeld's problem, ground wave I

Ground wave radiation of Low Frequency

Antennas is well known for the canonical case of a flat, azimuthally homogeneous, soil. In this paper, more complex kinds of environment are considered in the surroundings of the antenna or in the far field. Both of them are treated by FDTD. Some realistic environments are studied in the near field zone including forest and slopes. An original methodology is proposed for the far field zone. It rests on the use of the direct numerical integration of Sommerfeld's integrals which is presented also.

INTRODUCTION

Ground wave radiation of Low Frequency Antennas is well known for the canonical case of a flat, azimuthally homogeneous soil. In this paper, more complex kinds of environment are considered in the surroundings of the antenna or in the far field. Both of them are treated by The Finite-Difference in Time Domain (FDTD). Some realistic environments are studied in the near field zone including forest and slopes. An original methodology is proposed for the far field zone. It rests on the use of the direct numerical integration of Sommerfeld integrals which is presented also. FDTD method is a well-known method used for computational electromagnetics [START_REF] Taflove | The Finite-Difference Time Domain Method[END_REF]. For the aim of this paper, the main advantage of using this technique is the ability to treat materials with any electromagnetic properties and any shape. In a previous paper [2], 2D FDTD is applied to obtain the electric field in the direction of propagation.

In section II, low frequency antennas are placed in a forest and on the top of a hill, the Earth's soil being considered as an infinite homogeneous half-plane elsewhere. In the last part of this paper, the method for computing the effect of discontinuities are computed far from the antennas using numerical integration of Sommerfeld's problem.

II. NEARFIELD MODELING

A. Principles

An inverted L antenna is computed at 100 kHz, geometry of the system is shown on Fig. 1. The antenna is made up off two parts, a radiating wire modelled with Holland formalism [START_REF] Holland | Finite-Difference Analysis of EMP Coupling to Thin Struts and Wires[END_REF] and a square shaped metallic ground. At this frequency range, caution should be taken in extending infinitely the soil so UPML [START_REF] Gedney | An Anisotropic Perfectly Matched Layer Absorbing Medium for the Truncation of FDTD Lattices[END_REF] are the boundary conditions applied here.

B. Antenna in a large forest and medium wet soil

With FDTD method, a forest can be represented by a dielectric slab [START_REF] Tewari | Evaluation of Relative Permittivity and Conductivity of Forest Slab from Experimental Measured Data on Lateral Wave Attenuation Constant[END_REF] with correspondent dielectric constants ε rf and σ f . In our case, we chose the medium wet case (ε rf =1.065; σ f =10 -3 S.m -1 ) applied for a 25 m high and square shaped (width= 12 km) forest [START_REF] Tewari | Evaluation of Relative Permittivity and Conductivity of Forest Slab from Experimental Measured Data on Lateral Wave Attenuation Constant[END_REF]. Fig. 2 shows the results of an inverted L antennas placed at the center of the forest at the frequency f =100 kHz. The power is sampled at 9 km around the antenna, 1 m over the ground and then normalized by the input power antennas. One can see that the shape of the radiation pattern is not affected by the introduction of the forest but that significant losses occur.

C. Antenna on the top of a hill and dry soil

An inverted L antenna at the frequency f =100 kHz is placed on the top of a hill made up off 4 square shaped dielectric blocs. All of them are 25 m high, their widths are calculated so as to create an eleven percent average slope and their dielectric constants are the same as the soil. Fig. 3 shows the results of the computation, the power is sampled at 9 km around the antenna, 1 m over the ground and then normalized by the input power antennas. The shape of the radiation pattern is not modified but in this case the radiated power is enhanced.

III. FAR FIELD MODELING

A. Sommerfeld's problem Principles

In 1909 [START_REF] Sommerfeld | Propagation of Waves in Wireless Telegraphy[END_REF], Sommerfeld calculated the radiation of an infinitesimal vertical electric current element I = I 0 dlz situated at a height h above a homogeneous dielectric half-space. The XOY plane separates the two media, the medium 1 is free space in the upper half-space with constants μ 0 , ε 0 and soil is the medium 2 in the lower half-space with constants μ 1 = μ 0 , ε 1 and σ 1 . The problem is presented on Fig. 4. Using Maxwell's equation, one can write the relations between Hertz's potential vector Π and the electromagnetic field E and H. Symmetries and translations show that only the vertical z-component of Π is not equal to zero and depends on the distance r and the height z of the observation point M. This leads to calculate the components E ρ , E z and H ϕ . The problem is separated into two parts: the source (Π 0 ) and the boundary conditions between medium 1 and medium 2 (Π 1 and Π 2 ). The potential vector Π 0 created by an electric dipole in air equals to a spherical wave with the magnitude of (μ 0 I 0 dl/4π). It has been shown by Sommerfeld [6] that the spherical wave expanding from a current element I 0 l can be represented by a spectrum of cylindrical waves. Each of which travels as J 0 (λr) in the radial r-direction and as exp[-|z-h|(λ 2 -k 0 2 ) 1/2 ] in the z-direction with a magnitude of (μ 0 I 0 dl/4π)λdλ/(λ 2 -k 0 2 ) 1/2 . Here, J 0 refers to the Bessel function of the first kind and k 0 is the free space wavenumber given by ω(μ 0 ε 0 ) 1/2 . Potential Vectors Π 1 and Π can be written with the same expansion in (1).

2

Here, k 2 is the soil wavenumber given by ω(μ 1 ε 1 ) 1/2 (1+σ 1 /jωε 1 ) 1/2 , R = (r 2 +(z-h) 2 ) 1/2 is the distance between the dipole and the observation point, R λ and T λ are coefficients which depend on integration variable λ and calculated with the boundary conditions at height z=0 in (2). At this height, components E ρ and H ϕ between the two media have to be continuous. This leads to a two equations system (3) with two unknowns R λ and T λ .

The values of R λ and T λ are finally obtained in [START_REF] Gedney | An Anisotropic Perfectly Matched Layer Absorbing Medium for the Truncation of FDTD Lattices[END_REF]. Their expressions can be interpreted as a reflection coefficient R λ and a transmission coefficient T λ .

Using these formulas, one obtains in (5) all components of the electromagnetic field E and H in both media [START_REF] Wu | Radiowave Propagation Over Ground[END_REF].

In medium 1, decomposition leads to three other equations in [START_REF] Sommerfeld | Propagation of Waves in Wireless Telegraphy[END_REF] which can be physically interpreted.

Here, R' = (r 2 +(z+h) 2 ) 1/2 is the distance between the image element at distance h below the plane and the observation point on Fig. 5. In [START_REF] Sommerfeld | Propagation of Waves in Wireless Telegraphy[END_REF], the components are expressed as the sum of three terms. The first represents a wave travelling directly from the current element I 0 dl, the second is from a negative image -I 0 dl and the third is called a surface wave term.

B. Numerical integration

The integrands in all components in both media are made up of three parts. The first element is the quotient whose denominator is

(k 0 2 (λ 2 -k 2 2 ) 1/2 +k 2 2 (λ 2 -k 0 2 ) 1/2
). An oscillating term is created by spherical Bessel functions J 0 and J 1 and finally there is an exponential attenuation factor depending on the height of the dipole and the observation point.

Here, integral from z-directed electric field in ( 7) is studied but the same approach is correct for the other components of the electromagnetic field.

The denominator does not have special values because integration path of λ is the real axis. When the real and imaginary part of the integrand is plot for a distance r, one can observe four principle parts. The first one has low oscillations, the second is the most important because it has the most data. In fact, it is where the real part of This analyze shows that different integration methods can be used for each part of the integrand. Third degree Newton-Cotes method is chosen for part one and three integrations and Gauss-Legendre quadrature with more than ten points is chosen for the integration of part two.

An adaptative integration algorithm was realized and tested for vertical component E z over a medium wet soil, along the direction of propagation. Fig. 6. shows the comparison between three different methods, the numerical integration, the FDTD and integral equations computed with an industrial software FEKO. In the following, the integration algorithm will be named NISP which stand for Numerical Integration of Sommerfeld's Problem.

C. NISP -3D FDTD hybridisation

The main of this part is to compute the effect of a discontinuity in/on the soil with NISP -3D FDTD hybridization. An object is placed in the far field of an infinitesimal vertical electric dipole. Incident electric E inc and incident magnetic H inc fields are computed around a cube with numerical integration algorithm and added as an incident wave source in the 3D FDTD. The situation with a homogeneous flat ground is shown on Fig 7.

First the Total Field / Scattered-Field technique is used with an object in the air to observe the scattering and validate the method and then a homogenous flat ground will be introduced. This part will be detailed in the presentation.

IV. CONCLUSIONS

The effects of natural environments in the near field and in the far field of antennas have been studied. First, FDTD is used to find a configuration where the system is more efficient and hybridization is realized with numerical integration in order to obtain the scattering of objects far from the source. The common use of these techniques will be helpful to better design LF antennas on irregular ground with emphasis on ground wave.

Fig. 1 .

 1 Fig. 1. Geometry of the computed antennas The Earth's soil is characterized with the dielectric constant ε rs and the conductivity σ s , which depend of the humidity, see Tab. 1, and can be found in the literature.

Fig. 2 .

 2 Fig. 2. Simulations of an inverted-L antenna with and without the forest on a medium wet soil.

Fig. 3 .

 3 Fig. 3. Simulations of an inverted-L antenna with and without the hill on a dry soil.

Fig. 4 .

 4 Fig. 4. Vertical z-directed electric current element I 0 dl at height h above planar lossy gruound with electrical parameters (μ 1 , ε 1 ,σ 1 ).

Fig. 5 .

 5 Fig. 5. Vertical z-directed current element source at height h and image at height -h above perfectly conducting plane.

  appears (λ ≈ k 0 ). The third part is the effect of the Bessel Function because both real and imaginary values of integrand are oscillating. The last part is where the exponential is most important and it is negligible despite the oscillating behavior.

Fig. 6 .

 6 Fig. 6. Solution comparison of vertical electric component in Sommerfeld's problem.

Fig. 7 .

 7 Fig. 7. NISP -3D FDTD Hybridation. An infinitesimal vertical electric dipole creates electric and magnetic fields (E inc , H inc ).

TABLE I . DIECLECTRIC CONSTANT AND CONDUCTIVITY FOR DIFFERENCT SOILS Earth's soils ε rs σ s (S.m -1 )

 I

	Dry	5	10 -4
	Medium wet	15	10 -3
	Wet	30	10 -2


Emmanuel Anger

Monika Srebro

Nicolas Vanthuyne

Christian Roussel

Loïc Toupet

Jochen Autschbach

Régis Réau

Jeanne Crassous

Loïctoupet

Helicene-grafted vinyl-and carbene-osmium complexes: an example of acid-base chiroptical switching

published or not. The documents may come

Helicene derivatives have recently shown potential as molecular materials due to their inherent chirality, large magnitude chiroptical properties and -conjugated electronic structure. 1 For example, they have revealed fascinating optical, redox or chemical-stimuli induced chiroptical switching. 2,[START_REF] Anger | Helicene based chiroptical switches[END_REF] These compounds may find applications invarious domains such as in optical displays, telecommunications, or molecular electronics.The design of new types of functional chiral molecular switches is needed and in this context, the switching of chiral metal-based systems has great potential, 2c due to the possibility of reversibly modifying the molecule by changing the metal center's properties (geometry, redox state, ...). 2c,d In this communication, we describe the synthesis of an unprecedented enantiopure carbeneosmium helicene complex from an enantiopure helicene-grafted vinylosmium derivative. The potential use of these chiral organometallic systems as acid-base chiroptical switches is examined. The electronic interaction between the vinyl-osmium or carbene-osmium units and thecarbo [6]helicenecore is investigated by DFT calculations.

Enantiopure vinyl-osmium-helicene complexes M-and P-2 (Scheme 1) were prepared in good yield (65-70%) using a known procedure, 4 i.e. by hydro-osmylation of enantiopureM-(-)-and P-(+)-1 with 1 equivalent of OsHCl(P i Pr3)2 in CH2Cl2 at room temperature (r.t.). The NMR data fully support the proposed structure in accordance with previously reported analogues. 4 For instance, the 1 H NMR spectrum of 2 displays the fingerprint for trans-ethylene metal-substituted moieties with a doublet of triplets at 4.80 ppm ( 3 JH-H= 13.4 Hz, [START_REF] Anger | Helicene based chiroptical switches[END_REF] JP-H= 1.8 Hz) corresponding to the proton in the position of the osmium center.Furthermore, according to 31 P NMR spectroscopy the two P i Pr3 ligands are non-equivalent (23.96 and 23.92 ppm), indicating slightlydifferent chemical environments for the two fragments. The Xray structure of Scheme 1.Synthesis of enantiopurehelicene-grafted vinyl-osmium complex P-2 and osmium-carbeneP-3fromenantiopureP-1.i) HOs(CO)Cl(P i Pr3)2, CH2Cl2, r.t., Ar, 3 days, 70%.ii) HCl, toluene, 10 min., 86%.iii) NEt3, CH2Cl2.

complex 2 (Figure 1) reveals a square pyramidal pentacoordinated osmium center in which the two trans P i Pr3 ligands, Cl and CO form the basisofthe pyramidand the vinyl-helicene ligand is at the apical position. The Os-C1-C2-C3 bond lengths (1.987, 1.343 and 1.453 Å) are indicative of a trans vinyl-Os moiety that is almost coplanar with the helicene part (Os-C1-C2-C3, 176.53°), allowing an efficient interaction of the helicene ligand with theOscenter. The helical angle (angle between the terminal rings of the helix, 50.3°) is similar to reported values for other metal-helicene compounds. 3a Overall, complex 2 is isostructural to its Ru analogue 3a,4b but osmium complexes are usually less easy to handle due to a lower stability.

Complexes containing Os=C(H)R units have been prepared. 4a-e Such carbene complexes (also named alkylidene-osmium complexes) are an important class of compounds but to our knowledge their electronic and chiroptical properties have not yet been studied.The first enantiopure helicene-based carbene-osmium complexes M-and P-3 (Scheme 1) were straightforwardly prepared from vinyl-osmium complexes M-and P-2 by a simple addition of slight excess ofHCl in toluene at r.t. 4a and characterized by 1 H, 31 P{ 1 H} and 13 C{ 1 H} NMR, along with IR, UV-vis and ECD spectroscopies. For example, the 1 H NMR spectrum of monocarbene 3 in CD2Cl2 showed the presence of a triplet at 18.1 ppm ( 3 JHH= 6.6 Hz) which is typical of the CH=Os and two doublets of doublets at 2.85 and 3.13 ppm corresponding to the methylene group. In the 13 C{ 1 H} NMR spectrum, the sp 2 carbon appears as a doublet of doublet ( 2 JCP = 8,5 Hz) at 296.7 ppm. Furthermore, two singlets at 12.89 and 12.47 ppm are found in the 31 P{ 1 H} NMR spectrum. These are typical values for aryl-methylene-carbene-osmium species 4a and correspond to the two P atoms in different chemical environments. Note that in such complexes the osmium centers are octahedral with a coordination sphere completed by a chloride ligand (see optimized structure in Figure 2). The Os-C1-C2-C3 dihedral angle in the calculated structure of 3is 159.35°,which indicates a rotation of the helicene moiety around C1-C2 axis as compared to 2. See Supplementary Information, SI, for Cartesian coordinates of both calculated structures 2 and 3.

Vinyl-osmium helicene complex 2 displays two strong high-energy absorption bands (258 and 314 nm,  ~ 70•10 3 M -1 cm -1 and 45•10 3 M - 1 cm -1 , respectively, see SI) that originate from the -system of the helicene ligand, along with abroadlow-energy band at 405 nm of moderate intensity ( ~ 14•10 3 M -1 cm -1 ) that does not exist in ligand 1 and is due to efficient electronic interaction between the vinyl-osmium moiety and the helicenic-core (see below). Interestingly, the carbeneosmium complex 3 displays the same high-energy absorption bands as 2 for the * transition of the helicene-core while the low-energy band at 405 nm has disappeared and two new broad bands of low intensity are found at 500 nm ( ~ 5•10 3 M -1 cm -1 ) and 785 nm ( ~ 3•10 3 M -1 cm -1 ), see SI. In addition, enantiopure complexes M-andP-2display molar rotation (MR) values that are much higher than the correspondingligand1(   23 D  in °cm 2 /dmol: P-2: +18000 vs. P-1: +11030).The circular dichroism (CD) spectrum of the complex P-2 depicted in Figure 1 shows, aside from the regular strong bands present in the helicene core, new ECD-active bands between 380-480 nm as compared with the ethynyl-substituted precursor P-1.These new bands involve electronic interaction between the vinyl-metal-center and the helical -system as confirmed by density functional calculations performed at the BHLYP/SV(P) level (see SI for Computational Details). For example, calculated excitations n°3 (367 nm) and n°4 (354nm)manifest character of π to π* transitions within the helicene ligand (HOMO-LUMO: 65% and HOMO-LUMO+1: 60%, respectively; see Figure 2 +11030).The ECD spectrum of complex P-3shows a sharp negative band at 250 nm and a sharp positive one at 330 nm (Figure 1). Additionally, new weakly ECD-active bands at 500 and 785 nm (~ 3-14 M -1 cm -1 ) are also present in P-3. TDDFT calculations on the singlet osmium carbene complexreproduced well the high-energy strong positive ECD-active bands in P-3 originating from the purely helicenic transitions (excitation n°5, 309 nm, HOMO-LUMO+2: 50% and HOMO-1-LUMO: 35%, see Figure 2 and SI), calculated in the same energy range where the strongest rotatory strength excitation of it organic precursor 1 occurs. As depicted in Figure 2 and SI, for 3 these MOs are centered on the -helical core without visible involvement of the osmium d orbitals.Unlike complex2, it appears that in the carbene-osmium derivative 3 the electronic coupling between metal center and the helicoidal-system is brokendue to the presence of the CH2. ). Inserts: low-energy CD spectra of P-2,3.X-ray structure of complex 2 (one enantiomer shown). Despite of the correct description of the high-energy part of the CD spectrum of 3, the very weak low-energy ECD-active bands could not be reproduced by a variety of non spin-orbit calculations. Absorptions to triplet excited states are often observed in heavy metal complexes because spin-orbit coupling weakens the spin selection rules. [START_REF] Hsu | Selected examples of Os complexes with triplet state absorption[END_REF] Calculations did indeed show the presence of triplet excitations in the wavelength range of interest. However, the appearance of triplet transitions above 600 nm was dependent on various approximations 7 (see SI for more details). 8 The calculated MR values of 3 were found to be overestimated as compared to the experimental result (17557 with BHLYP and 14173 with LC-PBE0) and similar to the calculated MR of 1. The trend among the set of complexes is thus correctly reproduced by the calculations. Another interesting question pertaining the electronic structure of 3 is the presence or absence of hyperconjugation( conjugation between CH2 and C=Osor helicene)which may affect structure and bonding in metal carbene complexes 9 and leave signatures in UV-vis spectra. 10 The optimized BP/SV(P) geometry of 3 does however not exhibit structural features that usually would be identified with hyperconjugation. 9,4a The ability of the helicene-osmium system to behave as an acidbase switch was examined. The reversible formation of helicene-vinylosmium P-2 from a solution of osmium-carbene-helicene P-3 in CH2Cl2was performed by simply using NEt3 and the UV-vis and ECD spectra of P-2were almost recovered (see SI). Furthermore, the IR spectra of complexes 2 showed a characteristic CO stretching mode (CO) at~1895 cm -1 that disappears upon adding HCl in toluene, while a new bandgrew up at ~1932 cm -1 corresponding to 3. This shows astronger donating character of the helicene-vinyl ligand compared to the carbene.The reverse transformation to 2 could be conducted by using NEt3 as followed by their IR CO (Figure 3).This transformation is thought to involve a deprotonation of the acidic CH2 group by the base followed by a 1,3-shift leading to HCl abstraction. 4f As a consequence, this molecular system acts as an acid-base triggered chiroptical switch since the several ECD bands (330 nm, 800 nm, ...) are greatly modified upon 2 3 transformation. Although the system is not practically efficient due to the low stability of highly diluted solutions of carbeneosmium derivatives,  to our knowledge this is the first time that a carbene derivative is used as a chiroptical switch.

In conclusion, the first enantiopure carbene-osmium based helicene derivative (3) has been prepared from unprecedented vinyl-osmium helicene derivative (2). The examination of the reversible transformation from2 3enabled to set up a new kind of an acid-base chiroptical switch due to the modification of the electronic interaction between the -helical core and a grafted osmium moiety. This example illustrates the inspiration from simple and well-known organometallic chemistry to the design of new kinds of helicene-based molecular switching systems.  CO (3)

  and SI). Both excitations are at low energy / long wavelength, due to an enlargement of the helicene chromophore by conjugation through the vinyl bridge and osmium d orbital involvement. Yet, because of the chirality of the helicene moiety, the transitions carry significant rotatory strength and contribute predominantly to the positive CD intensity at long wavelengths. In addition, a metal-centered excitation (n°2, 373 nm, HOMO-4-LUMO+2: 78%) with moderate rotatory strength is present at lower energy. The resulting intense low-energy tail of the first positive ECD band is responsible for the increase in the molar rotation of the metalgrafted vinyl-helicene derivatives compared to the free ligands. Indeed, the calculated MR values of the P-isomers (2: 23915 with BHLYP and 17435 °cm 2 /dmol with LC-PBE0) reproduce well the strong increase relative to the helicene ligand (1: 17002 (BHLYP) and 13800 °cm 2 /dmol (LC-PBE0). See ref. 5 and SI for more details). The osmium-carbene-helicene derivative revealed   23 D  values close to those of the ligands (   23 D  in °cm 2 /dmol:P-3:+9600 vs. P-1:

Figure1.

  Figure1.CD spectra of P-1a (blue) vs.P-2 (red) andP-3 (black) in CH2Cl2 (C 1-5 10 -5

Figure 2 .

 2 Figure 2. a) Calculated CD spectra ofP-1(blue), P-2(red),andP-3 (black)BHLYP/SV(P). b)Isosurfaces (0.04 au) of frontier MOs of 2and 3.

Figure 3 .

 3 Figure 3. Reversibility test of the (carbene/vinyl)-osmium 2 3 transformation.IR spectra in CH2Cl2 solution in the 1800-200 cm -1 region (CO stretching band). Spectrum obtained after 3→2 transformation is indicated in dashed red lines.

Acknowledgements:

The authorsthank the Ministère de la Recherche et de l'Enseignement Supérieur, the CNRS, the Région Bretagne, Rennes Métropole, and the ANR (12-BS07-0004-METALHEL-01). The theoretical work has been supported by grants CHE 0952253 and1265833from the National Science Foundation (JA) and the Foundation for Polish ScienceHoming Plus programme co-financed by the European Regional Development Fund under the Operational ProgrammeInnovativeEconomy(MS).The Center for Computational Research (CCR, Buffalo)has provided computational resources.

† Electronic Supplementary Information (ESI) available: experimental procedures, X-ray data and theoretical calculations. See DOI: 10.1039/b000000x/ Complex 3 is stable for several months if kept under rigorously inert atmosphere. However, it was not practically feasible to perform several reversible cycles due to successive adding of reagents in the solution and therefore changes in concentrations. 1 Y. Shen, C. -F. Chen, Chem. Rev., 2012, 112, 1463.


Philippe Cuvillier

Arshia Cont

COHERENT TIME MODELING OF SEMI-MARKOV MODELS WITH APPLICATION TO REAL-TIME AUDIO-TO-SCORE ALIGNMENT

Keywords: Hidden Markov model, semi-Markov chains, alignment, score following

This paper proposes a novel insight to the problem of duration modeling for recognition setups where events are inferred from time-signals using a probabilistic framework. When a prior knowledge about the duration of events is available, Hidden Markov or Semi-Markov models allow the setting of individual duration distributions but give no clue about their choice. We propose two criteria of temporal coherency for such applications and prove they are fulfilled by statistical properties like infinite divisibility and log-concavity. We conclude by showing practical consequences of these properties in a real-time audio-to-score alignment experiment.

INTRODUCTION

Many natural phenomena exhibit a latent temporal structure and may be observed through a temporal signal, e.g. music, speech or text. They are often structured as time-contiguous events which generate specific observations. In music, basic events may be notes -pitched sounds -and silences.

To recognize the sequence of events that generates the observed signal, probabilistic models [START_REF] Murphy | Dynamic Bayesian Networks: Representation, Inference and Learning[END_REF] are relevant when statistical relationships between observation, events and any additional variables are known. In particular the Hidden Markov Models (HMM) [START_REF] Rabiner | A tutorial on hidden Markov models and selected applications in speech recognition[END_REF] assume that the signal is stationary on time-intervals and identify them with the occupancy of a hidden state. HMMs represent possible states as vertices of a graph whose transitions are possible evolutions of the state. Once this state-space is specified, the Bayesian inference can be readily computed to recognize the state-sequence.

Score alignment [START_REF] Cont | A coupled duration-focused architecture for real-time music-to-score alignment[END_REF] is a Music Information Retrieval (MIR) task consisting of synchronizing a musical performance with its symbolic score. Since ordering of events is known, recognition boils down to alignment. Among the numerous applications of HMM, music has an outstanding This work was partially funded by the French National Research Agency ANR-11-JS03-005-01. characteristics: a music score assigns to each event its nominal duration, i.e. a prior information on their likely duration.

A crucial and undermined question is about the modeling of the nominal duration. This investigation is built on the framework of hidden semi-Markov models (HSMM) as it provides explicit choice of the duration model. In section 2 we briefly introduce HSMM. This generalization of HMM involves many Bayesian priors whose tuning is a major issue. To this aim, most probabilistic models rely on learning with training datasets [START_REF] Rabiner | A tutorial on hidden Markov models and selected applications in speech recognition[END_REF][START_REF] Mitchell | Modeling duration in a hidden Markov model with the exponential family[END_REF][START_REF] Raphael | Automatic segmentation of acoustic musical signals using hidden Markov models[END_REF][START_REF] Orio | Score following using spectral analysis and hidden Markov models[END_REF]. This paper presents an alternative based on a theoretical study of the inference equations.

In section 3, we state our first coherency criterion with aggregates and show how it leads to the property of infinite divisibility [START_REF] Steutel | Infinite Divisibility of Probability Distributions on the Real Line[END_REF]. In section 4, we state our second criterion, leading to specific statistical properties like log-concavity [START_REF] Johnson | Univariate Discrete Distributions[END_REF]. Afterwards, we compare our prescriptions with the choices of duration models we have found in the MIR literature. In section 5, we make a comparative test to illustrate the practical benefits of taking into account these theoretical properties.

BACKGROUND & MOTIVATION

Background: semi-Markov models for alignment

Hidden semi-Markov models were introduced in [START_REF] Levinson | Continuously variable duration hidden Markov models for automatic speech recognition[END_REF] as a generalization of HMM that offer explicit duration probability densities. Both are defined with two stochastic processes [START_REF] Rabiner | A tutorial on hidden Markov models and selected applications in speech recognition[END_REF]. The process (S t ) t∈N * is a discrete-time homogeneous Markov chain on a finite state-space E def = {1, 2, . . . , J}. Since its realizations (s t ) t∈N * are not observable, they are called hidden states. The observation (o t ) t∈N * is considered as a realization of the second process (O t ) t∈N * . We denote Markov models were originally designed for recognition.

In alignment tasks the prior information of ordering is conveniently modeled by left-to-right topologies of transition probabilities p ij . We exclusively deals with the simplest topology, the linear semi-Markov chains: ∀i, j, p ij = δ i,i+1 .

Right-censored Forward inference for real-time estimation

Offline inference uses the Viterbi algorithm [START_REF] Rabiner | A tutorial on hidden Markov models and selected applications in speech recognition[END_REF] at final time T to decode the most likely state sequence ŝT 1 . But for realtime alignment several estimators compete. Systems like [START_REF] Orio | Score following using spectral analysis and hidden Markov models[END_REF][START_REF] Schwarz | Robust polyphonic midi score following with hidden Markov models[END_REF][START_REF] Montecchio | A discrete filter bank approach to audio to score matching for polyphonic music[END_REF] prefer using the Forward algorithm at each time t to estimate the most likely current state ŝt . We follow them and exclusively study this inference mode.

The right-censored Forward algorithm for semi-Markov states as proposed in [START_REF] Guédon | Estimating hidden semi-Markov chains from discrete sequences[END_REF] provides suitable inference mechanism for real-time applications. It suggests computing the probabilities of the events {S t = j} instead of the events {S t = j, S t-1 = j}. This detail turns out to be crucial for real-time applications since at the current time t one cannot know if the current event S t is over yet. The algorithm introduces the quantities F j (t)

def = P(S t = j | O t 1 = o t 1 ), F o j (t) def = P(S t+1 = j, S t = j | O t 1 = o t 1 ) and F i j (t) def = P(S t+1 = j, S t = j | O t 1 = o t 1 )
. To compute them, a key feature is the survivor distribution D j (u) def = P(S t+u t+2 = j | S t+1 = j, S t = j). It is induced by d j owing to the relationship D j (u) = v≥u d j (v). The Forward inference consists of the following recursions over time t

F j (t) = t u=1 b j (o t t-u+1 )F i j (t -u)D j (u) F o j (t) = t u=1 b j (o t t-u+1 )F i j (t -u)d j (u) (1)
while F i j (0) = π(j) and F i j (t) = i =j p ij F o i (t). The estimated current state is implemented as ŝt def = arg max j F j (t).

Motivation: modeling prior information of duration with HSMMs

Using semi-Markov models to decode event sequences requires a careful design of the Bayesian prior d j for each state.

Most systems tune these with statistical learning using solutions like the HSMM version of the Baum-Welch algorithm [START_REF] Guédon | Estimating hidden semi-Markov chains from discrete sequences[END_REF]. But estimating a distribution d j per state j would require huge training datasets. Consequently most systems learn the occupancy distributions with regularization techniques over parametric classes of pdfs [START_REF] Mitchell | Modeling duration in a hidden Markov model with the exponential family[END_REF]. This has led to ad hoc engineering choices to reduce the parameter space. This paper focuses on this tradition, and studies their mathematical behavior, application consequences to get insights on these choices.

Our study starts from an outstanding property of our MIR application: musical events are associated with a reference duration. Indeed a music score provides the prior tempo and prior durations for all notes. We denote this quantity the nominal duration l j . Although a few music alignment systems like [START_REF] Joder | An improved hierarchical approach for music-to-symbolic score alignment[END_REF] willingly discard this prior information, this work considers durations as an explicit element of modeling.

Accounting for nominal durations requires the following strong assumption: two events with identical nominal durations should get identical occupancy distributions. So the duration model consists of a set of durations L ⊂ R + and a duration-indexed family of pdfs (d l ) l∈L such that for all state j, l j ∈ L and d j = d lj . This framework sharpens the problematic: are there relevant choices of the mapping from nominal durations l to distributions d l ?

A review of literature shows that many modeling heuristics compete. Some of them are listed in section 4.5. Up to now no investigation has been done to compare such choices. The next two sections develop two criteria so as to justify or disqualify such heuristics. We introduce our first criterion of coherency and show that it leads to the notion of infinite divisibility and convolution semigroups. The motivation comes from aggregates of events. Figure 1 illustrates the concept with two toy sequences. The music score 1b informs of four events with equal nominal duration of 1 sec. In the music score 1a one silence of duration 2 replaces the two ones of duration 1. These two sequences would generate the same observations (O t ) since no physical signal could distinguish consecutive silences. The "one semi-Markov state per event" strategy would map 1a to the 3 states graphical model (A, B, C) and 1b to the 4 states one (A, B 1 , B 2 , C), so state B carries the same prior as the aggregate (B 1 , B 2 ). Coherency would ask that these two different graphical models infer the same values for (B 1 , B 2 ) and B respectively. To achieve this criterion, we study duration prior of linear aggregates. The duration spent in state j is a random variable L j whose law is d j . Intuitively the duration spent in two consecutive states i.e. the aggregate (j, j + 1) is the sum of their individual durations: L (j,j+1) = L j + L j+1 . As the Markovian assumption makes these two random variables independent, the law of an aggregate is the convolution of individual laws: d (j,j+1) = d j * d j+1 .

CRITERION 1: COHERENT TIME-AGGREGATES

Remark: the discrete convolution product * is defined as

[f * g](t) def = u∈Z f (u)g(t -u).
Actually, whether this intuition is valid depends on the inference mode. Simple algebraic computations prove it does not hold for the Viterbi algorithm whereas it does for the Forward one. Details are given in the following proposition. This proposition proves that the family (d l ) l∈L respects criterion 1 if

∀l 1 , l 2 ∈ L, l 1 + l 2 ∈ L (2) d l1+l2 = d l1 * d l2 (3)
It turns out that equation ( 2) is the exact definition of L being an additive subsemigroup of R + while equation ( 3) is the exact definition of (d l ) l∈L being a convolution semigroup [START_REF] Steutel | Infinite Divisibility of Probability Distributions on the Real Line[END_REF]. Now we study the question: do such families exist? The answer depends on the structure of L.

Case A. One could say that L = {l 0 , 2l 0 , 3l 0 , . . .} for some l 0 > 0. This base duration l 0 is called the tatum, or time quantum in the MIR literature [START_REF] Klapuri | Analysis of the meter of acoustic musical signals[END_REF]. To build the associated semigroup, one can choose d l0 as any valid pdf then compute its successive convolution powers.

Case B. One could say that L contains the subdivisions l 0 , l 0 /2, l 0 /3, . . . of some l 0 . Indeed music scores are written with rational subdivisions of a base duration and may virtually contain all rational values. This implies that

∀n ∈ N * , d l0 = d l0/n * d l0/n * . . . * d l0/n n times
. This latter property is the exact definition of d l0 being an infinitely divisible distribution [START_REF] Steutel | Infinite Divisibility of Probability Distributions on the Real Line[END_REF]. As this reference explains, such a distribution induces a convolution semigroup (d l ) l>0 on the "full" set L = R + . The discrete class of such distributions are known as compound Poisson distributions. Unfortunately none of them check the requirement d(0) = 0 to define a valid occupancy distributions so no valid family respects the exact condition [START_REF] Cont | A coupled duration-focused architecture for real-time music-to-score alignment[END_REF]. In practice we prescribe choosing a true convolution semigroup (d l ) l>0 and truncating the values at 0. This approximation of ( 3) is all the better as durations l are long since for all convolution semigroups,

d l (0) = (d 1 (0)) l → l→∞ 0.
Moreover, this problem of existence is due to time discretization. It could be solved by considering continuous-time HSMMs [START_REF] Levinson | Continuously variable duration hidden Markov models for automatic speech recognition[END_REF] but these models bring about implementation difficulties.

CRITERION 2: COHERENCY UNDER NON-DISCRIMINATIVE OBSERVATION

Non-discriminative observation

In this section we derive a second criterion by extrapolating the situation of "repeated events" that motivates previous section. What would happen if all events share the same observation probabilities? We call non-discriminative observation a model where b 1 = b 2 = . . . or equivalently b j ≡ 1. This assumption may hold in other realistic situations of Bayesian inference such as missing observations [START_REF] Yu | A hidden semi-Markov model with missing data and multiple observation sequences for mobility tracking[END_REF] and outliers. It also approximates observation probabilities b j being "too close" because of noisy training data for example. Consequently, we argue that the behavior of inference under non-discriminative observation assesses its robustness.

Inference in non-discriminative observation

We introduce our second criterion and show that it leads to new prescriptions inspired by precise notions of statistics. If the observation probabilities do not discriminate states, then the inference should respect the states ordering and their nominal duration as these are the only available information.

Coherency criterion 2. On a linear graphical model with non-discriminative observation, the inference successively decodes states 1, 2, . . . and assigns to each state j a duration which is equal to its nominal duration l j . This criterion is not validated by most existing systems and is worth a theoretical investigation. The Forward estimation respects criterion if and only if

∀j ∈ E, F j+1 (t) F j (t) ≤ 1 if t < l 1 + l 2 + . . . + l j > 1 else .
The assumption of non-discriminative observation leads to successive convolutions over Forward equations [START_REF] Murphy | Dynamic Bayesian Networks: Representation, Inference and Learning[END_REF]. A simple induction over states j proves that

F o j = d 1 * d 2 * . . . * d j F j = d 1 * d 2 * . . . * d j-1 * D j if j > 1 D 1 else
Now we are able to compare the Forward quantities between successive states F j (t), F j+1 (t). Since state 1 is the most likely one at first time step t = 1, we begin with a comparison between states 1 and 2.

The case of two-state chains

To achieve criterion 2, we successively present partial results on the evolution of the following quantity: F2(t) F1(t) = [d1 * D2](t) D1(t) . To begin with, the median gives a universal lower bound for the duration assigned to state 1 and this bound is tight.

Proposition 2. Let us denote m[d 1 ] def = max{t | D 1 (t) ≥ 1/2} the median of d 1 . Then for any distribution d 2 , t ≤ m[d 1 ] ⇒ F 1 (t) ≥ F 2 (t).
Reciprocally, there exists distributions

d 2 such that t > m[d 1 ] ⇒ F 1 (t) < F 2 (t).
Proof. The sufficient condition is straightforward. Since D 2 (t) ≤ 1 for all t,

t-1 u=1 d 1 (u)D 2 (t -u) ≤ t-1 u=1 d 1 (u) so F 2 (t) ≤ 1 -D 1 (t) and F 2 (t) -F 1 (t) ≤ 1 -2D 1 (t). Since D 1 is non-increasing, if t ≤ m[d 1 ] then D 1 (t) ≥ 1/2 and 1 -2D 1 (t) ≤ 0.
For the necessary condition, one may consider the trivial distribution

d 2 (t) = δ m[d1] (t).
Proposition 2 provides half of criterion 2: state 1 is more likely than state 2 while t < m[d 1 ]. But will state 2 be more likely afterwards? In the general case the answer is negative: there are distributions d 2 such that d 1 * D 2 (t) < D 1 (t) for all time t. Fortunately, we have found out that standard properties of statistics are relevant for this question.

Definition 1. A distribution d is Increasing Hazard Rate (IHR) if its hazard rate h(n) def = d(n) D(n) is non-decreasing.
Next proposition reveals the interest of IHR for our problem. Its proof is straightforward and is also a consequence of [START_REF] Lynch | Uniform stochastic orderings and total positivity[END_REF]Corollary 2.3].

Proposition 3. If d 1 is IHR then d1 * D2(t) D1(t)
is an increasing function of t for any distribution d 2 .

When this proposition holds, then by the monotone convergence theorem the limit f 1,2 def = lim t→∞ F2(t) F1(t) ∈ [0, +∞] exists. Two cases appear: If f 1,2 ≤ 1 then state 2 is never decoded and the criterion 2 is not fulfilled. Else there exists a unique time t 1 such that state 1 is more likely before and state 2 likewise after t 1 .

So now criterion 2 is almost reached. The remaining condition t 1 = l 1 can be numerically checked on the chosen couple (d 1 , d 2 ). We left a deeper investigation and only make two assertions. Firstly, f 1,2 = +∞ holds if d 1 has equal or faster tail decay than d 2 . Secondly, typical values of t 1 are located on the median m[d 1 ] or nearby for many couples (d 1 , d 2 ). In conclusion, we prescribe to calibrate the median of all distributions d l on l.

Generalization to N -state chains

The previous arguments cannot be generalized to more than two states without further assumptions. We present one sufficient condition which is informative although restrictive.

Definition 2. A discrete distribution d is log-concave if for all n in N * , d(n) 2 ≥ d(n -1) d(n + 1).
It is noteworthy that all log-concave distributions are IHR. Proof. Again we only sketch the proof. Firstly, for all j the quantities Fj+1(t) Fj (t) are non-decreasing. The proof is similar to proposition 3. Then these quantities all diverge to +∞ so state j +1 becomes more likely than state j after some time. Finally, they check the following property: F j+1 (t) < F j (t) implying F j+2 (t) < F j+1 (t). This ensures that each state j is decoded before state j +1.

Log-concavity plays an important role in many fields of statistics but has been scarcely studied on HSMMs. It has been highlighted in [START_REF] Tweed | Efficient hidden semi-Markov model inference for structured video sequences[END_REF] for improving computational efficiency of the Viterbi inference. Proposition 4 also prescribes this property for all d l as it provides theoretical coherency to the Forward inference.

Comparaison with pre-existing heuristics

Previous sections lead to prescriptions on the family (d l ) l∈L . We confront them with a list of duration models from the MIR literature. Refer to [START_REF] Johnson | Univariate Discrete Distributions[END_REF] for properties of the usual probability laws.

• [START_REF] Nakamura | Acoustic score following to musical performance with errors ad arbitrary repeats and skips for automatic accompaniment[END_REF] chooses geometric laws with mean fitted on l: d l ∼ G(1 -1/l). The family (d l ) l>1 respect criterion 2 but not criterion 1 as it is not a convolution semigroup;

• [START_REF] Cont | A coupled duration-focused architecture for real-time music-to-score alignment[END_REF] chooses exponential laws with mean fitted on l: d l ∼ E(1/l). The family (d l ) l>0 respect no criteria. It is not a convolution semigroup and the median of d l is not l but (l ln 2).

• [START_REF] Takeda | Rhythm and tempo analysis toward automatic music transcription[END_REF] chooses log-normal distributions with constant shape parameter σ > 0 and log-scale fitted on l: d l ∼ ln N (ln l, σ). This family is not a convolution semigroup and no d l is IHR.

• [START_REF] Joder | An improved hierarchical approach for music-to-symbolic score alignment[END_REF] uses normal distributions with mean equal to l and standard deviation proportional to l: d l ∼ N (l, l 2 σ 2 ) for some σ > 0.

[20] makes the same choice but sets the variance proportionally to l: d l ∼ N (l, lσ 2 ). Only the latter choice gives a convolution semigroup and it respects all prescriptions.

• [START_REF] Montecchio | A discrete filter bank approach to audio to score matching for polyphonic music[END_REF] uses negative binomial distributions with mean fitted on l: d l ∼ N B(l(1-p)/p, p) for some p ∈ (0, 1). This defines a convolution semigroup that do respect our prescriptions for long enough l. This might explain why the authors does observe their inference working well with repeated events.

Alternative approaches like [START_REF] Duan | Soundprism: An online system for score-informed source separation of music audio[END_REF][START_REF] Montecchio | A unified approach to real time audio-to-score and audio-to-audio alignment using sequential Monte Carlo inference techniques[END_REF] define the state-space E on the continuous real line and events as contiguous intervals. Inference is implemented by particle filtering methods. These models turn out to respect criterion 1. Indeed their occupancy distributions are implicitly defined as first-passage times of a diffusion process and such variables are always infinitely divisible [7, Section 7].

RESULTS & EXPERIMENTS

Accounting results of section 3 and 4, we recommend the following choice for duration model: the family of Poisson laws (d l ) l>0 ∼ (P o(l)) l>0 . This choice is the simplest compound Poisson semigroup and is optimal for all criteria exposed previously. Indeed, [START_REF] Watanabe | On the strong unimodality of Lévy processes[END_REF]Theorem 2] has proven that it is the only semigroup whose distributions are log-concave for all l > 0. Moreover the median of d l is close to l for all l.

To illustrate the benefits of our proposal, we expose comparative results of real-time audio-to-score alignment on an example that challenges state-of-the-art real-time systems. The example is a 1939 piano performance of Chopin's Mazurka Op. 17 No. 2 by A. Rubinstein from the Mazurka Project dataset [START_REF] Sapp | The Mazurka Project[END_REF]. For this experiment, we use exactly the same HSMM system proposed in [START_REF] Cont | A coupled duration-focused architecture for real-time music-to-score alignment[END_REF] (the Antescofo system which constitutes both practical and scientific state-of-the-art in this application domain). The default setting for d l consists of exponential laws and does not respect our prescriptions as discussed in section 4.5. In our proposal we set d l to Poisson laws but do not change anything else. Figure 2 compares the result of the alignment of this audio to its music score with these two settings1 . The default system gets lost after some seconds whereas the proposal manages the entire piece without losing track and with good precision. Similar enhancements have been observed on the whole mazurkas dataset but we leave deeper quantitative evaluation for further publications.

From this experiment, we conclude that the coherency of the duration model compensates the weakness of the observation model. Indeed usual performances of mazurkas make a deep use of the piano pedal that adds some pitches of a chord to the consecutive ones. These unexpected pitches "blur" probabilities between neighboring states. As a result the observation loses its the discriminative power and the inference rely more deeply on its duration model.

CONCLUSION & PERSPECTIVES

This papers introduces two criteria of time-coherent modeling in semi-Markov models for alignment. First condition is inference consistency with aggregates and subdvisions; second condition is estimation coherency under non-discriminative observation. We show that coherency is theoretically guaranteed if the chosen probability distributions has specific statistical properties and that respecting these prescriptions experimentally improves real-time alignments.

This short study calls for further theoretical and experimental developments. More necessary and sufficient conditions related to the criteria can be derived. The framework can be extended to other estimators such as Viterbi and Forwardbackward algorithms. Moreover the proposed prescriptions lead to constraints on the learning parameter space; adding these constraints in HSMM training algorithms would be an interesting issue. Finally, compound Poisson distributions have been extensively studied in statistics and econometrics. The consequences of their numerous properties for Bayesian inference would be worth deeper investigations.

=

  {0, 1, . . .} and N * def = {1, 2, . . .}. In probabilistic models the duration spent on a state j is a random variable L j . Its law is called the occupancy distribution d j (u) def = P(S t+u+1 = j, S t+u t+2 = j | S t+1 = j, S t = j). It is only defined for u ∈ N * so it must check d j (0) = 0. The explicit choice of d j sets HSMM apart from regular HMM modeling. For a Markovian state with self-transition p, d j implicitly follows a geometric law d j (u) = (1 -p)p u-1 . 978-1-4799-3694-6/14/$31.00 c 2014 IEEE Any semi-Markov chain consists of two additional choices per state j: the initial probability π(j) def = P(S 1 = j), and the transition probabilities p ij def = P(S t+1 = j | S t+1 = i, S t = i) with p ii = 0 for semi-Markov states. Hidden model paradigm describes how states (S t ) influence observations (O t ) using observation probabilities b j (o t+u t We make the usual Markovian assumption of conditional independence: ∀t, u, b j (o t+u t ) = u v=0 b j (o t+v ).

Fig. 1 :

 1 Fig. 1: Two equivalent event sequences and their graphical models. The silence B aggregates the short ones (B 1 , B 2 ).

Coherency criterion 1 .

 1 An aggregate of N states with nominal duration l 1 , l 2 , . . . and identical observation functions b 1 = b 2 = . . . = b induces the same inference quantities as a single state with duration l 1:N def = l 1 + . . . + l N and observation b.

Proposition 1 (

 1 Consistency with aggregates). For the Forward inference a linear aggregate (B 1 , B 2 , . . . , B N ) is equivalent to a single state B if and only if 1. the observation probabilities are identical: b B = b B1 = b B2 = . . . and check the Markovian assumption -see section 2.1, 2. the initial probabilities check π(B i ) = π(B) δ 1,i , 3. the occupancy distributions check d B = d B1 * d B2 * . . . * d B N .

Proposition 4 .

 4 If all states share the same occupancy distribution, i.e. d 1 = d 2 = . . . = d and if d is log-concave then all states are decoded in their ordering.

Fig. 2 :

 2 Fig. 2: Comparison of alignments. Left: d l are exponential laws. Right: d l are Poisson laws. Reference events are music beats. For each beat the error is the absolute time lapse between the estimate and true time of the beat. Alignment rate is the percentage of beats whose error is below 100 ms.

Files, video of the experiments and more details on the results are accessible on http://repmus.ircam.fr/mutant/mlsp14.


Z Wang
email: z.h.wang@hit.edu.cn

M Rodrigues
email: mickael.rodrigues@univ-lyon1.fr

D Theilliol
email: didier.theilliol@univ-lorraine.fr

Y Shen

Fault Estimation for A Class of Discrete-Time Nonlinear Systems ⋆

Keywords: Actuator fault estimation, nonlinear systems, Linear Parameter-Varying (LPV), observer design, Linear Matrix Inequality (LMI)

This paper deals with actuator fault estimation for a class of discrete-time nonlinear systems whose linear part is described by a Linear Parameter-Varying (LPV) form and nonlinear term is Lipschitz. In this paper, an augmented system is constructed by considering the fault as an auxiliary state vector. Then, a robust fault estimation observer is designed based on the augmented system. The proposed fault estimation observer is able to attenuate the effect of the fault variation and measurement noise and the observer design is formulated as a Linear Matrix Inequality (LMI) feasibility problem, which can be easily solved. Finally, a missile model is used to demonstrate the effectiveness of the proposed method.

INTRODUCTION

With the growing complexity of modern engineering systems and ever increasing demand for safety and reliability, fault diagnosis techniques have received great attention during the past decades. Generally, fault diagnosis approaches can be classified into three categories, i.e., modelbased methods, signal-based methods, and knowledgebased methods [START_REF] Dai | From model, signal to knowledge: A data-driven perspective of fault detection and diagnosis[END_REF]). Among these approaches, model-based fault diagnosis methods have been most intensively studied and various methods have been proposed in the literature, see e.g., [START_REF] Frank | Fault diagnosis in dynamic systems using analytical and knowledge based redundancy-A survey and some new results[END_REF], [START_REF] Chen | Robust Model-Based Fault Diagnosis for Dynamic Systems[END_REF], [START_REF] Isermann | Model-Based Fault-Detection and Diagnosis-Status and Applications[END_REF], [START_REF] Hwang | A survey of fault detection, isolation, and reconfiguration methods[END_REF], and the references therein.

As is known, most of the practical systems are nonlinear in nature. Therefore, there has been a lot of interest in fault diagnosis for nonlinear systems in the recent years. For example, [START_REF] Jiang | Fault accommodation for nonlinear dynamic systems[END_REF] and [START_REF] Zhang | Fast fault estimation and accommodation for dynamical systems[END_REF] studied fault diagnosis for Lipschitz nonlinear systems by using adaptive observer. In [START_REF] Polycarpou | Automated fault detection and accommodation: A learning systems approach[END_REF]; [START_REF] Polycarpou | Learning approach to nonlinear fault diagnosis: Detectability analysis[END_REF], the online approximators are used to develop fault diagnosis methods for nonlinear systems. In addition, Proportional Integral (PI) observer and descriptor observer are also be used to develop fault diagnosis methods for nonlinear systems in the recent years, see e.g., [START_REF] Zhang | Observer-based integrated robust fault estimation and accommodation design for discrete-time systemss[END_REF], [START_REF] Astorga-Zaragoza | Fault diagnosis for a class of descriptor linear parameter-varying systems[END_REF], [START_REF] Gao | Sensor fault reconstruction and sensor compensation for a class of nonlinear state-space systems via a descriptor system approach[END_REF] etc. However, most of these results are dedicated to a specific class of nonlinear systems, i.e. Lipschitz nonlinear systems. In fact, fault diagnosis in nonlinear systems is still a difficult problem because there does not exist a universal methodology that can be applied to all nonlinear systems. Since Linear parameter-varying (LPV) systems can be used to approximate nonlinear systems, the LPV representation provides a systematic and elegant way to deal with fault diagnosis for nonlinear systems. In recent years, a number of fault diagnosis methods have been proposed based on LPV models, e.g., [START_REF] Bokor | Detection filter design for LPV systems-A geometric approach[END_REF]; [START_REF] Alwi | Fault reconstruction using a LPV sliding mode observer for a class of LPV systems[END_REF]; [START_REF] Rodrigues | Sensor fault detection and isolation filter for polytopic LPV systems: A winding machine application[END_REF], etc. However, most of these methods are presented for continuous-time systems. In contrast to the continuous-time case, few results have been reported on fault diagnosis for discrete-time nonlinear systems. Nowadays, most of control systems are digitally implemented. Therefore, it is necessary to study the fault diagnosis design for discrete-time LPV systems.

On the other hand, when a nonlinear system with complex nonlinearity is considered, the corresponding LPV model will consists of a large number of linear models, which makes the fault diagnosis design complicated and also brings computation burden. Therefore, it is important to develop a more proper model to simplify the fault diagnosis design. In [START_REF] Boulkroune | H -/H ∞ fault detection filter for a class of nonlinear descriptor systems[END_REF], the authors considered robust fault detection for a class of nonlinear descriptor systems whose nonlinear term is globally Lipschitz and linear term is described by a LPV form. However, the methods in [START_REF] Boulkroune | H -/H ∞ fault detection filter for a class of nonlinear descriptor systems[END_REF] only address the fault detection problem but cannot be used to estimate the fault magnitudes. Here in our paper, the aim is not only to develop a fault detection method for discrete-time LPV systems with nonlinear terms but also to estimate actuator fault for such systems.

In this paper, a class of discrete-time nonlinear systems consists of a LPV term and Lipschitz nonlinearity are considered. The basic idea of this paper is to construct an augmented system by taking the fault as an auxiliary state vector. Then, a robust fault estimation observer is designed based on this augmented system. In this paper, the proposed fault estimation observer design is formulated as a set of Linear Matrix Inequalities (LMIs), which can be efficiently solved. The main contribution of this paper consists in the following aspects. First, a new fault estimation method is proposed for a class of discrete nonlinear systems which contains a LPV term and a Lipschitz nonlinear term. To the best of our knowledge, no work has been done on this subject. Moreover, fault variation and measurement noise are considered in the robust fault estimation observer design, which makes the proposed approach practical for real systems.

PROBLEM FORMULATION

Consider the following discrete-time nonlinear system

         x(k + 1) = h i=1 ρ i (θ(k)) A i x(k) + B i u(k) + F i f (k) +GΦ(θ(k), x(k)) y(k) = Cx(k) + Hv(k) (1) where x(k) ∈ R n is the state vector, u(k) ∈ R p is the input vector, f (k) ∈ R q is the fault vector, y(k) ∈ R m denotes the measured output and v(k) ∈ R m represents the measurement noise. Φ(θ(k), x(k)) ∈ R s is a known nonlinear function. A i ∈ R n×n , B i ∈ R n×p , F i ∈ R n×q , G ∈ R n×s , C ∈ R m×n
and H ∈ R m×n are known constant matrices. θ(k) ∈ R l is a bounded time-varying parameter vector measurable online and not affected by fault, which is a general assumption in the LPV framework, see e.g. [START_REF] Boulkroune | H -/H ∞ fault detection filter for a class of nonlinear descriptor systems[END_REF], [START_REF] Alwi | Fault reconstruction using a LPV sliding mode observer for a class of LPV systems[END_REF] and [START_REF] Rodrigues | Sensor fault detection and isolation filter for polytopic LPV systems: A winding machine application[END_REF]. ρ i (θ(k)), i = 1, 2, . . . , h are the weighting functions which satisfy the following properties for all k:

h i=1 ρ i (θ(k)) = 1, 0 ≤ ρ i (θ(k)) ≤ 1, i = 1, 2, . . . , h (2)
Herein, h is the total number of weighting functions.

Without loss of generality, it is assumed that rank

(C) = m, rank(F i ) = q, q ≤ m. Moreover, the nonlinear function Φ(θ(k), x(k)) is assumed to satisfy the Lipschitz condition, i.e. ||Φ(θ(k), x 1 (k)) -Φ(θ(k), x 2 (k))|| ≤ L g ||x 1 (k) -x 2 (k)||(3)
where L g > 0 denotes the Lipschitz constant.

The aim of this paper is to design an observer to obtain a fault estimation for f (k).

To this end, the LPV system (1) is firstly formulated as an augmented system representation by letting f (k) as a part of the augmented state vector, i.e.

x(k) = x(k) f (k) (4)
Then, an augmented system is constructed as follows

         x(k + 1) = h i=1 ρ i (θ(k)) Āi x(k) + Bi u(k) + D f (k) + ḠΦ(θ(k), C x x(k)) y(k) = C x(k) + Hv(k) (5) where Āi = A i F i 0 I q , Bi = B i 0 , D = 0 I q (6) Ḡ = G 0 , C x = [ I n 0 ] , C = [ C 0 ] (7) and f (k) = f (k + 1) -f (k) (8)
is the fault variation between two consecutive instants k and k + 1.

According to the definition of x(k) in ( 4), the actuator fault f (k) can be estimated if the estimation of the augmented state x(k) is obtained. Based on this idea, this paper proposes the following fault estimation observer

x(k + 1) = h i=1 ρ i (θ(k)) Āi x(k) + Bi u(k) +L i (y(k) -C x(k)) + ḠΦ(θ(k), C x x(k)) (9)
where x(k) ∈ R (n+q) is the estimate of the augmented state x(k). L i ∈ R (n+q)×m , i = 1, 2, . . . , h are matrices to be synthesized.

After the augmented state estimation x(k) is obtained, the fault estimation can be determined as follows

f (k) = C f x(k) ( 10
)
where

C f = [ 0 I q ] (11)
Now, the problem is to design observer (9) such that the fault estimation f (k) is robust to the fault variation f (k) and measurement noise v(k). This is an H ∞ design problem, which can be formulated as follows.

H ∞ fault estimation problem: Consider (5) and the fault estimation observer (9), the objective is to design that matrices

L i , i = 1, 2, . . . , h such that the fault estimation error e f (k) = f (k) -f (k) robustly converges torwards zero, i.e. e f 2 ≤ γ 2 1 f 2 2 + γ 2 2 v 2 2 + V (0) (12)
where • 2 represents the L 2 norm, γ 1 and γ 2 are the disturbance attenuation levels: γ 1 is related to the robustness against the fault variation f (k), γ 2 represents the attenuation level of the measurement noise and V (0) is a quadratic function which will be given latter.

FAULT ESTIMATION OBSERVER DESIGN

In this section, a solution to the H ∞ fault estimation problem is proposed.

First, a sufficient condition of the criterion ( 12) is given in the following Theorem. Theorem 1. The criterion ( 12) is fulfilled if there exists a Lyapunov function V (k) such that the following inequality holds for all k:

∆V (k)+e T f (k)e f (k)-γ 2 1 f T (k) f (k)-γ 2 2 v T (k)v(k) < 0 (13)
where

∆V (k) = V (k + 1) -V (k).
Proof. Define the following function

J = V (∞) -V (0) + e f 2 2 -γ 2 1 f 2 2 -γ 2 2 v 2 2 (14) Note that V (∞) ≥ 0, it is obvious that the inequality J < 0 implies e f 2 2 -γ 2 1 f 2 2 -γ 2 2 v 2 2 -V (0) < 0 (15)
i.e., the criterion ( 12) is satisfied.

Using the fact that ∆V (k) = V (k + 1) -V (k), we obtain

J = ∞ k=0 ∆V (k) + e f 2 -γ 2 1 f 2 -γ 2 2 v 2 = ∞ k=0 ∆V (k) + e T f (k)e f (k) -γ 2 1 f T (k) f (k) -γ 2 2 v T (k)v(k) (16)
It can be seen that J < 0 holds if ( 13) is fulfilled for all k. Then, it can be concluded that the criterion ( 12) is satisfied if ( 13) holds for all k. ✷ Based on Theorem 1, the following Theorem is proposed to design the matrices L i , i = 1, 2, . . . , h in observer (9). Theorem 2. For given scalars γ 1 > 0 and γ 2 > 0, if there exist a symmetric positive definite matrix P ∈ R (n+q)×(n+q) , a positive scalar η > 0 and matrices W i ∈ R (n+q)×m , i = 1, 2, . . . , h such that the following LMIs hold for all i = 1, 2, . . . , h:

     N 11 0 0 0 N 15 * -γ 2 1 I q 0 0 DT P * * -γ 2 2 I m 0 N 35 * * * -ηI s ḠT P * * * * -P      < 0 ( 17
)
where

N 11 = -P + C T f C f + η Lg , Lg = L g I n 0 0 0 (18) N 15 = ĀT i P -CT W T i , N 35 = -H T W T i ( 19
)
Then the fault estimation observer ( 9) is robust against the fault variation and measurement noise, i.e. satisfying the criterion (12). Moreover, if the LMIs in (17) are solved, the matrices L i , i = 1, 2, . . . , h are determined by

L i = P -1 W i , i = 1, 2, . . . , h (20)
Proof. Define the augmented state estimation error as

e(k) = x(k) -x(k) (21)
Using ( 1) and ( 9), the error dynamic is obtained as follows

e(k + 1) = h i=1 ρ i (θ(k)) ( Āi -L i C)e(k) -L i Hv(k) + D f (k) + Ḡ Φ(k) (22)
where Φ(k) is given by Φ

(k) = Φ θ(k), C x x(k) -Φ θ(k), C x x(k) (23)
Take the following Lyapunov function V (k) = e T (k)P e(k), P > 0 (24) Substituting ( 22) into the difference of V (k) gives

∆V (k) = e T (k)A T (k)P A (k)e(k) -e T (k)P e(k) +2e T (k)A T (k)P D f (k) -L (k)Hv(k) + Ḡ Φ(k) +2 f T (k) DT P -L (k)Hv(k) + Ḡ Φ(k) + f T (k) DT P D f (k) -2v T (k)H T L T (k)P Ḡ Φ(k) +v T (k)H T L T (k)P L (k)Hv(k)+ ΦT (k) ḠT P Ḡ Φ(k) (25)
where

A (k) = h i=1 ρ i (θ(k))(T Āi -L i C) (26) L (k) = h i=1 ρ i (θ(k))L i (27) Since Φ(θ(k), x(k)) satisfies the Lipschitz condition (3), it is easy to show that η ΦT (k) Φ(k) ≤ ηe T (k) Lg e(k) (28)
where η > 0 is a positive scalar and Lg is given by Lg = L g I n 0 0 0 (29) Using ( 25) and ( 28), we obtain

∆V (k) + e T f (k)e f (k) -γ 2 1 f T (k) f (k) -γ 2 2 v T (k)v(k) ≤ e T (k) A T (k)P A (k)e(k) -P + C T f C f e(k) +2e T (k)A T (k)P G (k) f (k) -L (k)v(k) + H Φ(k) + f T (k) G T (k)P G (k) -γ 2 1 f (k) +2 f T (k)G T (k)P H Φ(k) -L (k)v(k) +v T (k) L T (k)P L (k) -γ 2 2 v T (k)v(k) v(k) -2v T (k)L T (k)P H Φ(k) + ΦT (k) HT P H Φ(k) +ηe T (k) Lg e(k) -η ΦT (k) Φ(k) = ξ T (k)Ψξ(k) (30)
where

ξ(k) =     e(k) f (k) v(k) Φ(k)     , Ψ =    Ψ 11 Ψ 12 Ψ 13 Ψ 14 * Ψ 22 Ψ 23 Ψ 24 * * Ψ 33 Ψ 34 * * * Ψ 44    (31)
where

Ψ 11 = A T (k)P A (k) -P + C T f C f + η Lg Ψ 12 = A T (k)P G (k), Ψ 22 = G T (k)P G (k) -γ 2 1 I d Ψ 13 = -A T (k)P L (k), Ψ 23 = -G T (k)P L (k) Ψ 33 = L T (k)P L (k) -γ 2 2 I q , Ψ 14 = A T (k)P H Ψ 24 = G T (k)P H, Ψ 34 = -L T (k)P H Ψ 44 = HT P H -ηI s
From (30), it can be seen that ( 13) holds if Ψ < 0. From Theorem 1, it is obtained that the criterion ( 12) is fulfilled if Ψ < 0 is satisfied.

By using Schur complement Lemma (Boyd et al. (1994)), Ψ < 0 is equivalent to the following inequality

      Ψ ′ 11 0 0 0 A T (k)P * -γ 2 1 I q 0 0 ḠT P * * -γ 2 2 I m 0 -L T (k)P * * * -ηI s HT P * * * * -P       < 0 (32)
where Ψ ′ 11 = -P + C T f C f + η Lg . Using the property of ρ i (θ(k)) in (2), the above inequality (32) holds if the following matrix inequalities hold for all i = 1, 2, . . . , h:

     M 11 0 0 0 M 15 * -γ 2 1 I q 0 0 ḠT P * * -γ 2 2 I m 0 M 35 * * * -ηI s HT P * * * * -P      < 0 ( 33
)
where

M 11 = -P + C T f C f + η Lg , M 15 = ( Āi -L i C) T P, M 35 = -L T i P
By letting W i = P L i , the inequalities in (33) become the LMIs in ( 17). Consequently, if the LMIs in ( 17) are solved, the matrices L i , i = 1, 2, . . . , h are determined by (20). ✷ Remark 1. In fact, the index γ 1 is related to the robustness against the fault variation f (k) while γ 2 represents the attenuation level of the measurement noise. However, it should be noted that the robustness against fault variation and the insensitiveness to the measurement noise are conflicted. Therefore, it is necessary to take trade-off when choose γ 1 and γ 2 in fault estimation observer design.

SIMULATIONS

In this section, a simulation example is used to show the effectiveness of the proposed method. The system under consideration is a dynamic model of a missile from Lendek et al. (2010). The continuous-time model of the considered missile is given as follows

α = K α M C n (α, δ, M )cos(α) + q (34) q = K q M 2 C m (α, δ, M ) (35)
where M denotes the Mach number, which is an exogenous scheduling variable, α is the angle of attack, q is the pitch rate and δ is the action of the tail fin actuator. It is assumed that α and q can be measured. The coefficients in (34)-( 35) are expressed as

K α = (π/180)0.7P 0 S/(mv s ), K q = (π/180)0.7P 0 Sd/I y C n (α, δ, M ) = a n α 3 + b n |α|α + c n (2 -M/3)α + d n δ C m (α, δ, M ) = a m α 3 + b m |α|α + c n (-7 + 8M/3)α + d m δ
where P 0 = 973.3 lbs/ft 2 is the static pressure at 20000 ft, S = 0.44 ft 2 is the surface area, v s = 1036.4 ft/s is the speed of sound at 20000 ft, d = 0.75 ft is the diameter, I y = 182.5 slug • ft 2 is the pitch moment of inertia and m = 13.98 slugs is the mass of the missile. Other parameters in C n (α, δ, M ) and C m (α, δ, M ) are given as follows.

a n = 0.000103 deg -3 , b n = -0.00945 deg -2 c n = 0.1696 deg -1 , d n = -0.034 deg -1 a m = 0.000215 deg -3 , b m = -0.0195 deg -2 c m = 0.051 deg -1 , d m = -0.206 deg -1 Denote x = [ α q ]
T , u = δ, the dynamic equations ( 34)-( 35) are written as

ẋ = A c (M, α)x + B c (M, α)u + Φ c (M, α) (36)
where

A c (M, α) = K α M cos(α)c n (2 -M/3) 1 K q M 2 c m (-7 + 8M/3) 0 B c (M, α) = K α M cos(α)d n K q M 2 d m Φ c (M, α) = K α M (a n α 3 + b n |α|α)cos(α) K q M 3 (a m α 3 + b m |α|α)
Notice that K α is very small, then the above system (36) can be approximated by

x(k + 1) = Ãc (M )x(k) + Bc (M )u(k) + Φ c (M, α) (37)
where

Ãc (M ) = 0 1 K q M 2 c m (-7 + 8M/3) 0 Bc (M ) = 0 K q M 2 d m
By using the Euler approximation method (Arcak and Nešić (2003); [START_REF] Mao | Fault-tolerant control for a class of nonlinear sampled-data systems via a Euler approximate observer[END_REF]), the continuous system (37) is discretized as

x(k + 1) = A(M )x(k) + B(M )u(k) + Φ(M, α) (38)
where

A(M ) = t s Ãc (M ) + I 3 B(M ) = t s Bc (M ) Φ(M, α) = t s Φ c (M, α)
and t s = 0.2 is the sampling period. In this paper, the range of M is given by

M ∈ 2, 4
Moreover, it is assumed α is bounded as α ∈ -15, 15 .

Then, it can be derived that the nonlinear function Φ(M, α) satisfies the Lipschitz condition (3) with L g = 0.04.

In order to use the proposed method, the discretized system (38) should be written as an LPV representation with a nonlinear term. Using the methodology in [START_REF] Hamdi | Fault detection and isolation in linear parameter-varying descriptor systems via proportional integral observer[END_REF], the following polytopic system is constructed to approximate system (38):

x(k + 1) = 2 i=1 ρ i (θ(k)) A i x(k) + B i u(k) + F i f (k) +GΦ(θ(k), x(k)) (39)
where the parameter vector is θ(k) = M and the weighting functions ρ i (θ(k)), i = 1, 2 are described as follows

ρ 1 (θ(k)) = M -M M -M ρ 2 (θ(k)) = M -M M -M
Herein, M and M are the extreme values of the parameter M , i.e. M = 2, M = 4

The matrices in system (39) are given by A 1 = 1.0000 0.2000 0.0033 1.0000 , A 2 = 1.0000 0.2000 -0.0292 1.0000

B 1 = 0.0000 -0.0081 , B 2 = 0.0000 -0.0322 G = 1 0 0 1 , F i = B i , i = 1, 2
In addition, it is assumed that the measurement equation is described by

y(k) = Cx(k) + Hv(k) (40)
where

C = 1 0 0 1 , H = 0.01 0 0 0.01
Now, a nonlinear system in the form of ( 1) is obtained by combining ( 39) and ( 40). Then, the proposed method can be applied to design a fault estimation observer. Remark 2. Note that the fault distribution matrices F i , which are the same as B i , are very small. From the perspective of inverse problem, the fault estimation problem is not well-posed. This implies that the fault estimation will be severely effected by the measurement noise. Therefore, it is necessary to consider the influence of measurement noise and design a robust fault estimation method.

By choosing γ 1 = 7, γ 2 = 1.4 and solving the LMIs in (17), we obtain the following matrices in the fault estimation observer ( 9) In the simulation, the initial condition of the nonlinear system is x(0) = [ 10 2 ] T . The time-varying variable M is shown in Fig. 1 and the weighting functions are depicted in Fig . 2. In addition, the measurement noise used in the simulation is zero-mean Gaussian sequences with standard deviation σ = 0.01.

L 1 = 0.
First, consider an abrupt fault which is represented by

f (k) = 0 k < 150 15 k ≥ 150 (41)
The fault estimation result of the robust fault estimation observer is depicted in Fig. 3. Therein, the actual fault is depicted by solid line and the fault estimation is represented by the dash one. As shown in Fig. 3, the proposed A new actuator fault diagnosis technique capable of estimating actuator faults for a class of nonlinear systems has been proposed in this paper. The considered nonlinear system consists of an LPV part and a Lipschtiz nonlinear term. By considering the fault as an auxiliary state vector, an augmented system is first constructed. Based on the augmented system, a robust fault estimation observer is designed to attenuate the effect of the fault variation and measurement noise. The effectiveness of the proposed method is demonstrated by a missile model.

Fig. 2 .Fig. 3 .Fig. 4 .

 234 Fig. 1. The Mach number M in the simulation

⋆ This work was financially supported by the National Natural Science Foundation of China (Grant No. 61273162) and the Fundamental Research Funds for the Central Universities (Grant No. HIT.NSRIF.2015034).