J Salvi

M Roussey

F I Baida

M.-P Bernal

A Mussot

T Sylvestre

H Maillotte

D Van Labeke

A Perentes

I Utke

C Sandu

P Hoffmann

B Dwir

Annular aperture arrays: study in the visible region of the electromagnetic spectrum

Baida and Van Labeke recently proposed a structure that exhibits a supertransmission of light through an array of nanometric coaxial apertures in a metallic film that has been named an annular aperture array (

Ebbesen et al. 1 showed in 1998 that extraordinary light transmission could be achieved through a subwavelength array of cylindrical holes in metallic films. Following this work, Baida and Van Labeke recently proposed a structure that could exhibit a higher transmission. [2][3][4] This structure consists of an array of subwavelength coaxial apertures in a real metal (gold, for instance), and it has been named an annular aperture array (AAA). The coaxial cavities offer the advantage of showing a guided mode with a cutoff wavelength in the red region of the spectrum contrary to a simple cylindrical waveguide in the same metal. 5 This kind of nanostructure could have interesting applications (filters, flat displays, nanolithography, etc.). Finite-difference time-domain (FDTD) simulations have established that a 90% transmission could be reached in the visible region of the spectrum. 3 However, the designed structure has several difficult technical constraints. First, the metal is silver, which is not well adapted to nanofabrication. Second, the period and the coaxial radii must be very small (period 300 nm, inner radius 50 nm, and outer radius 75 nm). Recently, 80% transmission through an AAA structure was experimentally established in the near-infrared part of the spectrum ͑ = 1.97 m͒. 6 We plan to build an AAA structure to have a large transmission in the visible region, and we present here our first experimental results. So far, the AAA structure has been fabricated in gold by electron lithography and gold lift-off. We present images of the fabricated structure by electronic microscopy and atomic force microscopy (AFM). We experimentally record the structure's far-field spectral response by using a bright fibered light source that presents a continuum spectrum in the visible range.

The numerical simulation uses a homemade FDTD code already described elsewhere. [2][3][4] The 3D calculations use Berenger's perfectly matched layers (PMLs) boundary conditions to cancel the parasitic reflections on the edges of the calculation window. The latter technique treats only the homogeneous waves. 7 Therefore, to avoid unwanted reflections caused by evanescent waves, the object is placed at a distance from the edges larger than half of the highest value of the illuminating wavelength. In addition, to obtain accurate results, the spatial meshing of the object must be small, which implies a great number of spatial nodes. To decrease the number of nodes, we use a nonuniform spatial discretization of the metallic structure as described in Ref. 8.

In the present study, the object consists of a biperiodic (in the x and y directions) array of coaxial gold apertures. The metallic layer is deposited on a glass substrate that is also incorporated into our simulations. The calculation window consists of one array period (see Fig. 1, where four periods are presented).

The absorbing PMLs are therefore necessary only in the z direction. Periodic conditions are applied in the x and y directions. The spatial increment step is set to ⌬ x = ⌬ y = 10 nm. The nonuniform meshing is done in the z direction. A spatial step of ␦ z =10 nm is also chosen to mesh the metallic film, whereas ⌬ z = 25 nm for the rest of the window. This mesh variation must be carried out gradually to avoid artificial reflections. [START_REF] Taflove | Computational Electrodynamics: the Finite-Difference Time-Domain Method[END_REF] For the sample under study, the periodicity of the coaxial apertures is p = 600 nm and the gold-layer thickness is h = 150 nm. The incident field is a pulsed plane wave polarized in the x direction and propagating along the z axis (normal incidence from the glass substrate). It is centered on c = 650 nm, and it covers the spectral range between 530 and 900 nm. For real metal, like gold, silver, or aluminum, the dispersion of the dielectric permittivity must be taken into account. By use of Drude's analytical model, a simple computing scheme connecting the displacement vector to the electric field is obtained. It can be easily incorporated into the usual leap-frog FDTD scheme. [START_REF] Baida | [END_REF] For this work, we studied a broad spectral range, and therefore the dispersion of gold cannot be described by only one simple Drude's model. To bypass this problem, the spectral range is subdivided into small spectral intervals of ␦ = 10 nm. In every spectral domain, the tabulated experimental permittivities 11 are fit with an appropriate Drude's model and a FDTD calculation is performed.

The transmission is obtained by evaluating the zero-order intensity normalized by the incident intensity. In fact, because of the periodicity of the structure ͑p = 600 nm͒, only one diffracted order (zero order) exists for wavelengths larger than p. For wavelengths between = 532 nm (beginning of our continuum) and = 600 nm, higher diffracted orders appear ͑±1͒, but they are not experimentally recorded because they propagate far from the fiber axis (see below).

The 2D arrays of coaxial apertures are generated by electron-beam lithography and gold lift-off. The gold structures are made on a 20 nm thick solgelcoated transparent conductive oxide layer on a microscope slide. The resist is spin coated to 700 nm thickness to obtain the necessary high-aspect-ratio resist structures. The electron lithography equipment is based on a JEOL 6400 JSM scanning electron microscope (SEM) with a LaB 6 thermionic emitter. The electron beam is controlled with Nabity Systems patterning software (NPGS). The exposure dose is varied between 200 and 250 C/cm 2 . With 250 C/cm 2 the structures obtained correspond best to the calculated geometry (see Fig. 2).

After developing in the surfactant containing the developer, oxygen plasma treatment, and adhesionpromoting titanium evaporation coating (5 nm thick), 150 nm of gold is deposed by electron-beam evaporation. Finally, gold lift-off is achieved in an ultrasonic acetone bath.

SEM and AFM images of one of the best-looking arrays are presented in Fig. 2. The images consist of a large array of coaxial apertures with inner and outer diameters of 250 and 330 nm, respectively, and a grating period of p = 600 nm. Note that the apertures are well aligned and that the AAA does not present any important defects. The AFM line scan does not reveal the full depth of the cylindrical holes because of the cone angle of the AFM tip.

The experimental setup is shown in Fig. 1. A local illumination, under normal incidence, is obtained by use of a powerful white-light fiber source. 12 The laser output beam is focused into a 100 m dispersionshifted fiber (DSF) generating a broadband singlemode supercontinuum by stimulated Raman scattering and parametric four-wave mixing. To collect the light transmitted through the AAA, a cleaved-end multimode fiber (62.5 m core diameter) is brought to around 20 m close to the sample. The spectrum is then recorded by means of an optical spectrum analyser (S2000 Miniature Fiber Optic Spectrometer from Ocean Optics).

Figure 3 shows an image of a zone of the sample. It consists of a large area containing four AAAs and one triangular metal-free test zone. The image is obtained by means of a classical optical microscope used in transmission with white light. The color modification reflects the variations of the spectral properties of AAA structures when changing its geometric parameters. To obtain the spectral response of one matrix of AAA (see Fig. 2), we have to record the source spectrum through the metal-free test zone.

Figure 4 shows the experimental spectral response of an AAA zone presented in Fig. 2 compared to the theoretical spectral response obtained by the FDTD calculation. We would like to point out that, apart from a background subtraction, no image treatment was performed on the curve. The noise appearing in the experimental spectrum is probably induced by temporal instability of the laser and mechanical instabilities in the injection and/or detection processes.

One can see from Fig. 4 that there is good agreement between experience and theory. For both curves a transmission peak around 700 nm (theoretically around 16%) is observed. The experimental curve is slightly translated toward higher wavelengths. This shift is probably linked to the polarization direction of the incident beam (see Fig. 7 in Ref. 2). In the theoretical simulation, the polarization was set along the x or y direction. The experimental setup could not allow us to control this parameter. Moreover, the fabrication parameters slightly differ from one coax with respect to the other and from the theoretical ones. Finally, the calculation does not take into account the very thin titanium layer (5 nm thick) that is used to improve the adhesion of gold onto the substrate.

In this Letter we have presented the first study, to our knowledge, of transmission through an AAA structure in the visible region of the spectrum. A FDTD calculation reproduces quite well the experimental spectrum. Nevertheless, in this first attempt, the physical parameters of the structure only lead to a maximum of transmission around 17% for = 700 nm. In fact, more than 60% of transmission for the fabricated structure is theoretically found at = 1330 nm (see the inset in Fig. 4). The peak at = 700 nm is then the first harmonic of the TE 11 mode as demonstrated in Ref. 5.

The fabrication of other structures is in progress at the École Polytechnique Fédérale de Lausanne. The challenge is difficult as it needs a strong enhancement of our technological process to reduce all the dimensions (period and radii) by a factor of 2. Moreover, the theoretical study of the transmission of light by AAA structures remains an interesting subject. For instance, the transmission spectrum obtained in this Letter is different from those published in Refs. 2 and 3. This difference has to be studied in more detail and explained in order to design new structures and to optimize their parameters.

Fig. 1 .

 1 Fig. 1. Scheme of the studied structure and experimental setup.

Fig. 2 .

 2 Fig. 2. (a) SEM and (b) AFM ͑10 m ϫ 10 m͒ images of the AAA under study. (c) Cross section (in red) made over four annular apertures shown by the red line in (b).

Fig. 4 .

 4 Fig. 4. Theoretical (red) and experimental (blue) spectral responses of the studied AAA structure located by the white arrow in Fig. 3. Inset, theoretical transmission in the nearinfrared region showing a maximum of transmission at = 1330 nm.

Fig. 3 .

 3 Fig. 3. Far-field image of one sample obtained with a conventional microscope. The white arrow indicates the AAA for which the spectral response has been studied.

0146-9592/05/131611-3/$15.00 © 2005 Optical Society of America

This research was supported by the Laboratoire Européen Associé en Microtechnique Franco-Suisse, and we also thank the Swiss National Science Foundation (N 200021-103333) for partial financial support of this work. F. I. Baida's e-mail address is fbaida@univ-fcomte.fr.


Mathieu Pouliquen
email: mathieu.pouliquen@unicaen.fr

Eric Pigeon

Olivier Gehan

Recursive Identification for Hammerstein Systems with Bounded Noise

In this paper, we consider the identification of Hammerstein systems in presence of bounded disturbances. The proposed identification algorithm is a recursive algorithm. It ensures that the estimation error is bounded and allows an online separation of the nonlinear part and linear part. Furthermore, some stability and convergence results are demonstrated and a simulation example is given to illustrate the performances of the presented method.

I. INTRODUCTION

The Hammerstein systems are the nonlinear systems consisting of a static nonlinearity block followed by a linear dynamic system. The identification of such systems has received much interest for the last decades ( [START_REF] Narendra | An iterative method for the identification of nonlinear systems using a hammerstein model[END_REF], [START_REF] Chang | A noniterative method for identification using hammerstein models[END_REF], [START_REF] Haist | Nonlinear identification in the presence of correlated noise using a hammerstein model[END_REF], [START_REF] Greblicki | Identification of discrete hammerstein systems using kernel regression estimates[END_REF], [START_REF] Bai | An optimal two-stage identification algorithm for hammersteinwiener nonlinear systems[END_REF], [START_REF] Bai | Identification of linear systems with hard input nonlinearities of known structure[END_REF], [START_REF] Cerone | Parameter bounds for discrete-time hammerstein models with bounded errors[END_REF], [START_REF] Ding | Identification of hammerstein nonlinear armax systems[END_REF], [START_REF] Zhao | Recursive identification for hammerstein system with arx system[END_REF], [START_REF] Giri | Identification of hammerstein systems in presence of hysteresis-backlash and hysteresisrelay nonlinearities[END_REF], [START_REF] Zhao | Adaptive tracking and recursive identification for hammerstein systems[END_REF], [START_REF] Bai | Convergence of the iterative algorithm for a general hammerstein system identification[END_REF], [START_REF] Giri | Block Oriented Nonlinear System Identification[END_REF], [START_REF] Le | Recursive identification of hammerstein systems with application to electrically stimulated muscle[END_REF], [START_REF] Han | Hammerstein system identification using nuclear norm minimization[END_REF] and the four sessions devoted to block oriented nonlinear identification at SYSID 2012). The reason is that these systems take into account nonlinearities commonly encountered in practice, generated by the technological limitations necessary for the proper functioning of the system (saturation, the limit stops, etc.).

For Hammerstein systems or not, all practical identification algorithms have to deal with measurements corrupted by noise and most of the previous works consider the stochastic noise assumption (except [START_REF] Cerone | Parameter bounds for discrete-time hammerstein models with bounded errors[END_REF] discussed later). In many cases (unknown probability distribution of the disturbances, modeling inaccuracy) such assumption cannot be done and the bounded noise assumption seems most appropriate. To address this identification problem we will use in this paper Set Membership Identification (SMI) algorithms (see [START_REF] Milanese | Estimaton theory and uncertainty intervals evaluation in presence of unknown but bounded errors. linear families of models and estimators[END_REF], [START_REF] Lozano-Leal | Reformulation of the parameter identification problem for systems with bounded disturbances[END_REF], [START_REF] Vicino | Sequential approximation of parameter sets for identification with parametric and nonparametric uncertainty[END_REF], [START_REF] Garulli | Error bounds for conditional algorithms in restricted complexity set membership identification[END_REF] and [START_REF] Vicino | Conditional central algorithms for worst-case set membership indentification and filtering[END_REF]) and particularly the Optimal Bounding Ellipsoid (OBE) algorithms. These algorithms have a low computational complexity and they are appropriated to handle the identification problem in presence of bounded disturbances. Their principle consists in the estimation of a feasible set of parameters which must be consistent with the measurement data and the model structure. Some important contributions have been presented in [START_REF] Fogel | On the value of informaton in system identification -bounded noise case[END_REF], [START_REF] Dasgupta | Asymptotically convergent modified recursive least square with data-dependent updating and forgetting factor for systems with bounded noise[END_REF], [START_REF] Canudas-De-Wit | A modified EW-RLS algorithm for systems with bounded noise[END_REF] and [START_REF] Tan | Identification for systems with bounded noise[END_REF].

The aim of this paper is to propose such recursive identification algorithm for Hammerstein systems subjected to bounded noises. The proposed algorithm is a two stages algorithm: the first stage updates a parameters vector from the new measurements, the second stage separates recursively the linear and nonlinear part. This separation induces an internal disturbance term and a great difficulty is to demonstrate the stability of the algorithm even with this internal disturbance term. This is one of the main contributions of this paper.

Briefly, the paper is organized as follows: section II introduces the model structure and states the objective of the study. In section III, the adaptation algorithm is proposed and some stability and convergence properties are addressed. Some simulation results are given in section IV. Finally, section V concludes the paper. Appendices contain most of the proofs.

II. PROBLEM STATEMENT

The Hammerstein system presented in figure 1

is described by A(q)y t = q -d B(q)x t + v t (1)
with

B(q) = b 0 + b 1 q -1 + • • • + b n b q -n b A(q) = 1 + a 1 q -1 + • • • + a n a q -n a and x t = h(u t ) (2)
where u t and y t are respectively the system input and the system output. x t is the internal signal and v t is an unknown but bounded disturbing term (noise measurement, state disturbances, modeling inaccuracy, etc.). h(.) is a static nonlinear function without memory. It is represented by a sum of the known nonlinear functions h 1 (.), h 2 (.), ..., h n h (.):

h(u t ) = µ 1 h 1 (u t ) + µ 2 h 2 (u t ) + • • • + µ n h h n h (u t ) (3)
✲ h(.) vt yt xt q -d B(q) ✲ 1/A(q) ✲ ❄ ✲ ✲ ut + + Fig. 1. Studied Hammerstein system
The objective: Given the delay d, the degrees n a , n b and n h , estimate recursively the parameters {a i }, {b i } and {µ i } of the model based on the input-output measurements of the system described by (1), ( 2) and (3).

Remark 1: In [START_REF] Cerone | Parameter bounds for discrete-time hammerstein models with bounded errors[END_REF] a similar problem is treated. The difference lies in the fact that in [START_REF] Cerone | Parameter bounds for discrete-time hammerstein models with bounded errors[END_REF] the linear part is a OE model, the algorithm is not a recursive algorithm, several experiments are needed to identify some points of h(.) and the separation of the linear and nonlinear part is not realized in the same time.

In this paper, so as to be able to identify the parameters, some assumptions are needed: A1 The delay d, the orders of the linear part and of the nonlinear part (n a , n b and n h ) are supposed to be available; A2 The coefficient b i are such that b 0 = 1; A3 Regardless of the disturbance nature, an upper bound is supposed to be available on the magnitude of v t :

|v t | ≤ δ v t ( 4
)
Remark 2:

The assumption [A2] is necessary so as to satisfy the identifiability condition. Without this assumption we are not able to distinguish the solution G(q)h(u t ) from the solution (aG(q))( 1a h(u t )) with a = 0.

Remark 3: About [A3], it is not necessary to know exactly the value of δ v t . When identifying a system for which this bound is unknown, it is possible to process by dichotomy to refine the results of the estimation.

III. IDENTIFICATION ALGORITHM AND ANALYSIS

From (1), ( 2) and (3) the output of the Hammerstein system is expressed as follows:

y t = n b ∑ k=0 b k n h ∑ i=1 µ i h i (u t-k-d ) - n a      h 1 (u t-d ) . . . h 1 (u t-d-n b ) h 2 (u t-d ) . . . h 2 (u t-d-n b ) . . . h n h (u t-d ) . . . h n h (u t-d-n b ) -y t-1 . . . -y t-na                                       (5)
the system is rewritten in the compact form

y t = φ T t θ * + v t (6)
θ * ∈ R n is the unknown parameters vector to be identified with n = (1 + n b )n h + n a the number of parameters. The first (1 + n b )n h terms of the vector θ * can be inserted in a (1 + n b ) × n h matrix as follows

M θ * =    θ * (1) • • • θ * ((1 + n b )(n h -1) + 1) . . . . . . θ * (1 + n b ) • • • θ * ((1 + n b )n h )   
which can be rewritten as

M θ * =      1 b 1 . . . b n b      µ 1 • • • µ n h (7)
Such structure will be used to separate the parameters {b i } and {µ i }. By now, our aim in what follows is to

• propose an estimation algorithm based on ( 6) and (4);

• estimate a parameters vector structured as [START_REF] Cerone | Parameter bounds for discrete-time hammerstein models with bounded errors[END_REF];

• give a proof of stability and convergence of this algorithm.

The algorithm considered in this paper recursively estimates the parameters θ * . For this purpose the one-step-ahead output predictor is defined as follows ( [START_REF] Ljung | System identification: theory for the user[END_REF]):

ŷt = y t -v t = φ T t θ *
Before introducing the parameter adaptation algorithm, let define the a priori and a posteriori predictors as ŷt, θt-1 = φ T t θt-1 ŷt, θt = φ T t θt θt represents an estimation of the parameters vector at the actual time t. Let define also θ s t and θ us t as follows: • θ s t is an estimation which satisfies (7); • θ us t is an estimation which does not satisfy [START_REF] Cerone | Parameter bounds for discrete-time hammerstein models with bounded errors[END_REF]. Subscripts s and us respectively stand for "structured" and "un-structured".

A. The adaptation algorithm.

To address the identification problem we will take inspiration from the OBE identification methods presented in [START_REF] Canudas-De-Wit | A modified EW-RLS algorithm for systems with bounded noise[END_REF], [START_REF] Tan | Identification for systems with bounded noise[END_REF] and [START_REF] Sun | Comments on "identification for systems with bounded noise[END_REF]. As depicted on Fig. 2, at the time t the proposed algorithm is divided in two successive stages: stage 1 Update of the estimated parameters vector taking into account the new observation vector φ t and the new output y t : estimation of θ us t from θ s t-1 ;

stage 2 Structuring of the estimated solution: estimation of θ s t from θ us t .

The estimated parameters vector must satisfy [START_REF] Cerone | Parameter bounds for discrete-time hammerstein models with bounded errors[END_REF] and has to maintain the prediction error below a bound defined from the upper bound on the disturbance v t . Under some conditions describes below, the following estimation algorithm provides such estimation:

stage 1            θ us t = θ s t-1 + Γ t ỹt, θ s t-1 Γ t = P s t-1 φ t σ t λ +φ T t P s t-1 φ t σ t ỹt, θ s t-1 = y t -ŷt, θ s t-1 (P us t ) -1 = λ (P s t-1 ) -1 + φ t σ t φ T t (8) stage 2 θ s t = S( θ us t ) P s t = 1 1-ζ P us t + 1 ζ W t (9)
Some useful insights on the algorithm are given below. about stage 1: θ s t-1 and P s t-1 are computed at the previous iteration. ỹt, θ s t-1 is the a priori prediction error, 0 < λ ≤ 1 is the forgetting factor and σ t is a switching flag given by:

σ t =            λ φ T t P s t-1 φ t ỹt, θ s t-1 δ t -1 if ỹt, θ s t-1 δ t > 1 and φ T t P s t-1 φ t > 0 0 otherwise (10) Its role is specified below. From (8) the unstructured a posteriori prediction error ỹt, θ us t = y t -ŷt, θ us t can be written as ỹt, θ us t = (1 -φ T t Γ t ) ỹt, θ s t-1
which can also be rewritten, by using the expression for Γ t , as:

ỹt, θ us t = λ λ + φ T t P s t-1 φ t σ t ỹt, θ s t-1 (11)
Using the value of σ t for σ t = 0 (11) yields:

ỹt, θ us t δ t = 1
This clearly shows that, provided that φ T t P t-1 φ t > 0, the introduction of σ t in the first part of the algorithm ensures the following property:

       if ỹt, θ s t-1 δ t > 1 then σ t = 0 and ỹt, θ us t δ t = 1 if ỹt, θ s t-1 δ t ≤ 1 then σ t = 0 and ỹt, θ us t δ t ≤ 1
δ t is then a bound on the magnitude of the unstruc- tured a posteriori prediction error. This bound has to be specified taking into account the bound δ v t .

about stage 2:

The second part of the algorithm as to realize the structuring of θ us t . This is done by a function S(.): θ s t = S( θ us t ). This function will be detailed in subsection III-C (see table I page 5).

The interest of matrix W t is to maintain the stability of the algorithm. In fact, the use of that structuring step introduce a correction term w t in the following manner w t = θ s t -θ us t and that correction term disturbs the algorithm behavior. So as to take into account this disturbance W t is designed in the following manner:

W t : W t = 0 if w t = 0 W t such that w T t W -1 t w t ≤ 1 otherwise (12) 0 < ζ < 1 is
a design parameter whose role will be discussed later.

Let notice that without the second stage, the above algorithm is similar to OBE algorithms presented in [START_REF] Canudas-De-Wit | A modified EW-RLS algorithm for systems with bounded noise[END_REF], [START_REF] Tan | Identification for systems with bounded noise[END_REF] and [START_REF] Boutayeb | Recursive identification of linear multivariable systems with bounded disturbances[END_REF]. These algorithms can be apply to any processes satisfying

y t = φ T t θ + v t with |v t | ≤ δ v t .
The main interest of the proposed algorithm is the introduction of the correction term w t in the second stage. In this paper this correction term corresponds to a structuring procedure of the estimated parameters vector. It is possible to imagine an other identification problem in which the correction term will have a different role. In the next subsection it is shown that stability and convergence properties can be demonstrated even in presence of that correction term.

B. Stability and convergence properties.

The analysis of the estimation algorithm is performed via the analysis of the estimation error θ s t defined by: θ s t = θ * -θ s t

In the first part of the analysis some geometrical interpretations of the adaptation algorithm are provided. To this end let us define for each time t the observation set S t :

S t = θ ∈ R n , y t -φ T t θ δ t ≤ 1 (13)
Given y t and φ t , S t is the set of all possible θ which are consistent with the chosen threshold δ t . An important property of that observation set is given in the following propertie.

Propertie 1: Assume that δ t ≥ δ v t , then

θ * ∈ S t ( 14
)
S t is then a strip in parameter space bounded by two parallel hyperplanes and which contains θ * . The principe of OBE type algorithms is to find a parameters vector included in t i S i . Let define the un-structured ellipsoid E us t :

E us t = θ ∈ R n , (θ -θ us t ) T P us -1 t (θ -θ us t ) ≤ 1 and the structured ellipsoid E s t : E s t = θ ∈ R n , (θ -θ s t ) T P s -1 t (θ -θ s t ) ≤ 1 (15)
The following theorem shows that the proposed identification algorithm provides rules for computing θ s t and P s t so that E s t contains the intersection (S t E s t-1 ) (Fig. 3).

S t

E us t E s t-1 E s t Fig. 3. (S t ∩ E s t-1 ) ⊂ E us t ⊂ E s t
Theorem 1: Consider the class of systems defined in section II and the identification algorithm given by ( 8), ( 9), [START_REF] Ding | Identification of hammerstein nonlinear armax systems[END_REF] and [START_REF] Garulli | Error bounds for conditional algorithms in restricted complexity set membership identification[END_REF]. If δ t ≥ δ v t hold and suppose that θ * ∈ E s t-1 [START_REF] Haist | Nonlinear identification in the presence of correlated noise using a hammerstein model[END_REF] then for all initial conditions

• a-An outer bounding ellipsoid of (S t ∩ E s t-1 ) is given by

E s t ; • b- θ * ∈ E s t ( 17
)
If furthermore, {φ t } is a persistently exciting sequence of order o e ≥ n, i.e there exist α > 0 and β > 0 such that for all t ≥ o e

αI n ≤ o e -1 ∑ i=0 φ t-i σ t-i φ T t-i ≤ β I n ( 18
)
and if there exists β W > 0 such that for all t

0 ≤ W t ≤ β W I n ( 19
)
Then the following property holds:

• c-The volume of the ellipsoid E s t is bounded.

A sketch of proof is given in appendix A page 7.

A consequence of this result is the fact that if in the beginning (i.e. at the time t = 0) θ * is included in E s 0 (this is true with for example θ s 0 = 0 and P s 0 = µI n with µ large enough) then θ * belongs to E s t for all t.

The user parameter ζ determines the shape of the ellipsoid E s t . Different policies can be considered so as to obtain a particular effect on this ellipsoid (to minimize the product of the length of its axes or to minimize the sum of the squares of the lengths of the axes). It can also be chosen so as to adapt the dynamic of θ s t , this will be shown in the next theorem (see [START_REF] Lozano-Leal | Reformulation of the parameter identification problem for systems with bounded disturbances[END_REF] and ( 22)).

The second part of the analysis is performed via the analysis of the estimation error θ s t which is a non linear function of θ s t-1 , w t and v t . Theorem 2 focuses on the behavior of θ s t 2 . It is shown that under some conditions the center θ s t of the ellipsoid E s t converges next to the solution θ * .

Theorem 2: Consider the class of systems defined in section II and the identification algorithm given by ( 8), ( 9), [START_REF] Ding | Identification of hammerstein nonlinear armax systems[END_REF] and [START_REF] Garulli | Error bounds for conditional algorithms in restricted complexity set membership identification[END_REF]. If δ t ≥ δ v t hold and if conditions ( 18) and (19) in Theorem 1 hold, then for all initial conditions

• a-for all t ≥ o e θ s t 2 ≤ γ 1 (1 -ζ ) t λ t θ s 1 2 + γ 2 1 -(1 -ζ ) t λ t 1 -(1 -ζ )λ (20)
with γ 1 and γ 2 two finite scalars.

If, furthermore, there exists a time t s such that for all t ≥ t s w t = 0 [START_REF] Milanese | Estimaton theory and uncertainty intervals evaluation in presence of unknown but bounded errors. linear families of models and estimators[END_REF] the following properties hold

• b-for all t ≥ t s + o e θ s t 2 ≤ γ 3 (1 -ζ ) t-t s λ t-t s ( 22
)
with γ 3 a finite scalar.

• c-one has

lim t→∞ ỹt, θ s t-1 δ t ≤ 1 (23)
A sketch of proof is given in appendix B page 8.

Remark 4: Note that if δ t = δ v t and w t = 0 then from (22) we have lim

t→∞ θ s t = θ
where θ parametrizes a model such that:

A(q)y t = q -d B(q) x t + ỹt x t = h(u t ) | ỹt | ≤ δ v t ( 24
)
which is coherent with the objective stated in section II.

C. Design of function S(.).

The aim of this subsection is to propose a procedure for the design of a structured estimated parameters vector θ s t from θ us t in the case where θ us t is not structured as [START_REF] Cerone | Parameter bounds for discrete-time hammerstein models with bounded errors[END_REF]. Such procedure corresponds to function S(.): The proposed function S(.) corresponds to the 5 steps procedure given in table I.

θ s t = S(
The interest of this procedure is twofold:

• It provides structuring of θ s t . • θ s t is such that the magnitude of the a posteriori prediction error ỹt, θ s t is bounded by δ t .

The condition on the existence of a finite scalar β W is difficult to guarantee in a such recursive algorithm. The procedure of table I is ill-conditioned if the parameters {µ i } extracted in step 1 are close to zero. This is the reason why we advice the use of the above procedure and the use of the second stage of the algorithm after a few times. Typically the implementation of the procedure must be realized only after n a + n b + n h samples (this corresponds to the number of parameters to estimate). Numerous simulations have shown the efficiency of that procedure. This is illustrated in the next section.

Function S(.) -Estimation of a structured parameters vector step 1:Extract the parameters {a i } and {µ i } from θ us t .

step 2:If ∑ n h i=1 µ 2 i = 0, compute a first estimation of P =    b 1 . . . b n b    from θ us t .
This estimation is denoted P (1) . If ∑ n h i=1 µ 2 i = 0, take P (1) = 0.

step 3:Compute the predictor ŷt = ϕ T t P (1) + z t .

step 4:Realize a second estimation P (2) as follows:

P (2) = P (1) + ϕ t σ ′ t 1 + ϕ T t ϕ t σ ′ t (y t -ŷt ) with σ ′ t =        1 ϕ T t ϕ t y t -ŷt δ t -1 if y t -ŷt δ t > 1 and ϕ T t ϕ t > 0 0 otherwise ( 25
)
step 5:Finally θ s t is the parameterized vector built from the estimates for {a i } (step 1), {b i } (step 4) and {µ i } (step 1).

TABLE I

Remark 5: It is possible to make use of a Singular Values Decomposition (SVD) in step 2. The algorithm will have a higher computational complexity, however it will be numerically more robust when dealing with ill-conditioned data.

IV. NUMERICAL EXAMPLE Numerical data have been generated according to (1), ( 2), ( 3) and ( 4). The following system has been considered in [START_REF] Bai | Convergence of the iterative algorithm for a general hammerstein system identification[END_REF]:

q -d B(q) = q -1 (1 + 0.5q -1 + 0.5q -2 ) A(q) = 1 -0.3q -1 -0.2q -2 + 0.3q -3 h(u t ) = 0.9545h 1 (u t ) + 0.2983h 2 (u t )
with h 1 (u t ) = u t and h 2 (u t ) = 2u 2 t . The input u t is white noise uniformly distributed in [-1; 1] of length N = 4000. v t is a noise generated in the following manner:

v t = δ v t 1 2 (e t + r t )
where r t is a binary square signal with a period of 20 samples, e t is a white noise uniformly distributed in [-1; 1] and δ v t is close to 0.7 in order to have a signal to noise ratio close to SNR= 5dB. The parameters {a i }, {b i } and {µ i } are estimated with three procedures:

• Procedure 1: this first procedure corresponds to the identification of an extended ARX structure using (6) and the ordinary least square algorithm. In order to separate the parameters {b i } and {µ i } a SVD is applied at each iterations.

• Procedure 2: this second procedure corresponds to our identification algorithm without stage 2. However the parameters {b i } and {µ i } are separated using the procedure presented in subsection III-C. In fact this amounts to not use the structured parameters vector estimated at time t for time t + 1. The threshold δ t is taken equal to δ v t . • Procedure 3: this third identification procedure is our identification algorithm using stages 1 and 2. Here again δ t = δ v t . In Fig. 4, 5 and 6 the thick solid lines denote the true parameters and the other lines the parameters estimates. These figures present the estimations for a 1 , b 1 and µ 1 . The other components have similar behavior. In Fig. 7, 8 and 9 the solid lines denote thresholds ±δ v t , dotted lines correspond to the a posteriori prediction error for each identification procedure.

From these figures we conclude the following:

• Procedure 1 is the only procedure for which the a posteriori prediction error is not bound by δ t . • Procedure 1 clearly delivers biased estimates (see Fig. 4 and5). This is due to the fact that the disturbing term v t is not a white noise, the ordinary least square algorithm is not adapted to the considered identification problem. • The estimations computed with procedures 2 and 3 converge towards a neighborhood of the true parameters. The difference between procedures 2 and 3 lies in the introduction of the correction term w t in procedure 3. This correction term affects the behavior of the algorithm without loss of the stability property as demonstrated in subsection III-B. • An interesting point is the fact that the convergence towards the true parameters is faster with procedure 3. This is a consequence of taking into account the specific structure of the parameters vector i.e. the use of stage 2 (here again the effect w t and W t ).

V. CONCLUSION AND FUTURE WORKS

A recursive algorithm for Hammerstein systems identification has been proposed. This algorithm is adapted for the bounded noise case. It consists in two successive stages: the first stage realizes an update of the estimate from new measurements, the second stage allows an online separation of the nonlinear part and the linear part. Stability and convergence properties have been investigated. The second stage of the algorithm consists in the introduction of a "disturbance term" on the parameters computed in the first stage. This "disturbance term" is used here in order to separate the nonlinear part and the linear part, however it could be used so as to impose other constraints on the parameters (impose some parameters to be positive for example). Future researches will concentrate on the use of that two stages algorithm on some other applications. In this paper the linear part is an ARX process, we want also to generalize our results on systems such that the linear part is parametrized by an OE process.

A-PROOF OF THEOREM 1 For lack of space we do not detail all the proof, we will simply state the essential points.

• Sketch of proof for S t ∩ E s t-1 ⊂ E us t . Consider θ such that θ ∈ E s t-1 and θ ∈ S t . From ( 13) and [START_REF] Greblicki | Identification of discrete hammerstein systems using kernel regression estimates[END_REF] this gives

λ (θ -θ s t-1 ) T P s -1 t-1 (θ -θ s t-1 ) + (y t -φ T t θ ) 2 σ t ≤ λ + ỹ 2 t, θ s t-1 σ t δ t ỹt, θ s t-1 2
This can be rewritten as: • Sketch of proof for S t ∩ E s t-1 ⊂ E us t . This result is based on the following lemma from [START_REF] Schweppe | Uncertain dynamic systems[END_REF]: Lemma 1: [START_REF] Schweppe | Uncertain dynamic systems[END_REF] Consider two ellipsoids

(θ -θ s t ) T P s -1 t (θ -θ s t ) ≤ λ + ỹ 2 t, θ s t-1 σ t δ t ỹt, θ s t-1 2 - λ σ t λ +φ T t P s t-1 φ t σ t If σ t = 0 then (θ -θ s t ) T P s -1 t (θ -θ s t ) = λ
E i = θ ; (θ -θ ci ) T P -1 i (θ -θ ci ) ≤ 1 and their sum E 1 ⊕ E 2 = {θ = θ 1 + θ 2 ; θ 1 ∈ E 1 and θ 2 ∈ E 2 }. Then ∀ζ such that 0 < ζ < 1, we have E 3 ⊇ E 1 ⊕ E 2 with E 3 = θ ; (θ -θ c3 ) T P -1 3 (θ -θ c3 ) ≤ 1 , θ c3 = θ c1 + θ c2 and P 3 = P 1 1-ζ + P 2 ζ .
Here we take P 1 = P us t , P 2 = W t and P 3 = P s t and θ c3 = θ c1 + θ c2 corresponds to θ s t = θ us t + w t .

• From the two previous results it is easy to conclude

S t ∩ E s t-1 ⊂ E us t ⊂ E s t , this corresponds to proposition a-.
Using assumption [START_REF] Haist | Nonlinear identification in the presence of correlated noise using a hammerstein model[END_REF] and the fact that δ t ≥ δ v t , propo- sition b-is a consequence of proposition a-.

P s -1 t ≤ (1 -ζ )λ P s -1 t-1 + (1 -ζ )φ t σ t φ T t P s -1 t ≤ (1-ζ ) t λ t P s -1 1 +(1-ζ ) t-1 ∑ i=0 (1-ζ ) i λ i φ t-i σ t-i φ T t-i
(26) We have 0 < ζ < 1 and 0 < λ ≤ 1, consequently there exist s sup such that P s -1 t ≤ s sup I n .

• Sketch of proof that there exists a scalar s in f = 0 such that P s -1 t ≥ s in f I n . From (9) we have

P s t = 1 1 -ζ P us t I n + 1 -ζ ζ P us -1 t W t Let us denote J t = I n + 1-ζ ζ P us -1 t W t , this gives together with (8) P s -1 t = (1 -ζ )λ J -1 t λ P s -1 t-1 + (1 -ζ )J -1 t φ t σ t φ T t
After o e iterations we get 26) and ( 19) P s -1 t W t ≤ s sup β W I n , then there exists a scalar s J such that I n ≥ J -1 t ≥ s -1 J I n . Together with (27) we show

P s -1 t = (1 -ζ ) o e λ o e ∏ o e -1 i=0 J -1 t-i λ P s -1 t-o e + ∑ o e -1 i=0 (1 -ζ ) i λ i ∏ i k=0 J -1 t-k (1 -ζ )λ φ t-i σ t-i φ T t-i (27) From (
P s -1 t ≥ (1 -ζ ) o e -1 ∑ i=0 (1 -ζ )λ i s -i J φ t-i σ t-i φ T t-i
It can be shown that there exists a scalar η such that

P s -1 t ≥ η o e -1 ∑ i=0 φ t-i σ t-i φ T t-i
together with condition (18) we get P s -1 t ≥ ηαI n = 0, then there exists a scalar s in f = 0 such that P s -1 t ≥ s in f I n .

• Proposition c-is a consequence of the fact that P s t is bounded.

B-PROOF OF THEOREM 2

For lack of space we do not detail all the proof, we will simply state the essential points. Consider the following Lyapounov function: (29)

V s t =
• Sketch of proof proposition a-From ( 28) and (29) we obtain

V s t ≤ (1 -ζ ) t λ t V s 1 + ζ t-1 ∑ i=0 (1 -ζ ) i λ i
We have shown that 0 < s in f I n ≤ P s -1 t ≤ s sup I n , then

s in f θ s t 2 ≤ (1 -ζ ) t λ t V s 1 + ζ t-1 ∑ i=0 (1 -ζ ) i λ i
Proposition a-is a consequence of this inequality with γ 1 = µ max (P s -1 1 ) s in f and γ 2 = ζ s in f .

• Sketch of proof proposition b-If there exists a time t s such that for all t ≥ t s w t = 0, it can be shown that

V s t ≤ (1 -ζ ) t-

Fig. 2 .

 2 Fig. 2. A two stages algorithm

Fig. 4 .

 4 Fig. 4. Simulation results component a 1

Fig. 5 .Fig. 6 .Fig. 7 .Fig. 8 .Fig. 9 .

 56789 Fig. 5. Simulation results component b 1

  s in f . With 0 < λ ≤ 1 and 0 < ζ < 1, proposition c-follows from proposition b-.

  W such that W t ≤ β W I n is needed in Theorem 1 and Theorem 2.

	W t :	W t = 0 if w t = 0 W t such that w T t W -1 t w t ≤ 1 otherwise
	The existence of a scalar β • θ s t must satisfy ỹt, θ s t δ t ≤ 1.

θ us t ) This function S(.) must satisfy two conditions:

• The correction term w t = θ s t -θ us t must be such that

•

  Sketch of proof that there exists a scalar s sup such that P s -1

	≤ s sup I n . From (9) we have P s t t = 1 1-ζ P us t + 1 ζ W t and W t ≥ 0, thus P s t ≥ 1 1-ζ P us t From (8) this gives

  Sketch of proof on the relation between V us t and V s t-1 . Using the fact that ỹt, θ us t ≤ δ t and δ v t ≤ δ t it can be

				θ s T t	P -1 t	θ s t
	with θ s t = θ * -θ s t .
	• Sketch of proof on the relation between V s t and V us t . V s t can be rewritten as
	V s t = θ us t -w t	T	1 1 -ζ	P us t +	1 ζ	W t	-1	θ us t -w t
	We have θ us T t	P us -1 t	θ us t ≤ 1 and w T t W -1 t w t ≤ 1 this
	gives
	V s t ≤ (1 -ζ )V us t + ζ w T t W -1 t w t	(28)
	shown that
			V us t ≤ λV s t-1

•


Mathieu Pouliquen
email: mathieu.pouliquen@unicaen.fr

Fouad Giri

Olivier Gehan

Eric Pigeon

Miloud Frikel

Boubekeur Targui

Subspace Identification of Hammerstein Systems with Nonparametric Input Backlash and Switch Nonlinearities

In this paper, a subspace identification algorithm for a class of Hammerstein systems is developed. We consider dynamical systems subject to input backlash or switch nonlinearities. The idea is to use a specific input signal allowing the estimation of the nonlinear part and the estimation of a state space model for the linear part. The identification algorithm is a subspace type algorithm. A simulation example is given to illustrate the performances of the present method.

I. INTRODUCTION

The Hammerstein systems are the nonlinear systems consisting of a nonlinearity block followed by a linear dynamic system. The identification of such systems has received much interest for the last decades as evidenced by the many following references: [START_REF] Narendra | An Iterative Method for the Identification of Nonlinear Systems Using a Hammerstein Model[END_REF], [START_REF] Chang | Luus A Noniterative Method for Identification Using Hammerstein Model[END_REF], [START_REF] Bai | Identification of linear systems with hard input nonlinearities of known structure[END_REF], [START_REF] Chen | Strong Consistency of Recursive Identification for Hammerstein Systems With Discontinuous Piecewise-Linear Memoryless Block[END_REF], [START_REF] Zhao | Recursive Identification for Hammerstein System With ARX Subsystem[END_REF], [START_REF] Giri | Identification of Hammerstein systems in presence of hysteresis-backlash and hysteresisrelay nonlinearities[END_REF], [START_REF] Wang | Revisiting Hammerstein system identification through the Two-Stage Algorithm for bilinear parameter estimation[END_REF], [START_REF] Zhao | Adaptive tracking and recursive identification for Hammerstein systems[END_REF], [START_REF] Bai | Convergence of the iterative algorithm for a general Hammerstein system identification[END_REF], [START_REF] Giri | Block Oriented Nonlinear System Identification[END_REF], [START_REF] Rochdi | Identification of block-oriented systems in the presence of nonparametric input nonlinearities of switch and backlash types[END_REF], [START_REF] Voros | Modeling and identification of systems with backlash[END_REF], [START_REF] Chen | Recursive Identification for MIMO Hammerstein Systems[END_REF], [START_REF] Giri | Parameter identification of Hammerstein systems containing backlash operators with arbitrary-shape parametric borders[END_REF], [START_REF] Le | Recursive Identification of Hammerstein Systems with application to Electrically Stimulated Muscle[END_REF], [START_REF] Han | Callafon Hammerstein system identification using nuclear norm minimization[END_REF] and the four sessions devoted to block oriented nonlinear identification at SYSID 2012. This interest is justify by the fact that these systems take into account nonlinearities commonly encountered in practice, generated by the technological limitations necessary for the proper functioning of the system (saturation, the limit stops, etc.). In this paper we consider the identification in presence of some particular input nonlinearities: backlash and switch nonlinearities.

Few contributions deal with the identification of system with input backlash or switch nonlinearity. The approaches developed in [START_REF] Bai | Identification of linear systems with hard input nonlinearities of known structure[END_REF], [START_REF] Cerone | Regruto Bounding the parameters of linear systems with input backlash[END_REF], [START_REF] Giri | Identification of Hammerstein systems in presence of hysteresis-backlash and hysteresisrelay nonlinearities[END_REF] and [START_REF] Cerone | Bounded-error identification of linear systems with input and outmut backlash[END_REF] are devoted to nonlinearities bordered by straight lines. The method introduced in [START_REF] Rochdi | Identification of block-oriented systems in the presence of nonparametric input nonlinearities of switch and backlash types[END_REF] is adapted to nonlinearities with more general borders however the proposed method is a two experiments method. Moreover it uses Least Squares algorithm coupled to an overparameterization, which leads to the need for a high number of data. An other method is also depicted in [START_REF] Giri | Parameter identification of Hammerstein systems containing backlash operators with arbitrary-shape parametric borders[END_REF]. Nevertheless this method is only adapted to backlash nonlinearity and it is a two stages method using two independent experiments.

In this paper, we propose an approach which alleviates some of the issues present in the previous methods. First, our method is devoted to backlash and switch nonlinearities with general borders. Second, our method is a single stage method which allows to estimate both the linear subsystem and the nonlinearity. Let us notice that the borders of the nonlinearities are nonparametric, they will be piecewise parameterized as in [START_REF] Van Pelt | nonlinear system identication using Hammerstein and nonlinear feedback models with piecewise linear static maps[END_REF].

The proposed identification method is based on two key idea: (1) A specific input sequence is design in order to allow a linear parameterization of the identification problem and in order to satisfy the persistent excitation condition. ( 2) Subspace identification methods are used for the estimation. The use of such methods is mainly motivated by a set of interesting properties: the simplicity, the intrinsic numerical robustness and their straightforward application to MIMO systems ( [START_REF] Van Overschee | Subspace identification for linear systems. Theory, implementation applications[END_REF], [START_REF] Verhaegen | Subspace model identification part 3: analysis of the ordinary output error state space model identification algorithm[END_REF], [START_REF] Verhaegen | Identification of the deterministic part of MIMO state space models given in innovations form from input output data[END_REF], [START_REF] Favoreel | Subspace state space system identification for industrial processes[END_REF], [START_REF] Katayama | Subspace methods for system identification[END_REF], [START_REF] Qin | A novel subspace identification approach with enforced causal models[END_REF], [START_REF] Qin | An overview of subspace identification[END_REF], [START_REF] Chiuso | The role of vector auto regressive modeling in predictor based subspace identification[END_REF]). Some subspace methods adapted to Hammerstein systems have been introduced ([15], [START_REF] Goethals | Subspace intersection identification of Hammerstein-Wiener systems[END_REF], [START_REF] Bako | Recursive Subspace Identification of Hammerstein Models Based on Least Squares Support Vector Machines[END_REF], [START_REF] Yana | Subspace identification methods for Hammerstein systems: rank constraint and dimension problem[END_REF], [START_REF] Van Wingerden | Subspace identification of Hammerstein-Wiener systems operating in closed loop[END_REF], [START_REF] Jalaleddini | An identification algorithm for Hammerstein systems using subspace method[END_REF], [START_REF] Naitali | Persistent Excitation by Deterministic Signals in Subspace MISO Hammerstein System Identification[END_REF]) but, to the best of our knowledge, these methods aren't adapted to the case of backlash or switch nonlinearity.

The outline of this paper is as follows: in section II the model structure is presented and the specific input signal is designed. In section III the subspace identification algorithm is applied. Some simulation results are given in section IV. Finally, section V concludes the paper.

II. PROBLEM STATEMENT AND PARAMETRIZATION

In this section we first state the identification problem. Then we described the design of the excitation signal. The association of that signal with a parametrization of the nonlinearity yields a linear identification problem.

A. Identification problem statement

Consider the identification of the Hammerstein system shown in Fig. 1 where u(t) and y(t) are respectively the input and the output of the plant which is expressed as:

y(t) = G(q)H[u(t)] + v(t) (1)
G(q) is the proper transfer function of the deterministic linear part in the Hammerstein system and v(t) is an additive noise on the output. v(t) is wide sense stationary zero mean process uncorrelated with u(t). The nonlinear part in the Hammerstein system is described by H [.]. This function can be a switch or a backlash element. Its borders are characterized by the descendant and ascendent continuous borders H d [.] and H a [.]. These two nonlinearities are illustrated on Fig. 2 and Fig. 3. The identification problem treated in this paper is stated as: based on the input-output measurements, find the minimal order n of the system, a realization of state space matrices (A, B,C, D) for G(q) and a nonparametric description of the descendant and ascendent borders H d [.] and H a [.]. The identification of the nonlinearity will be performed on a working interval: u(t) ∈ [u min ; u max ].

✲ y(t) v(t) u(t) ✲ H[.] ✲ G(q) ✲ ❄ + + Fig. 1. Considered Hammerstein system ✻ ✲ u(t) H[u(t)] ✒ ✛ ✛ ✛ ✛ ✠ H d [.] H a [.]
So as to be able to solve the identification problem we make some assumptions throughout the paper: A.1 In the case of the switch nonlinearity we have

H a [u min ] = H d [u min ]; A.2 The static gain of G(q) is unitary: G(1) = 1;

Remark 1:

The assumption [A.2] is necessary so as to satisfy the identifiability condition. Without this assumption we are not able to distinguish the solution G(q)H[u(t)] from the solution (aG(q))( 1a H[u(t)]) with a = 0.

B. Design of the excitation signal and parametrization of the identification problem

In this subsection first we describe the suggested excitation signal, then we use it with a piecewise parameterization of the nonlinearity.

The design of the input signal u(t) is the corner stone of the identification method proposed in this paper. A major difficulty is to select an excitation satisfying the persistent excitation condition and allowing a characterization of the

✻ ✲ u(t) H[u(t)] H a [.] H d [.] ✒ ✠ ✻ ✻ ✻ Fig. 3. Switch nonlinearity
Input design procedure

• The initial value is u(0) = u min
• For t ≥ 1, follow the following cycle: step 1 u(t) equal to a random value in ]u min ; u max [ (we take the ascending path) then t = t + 1; step 2 u(t) = u max (we continue the ascending path until u max ) then t = t + 1; step 3 u(t) equal to a random value in ]u min ; u max [ (we take the descending path) then t = t + 1; step 4 u(t) = u min (we continue the descending path until u min ) then t = t + 1 and go back to step 1. ascendant and descendant paths. Such an excitation is illustrated on Fig. 4 and table I presents its design procedure.

This excitation signal has two key properties. Firstly it produces persistent excitation on the input of the system thanks to the random value one sample on two. Secondly u(t) browses continuously in their entirety the ascending and descending paths. Such an excitation ensures cycles described on Figs. 2 and 3 to be driven always in the same manner between u min and u max . If u min and u max are constant (we assume they are), these paths are unique for the backlash nonlinearity and the switch nonlinearity. Let us denote respectively P a [.] and P d [.] the ascending and descending path i.e. the extension of the borders H a [.] and H d [.] in the following manner:

• if u(t) is increasing then H[u(t)] = P a [u(t)]; • if u(t) is decreasing then H[u(t)] = P d [u(t)].
At any time t it is possible to know which path uses H[u(t)]. To this end, let define a direction signal δ (t) which

[1] [2] [1] [2] [1]
Fig. 4. The excitation signal is divided into two phases: the first [START_REF] Bai | Identification of linear systems with hard input nonlinearities of known structure[END_REF] in which we take the ascending path, the second [START_REF] Bai | Convergence of the iterative algorithm for a general Hammerstein system identification[END_REF] for which we take the descending path.

indicates if u(t) is increasing or decreasing:

δ (t) =    1 if u(t) > u(t -1) δ (t -1) if u(t) = u(t -1) 0 if u(t) < u(t -1)
From this definition, H[u(t)] can be formulated as follows

H[u(t)] = δ (t)P a [u(t)] + (1 -δ (t))P d [u(t)] (2)
So as to characterize the functions P a [u(t)] and P d [u(t)], we approximate them by continuous piecewise linear functions as in [START_REF] Van Pelt | nonlinear system identication using Hammerstein and nonlinear feedback models with piecewise linear static maps[END_REF]. These approximations consist in a sum of functions {P a i (.)} for the ascending path and a sum of functions {P d i (.)} for the descending path in the following manner:

P a [u(t)]=µ a 0 +µ a 1 P a 1 (u(t))+µ a 2 P a 2 (u(t))+•••+µ a n P P a n P (u(t)) P d [u(t)]=µ d 0 +µ d 1 P d 1 (u(t))+µ d 2 P d 2 (u(t))+•••+µ d n P P d n P (u(t)) (3)
with µ a 0 = P a (u min ), µ d 0 = P d (u min ) and

P a i (u(t)) =    0 if u(t) ≤ n a i-1 u(t) -n a i-1 if n a i-1 ≤ u(t) ≤ n a i n a i -n a i-1
if n a i ≤ u(t) Fig. 5 illustrates such function P a i (u(t)). n P + 1 represents the number of nodes between u min and u max . n P + 1 is a user defined index which depends on the application considered as the distributions of the nodes {n a i } and {n d i }. Theses nodes can be arbitrarily selected by the user on the interval [u min ; u max ].

Remark 3: For ease of notations we have chosen the same number of knots for the ascending path and the descending path, however they could be different.

From [START_REF] Bai | Convergence of the iterative algorithm for a general Hammerstein system identification[END_REF] and (3) we get

H[u(t)]=µ a 0 δ (t)+µ d 0 (1-δ (t))+∑ n P i=1 µ a i P a i (u(t))δ (t)+∑ n P i=1 µ d i P d i (u(t))(1-δ (t)) ✻ ✲ u(t) P a i (u(t)) • • • n a i n a i-1 • • • Fig. 5.
The function P a i (u(t)) for the approximation on the segment

[n a i-1 ;n a i ]
and according assumption A.1 whatever the nonlinearity we have µ a 0 = µ d 0 , this gives

H[u(t)] = µ a 0 + n P ∑ i=1 µ a i P a i (u(t))δ (t) + n P ∑ i=1 µ d i P d i (u(t))(1 -δ (t))
and then

H[u(t)] = µz(t) + µ a 0 ( 4
)
with

z(t) =               P a 1 (u(t))δ (t) P a 2 (u(t))δ (t) . . . P a n P (u(t))δ (t) P d 1 (u(t))(1 -δ (t)) P d 2 (u(t))(1 -δ (t))
. . .

P d n P (u(t))(1 -δ (t))               and µ = µ a 1 µ a 2 • • • µ a n P µ d 1 µ d 2 • • • µ d n P
From assumptions A.2 and making use previous notations in (1) yields

y(t) -µ a 0 = G(q)µz(t) + v(t) (5)
Let consider here a minimal state space model for G(q), this gives

x(t + 1) = Ax(t) + Bµz(t) y(t) -µ a 0 = Cx(t) + Dµz(t) + v(t) (6)
where x(t) ∈ R n is the state vector with n the minimal order,

A ∈ R n×n , B ∈ R n , C ∈ R n and D ∈ R. (A,C
) is assumed to be observable and (A, B) is assumed to be controllable.

From [START_REF] Chang | Luus A Noniterative Method for Identification Using Hammerstein Model[END_REF] it is possible to define an over-parameterized state space model in the following manner:

x(t + 1) = Ax(t) + Bz(t) y(t) -µ a 0 = Cx(t) + Dz(t) + v(t) (7)
with B = Bµ and D = Dµ. This corresponds to a linear state space model with the input z(t) ∈ R 2n P and the output y(t)µ a 0 . In the next section this over-parameterized model is used in the identification step of the approach.

III. SUBSPACE IDENTIFICATION ALGORITHM

As in [START_REF] Yana | Subspace identification methods for Hammerstein systems: rank constraint and dimension problem[END_REF], we use here standard subspace algorithms for the identification of the state space matrix of [START_REF] Chen | Strong Consistency of Recursive Identification for Hammerstein Systems With Discontinuous Piecewise-Linear Memoryless Block[END_REF]. Input and output are respectively z(t) and y(t) -µ a 0 , consequently µ a 0 is supposed to be known. To this end we recommend the application of an input sequence equal to u min for a few moments before the application of the excitation sequence described in subsection II-B. The mean of the corresponding output gives µ a 0 . In the following, let us denote ȳ(t) = y(t)µ a 0 .

Standard subspace algorithms split the available data into block matrices: output Hankel matrices Y p and Y f are defined as

Y p Y f =             ȳ(t-p)
ȳ(t-p+1) . . . ȳ(t-p+ j-1) ȳ(t-p+1) ȳ(t-p+2) . . . . . .

ȳ(t+ j-1) ȳ(t+1) ȳ (t+2)
. . .

ȳ(t+ j) . . . . . . . . . . . . ȳ(t+ f -1) ȳ(t+ f ) . . . ȳ(t+ f + j-2)            
Similar definitions hold for V f , Z f and Z p . Subscripts f and p respectively stand for "future" and "past" and are user-defined indexes. j is the number of columns in the Hankel matrices. In the subspace identification framework it is commonly assumed that there are long time series available i.e. j ≫ 1.

Using these Hankel matrices, the over-parameterized state space model ( 7) can be represented by the following basic subspace equations:

X t = A p X t-p + ∆ p Z p Y f = Γ f X t + Ξ f Z f + V f (8)
where X t-p and X t are respectively the past and future state vectors of the plant:

X t-p = x(t -p) x(t -p + 1) . . . x(t -p + j -1) X t = x(t) x(t + 1) . . . x(t + j -1) Γ f =     C CA • • • CA f -1     and ∆ p = A p-1 B • • • AB B
are respectively the extended observability matrix and reversed extended controllability matrix of the model [START_REF] Chen | Strong Consistency of Recursive Identification for Hammerstein Systems With Discontinuous Piecewise-Linear Memoryless Block[END_REF] . Ξ f is defined by

Ξ f =     D 0 • • • 0 CB D 0 • • • • • • • • • • • • • • • CA f -2 B • • • CB D    
and it corresponds to a lower triangular Toeplitz matrix containing the first f Markov parameters of [START_REF] Chen | Strong Consistency of Recursive Identification for Hammerstein Systems With Discontinuous Piecewise-Linear Memoryless Block[END_REF].

From [START_REF] Chen | Recursive Identification for MIMO Hammerstein Systems[END_REF], Y f can be rewritten as

Y f = L Z f Z f + L Z p Z p + Γ f A p X t-p + V f with L Z p = Γ f ∆ p L Z f = Ξ f
This becomes

Y f = LΦ + Γ f A p X t-p + V f (9) with L = L Z p L Z f Φ = Z p Z f ( 10
)
Most of classical subspace identification methods are based on the estimation of the extended observability matrix Γ f ([26]). This can be done making use the projection of

Y f Φ T ΦΦ T -1 Φ onto the orthogonal complement of the row space of Z f : O f = Y f Φ T ΦΦ T -1 ΦΠ ⊥ Z f (11) with Π ⊥ Z f = I j -Z T f (Z f Z T f ) -1 Z f .
In [START_REF] Van Overschee | Subspace identification for linear systems. Theory, implementation applications[END_REF] it is shown that if f , p > n, if the input z(t) is uncorrelated with the additive noise v(t) and if z(t) is persistently exciting of order ( f + p)2n P then the singular values decomposition (SVD) of lim j-→∞ O f is given by lim

j-→∞ O f = U 1 U 2 S 1 0 0 0 V T 1 V T 2 (12) lim j-→∞ O f = U 1 S 1 V T 1
where S 1 ∈ R n×n . It follows that the (unknown) order of the system is equal to the number of singular values different from zero and the observability matrix Γ f can be taken equal to

Γ f = U 1 T ( 13
)
where T is a similarity transformation.

Remark 4: Some weighting matrices are often added in [START_REF] Giri | Block Oriented Nonlinear System Identification[END_REF] in order to include some well known subspace methods such that N4SID, PO-MOESP, PI-MOESP, CVA, IVM, etc.. Here [START_REF] Giri | Block Oriented Nonlinear System Identification[END_REF] corresponds to the PI-MOESP method ([29], [START_REF] Verhaegen | Identification of the deterministic part of MIMO state space models given in innovations form from input output data[END_REF]).

Once the order and the extended observability matrix estimated (with ( 13)), we have to find the state matrices (A, B,C, D) and the parameters {µ a i } and {µ d i }. The following three steps procedure can be considered: 1) From (13), the system matrices A and C are obtained as follows 1 C = Γ f (1, :)

A = Γ T f Γ f -1 Γ T f Γ f with Γ f = Γ f (1 : f -1, :) and Γ f = Γ f (2 : f , :).
2) B and D are determined using the following equation

Y f Φ T ΦΦ T -1 Φ = LΦ + Γ f A p X t-p Φ T ΦΦ T -1 Φ which gives Γ f ⊥ Y f Φ T ΦΦ T -1 Φ = Γ f ⊥ Ξ f (B, D)Z f
Then the procedure given in [START_REF] Van Overschee | Subspace identification for linear systems. Theory, implementation applications[END_REF] is used to obtain B and D.

3) B and µ are extracted from B using a second SVD:

B = u ′ 1 U ′ 2 s ′ 1 0 0 0 v ′ T 1 V ′ T 2 (14) This gives B = γu ′ 1 s ′ 1 µ 1 µ 2 • • • µ n P = 1 γ v ′ T 1
γ is a parameter adjusted so as to ensure a unitary static gain (assumption A.2).

Remark 5: To improve the numerical robustness, once A and C are known, a method is proposed in [START_REF] Yana | Subspace identification methods for Hammerstein systems: rank constraint and dimension problem[END_REF] in order to estimate the matrices B and µ separately. This method could be used here.

Remark 6: In this paper only the deterministic part of the plant is identified. It is possible to consider a state space model with noise model which can be rewritten in an innovation form as

x(t + 1) = Ax(t) + BH[u(t)] + Ke(t) y(t) = Cx(t) + DH[u(t)] + e(t)
The stochastic part is described by the matrix K ∈ R n and the covariance σ 2 e = E e(t) 2 . These parameters could easily be estimated using a method proposed in [START_REF] Knudsen | Subspace identification -reducing uncertainty on the stochastic part[END_REF].

IV. SIMULATION RESULTS

To investigate the interest of the method, two numerical example are proposed. On these examples we consider a plant whose linear part is described by the following two order transfert function:

G(q) = q -1 0.1994 -0.1804q -1 1 -1.8858q -1 + 0.9048q -2 1 MATLAB c notations v(t) is such that v(t) = G v (q)e(t)
where e(t) is a gaussian white noise and the noise model is given by

G v (q) = 1 + 0.5q -1 1 -0.85q -1
The identification method is applied on the basis of Monte Carlo simulation of 100 experiments. The variance of e(t) is chosen such that signal to noise ratio is SNR = 10dB. The number of columns in Hankel matrices is j = 4000.

A. Case 1: a switch nonlinearity

In this first case the nonlinear part H[.] is a switch nonlinearity characterized by the following functions H a [.] and H d [.]:

H d [u(t)] =    -1 if u(t) ≤ -2 arctan(u(t) + 1) -H 0 if -2 < u(t) ≤ 2 1 if 2 < u(t)
and

H a [u(t)] =    -1 if u(t) ≤ -2 arctan(u(t) -1) + H 0 if -2 < u(t) ≤ 2 1 if 2 < u(t)
with H 0 = 1 2 (arctan(3)arctan( 1)). These functions are approximated on the interval [-3; 3] and we choose a low number of nodes n P = 7. These nodes are uniformly distributed on [-3; 3] and they are the same for the descendant and ascendent paths.

The estimated poles are shown in Fig. 6, the crosses indicate the true pole locations. This figure shows that the proposed method gives good estimation: no bias and low standard-deviation. Such good estimation quality is confirmed by Fig. 7 which shows the frequency response of the true system and the average of the estimated frequency responses. The estimated average nonlinearity is depicted on Fig. 8. The same conclusion as for Fig. 7 holds.

B. Case 2: a backlash nonlinearity

In this second case the nonlinear part H[.] is a backlash nonlinearity characterized by the following functions H a [.] and H d [.]:

H d [u(t)] = arctan(u(t) + 1) and H a [u(t)] = arctan(u(t) -1)
These functions are again approximated on the interval [-3; 3] with the same number of nodes: n P = 7.

The results are given in Figs. 9, 10 and 11. They confirm the previous conclusion: the proposed method works well. It allows an approximation of the nonlinearity and an unbiased estimation of the linear part. In this paper, an identification method is proposed for the identification of linear systems with input switch or backlash nonlinearity. The key step of the algorithm is the design of a specific excitation signal allowing a suitable parametrization of the identification problem. This choice allows the use of subspace identification methods and then the estimation of a model for both the linear part and the non linear part. A numerical example shows the effectiveness of the proposed algorithm. In this paper only the open loop identification problem is considered, we believe that the method developed here can be extended to the closed loop identification problem.

Fig. 2 .

 2 Fig. 2. Backlash nonlinearity

Remark 2 :

 2 In order to characterize homogeneously functions H a[.] and H d[.] on [u min ; u max ], we recommend choosing a random value uniformly distributed over this interval.

  .............................................. .................................................

Fig. 6 .

 6 Fig.6. Switch nonlinearity: the eigenvalues of estimated A matrix: * estimated pole, + system pole.

Fig. 7 .Fig. 8 .Fig. 9 .

 789 Fig. 7. Switch nonlinearity: the Bode magnitude plot of the system and of the estimated average model

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Backlash nonlinearity: the Bode magnitude plot of the system and of the estimated average model

TABLE I

 I