S Berthet

M Leriche

J.-P Pinty

J Cuesta
email: cuesta@lmd.polytechnique.fr

G Pigeon

SCAVENGING OF AEROSOL PARTICLES BY RAIN IN A CLOUD RESOLVING MODEL

Keywords: Below Cloud Scavenging, Collision efficiency, Gauss quadratures, aerosol transport simulations

We describe a below-cloud scavenging module of aerosol particles by raindrops for use in a three-dimensional mesoscale cloud resolving model. The rate of particle removal is computed by integrating the scavenging efficiency over the aerosol particle and the drop size distributions. Here the numerical integration is performed accurately with a Gauss quadrature algorithm. The efficiency of the scavenging module is partially confirmed with experimental data. More interestingly, it is illustrated by two numerical experiments: the simulation of a forced convective circulation in a tropical cloudy boundary layer and a two-dimensional simulation of an African squall line. The results show a very selective wet removal of the aerosol particles which clearly depends on the mode radius, the width and the vertical profile of concentration.

Furthermore, the squall line case shows the importance of resolving internal

Introduction

Wet deposition is the most important sink of aerosol particles (AP) in the troposphere [START_REF] Pruppacher | Microphysics of Clouds and Precipitation[END_REF]. It involves complex microphysical interactions between AP and hydrometeors which are categorized as nucleation scavenging and impaction scavenging processes. The nucleation scavenging leads to the formation of cloud droplets and ice crystals on cloud condensation nuclei and on ice forming nuclei by heterogeneous nucleation, respectively. The impaction scavenging applies to all the dynamical and physical processes leading to the capture of AP by cloud droplets, ice crystals and mostly by falling hydrometeors. It includes Brownian diffusion, interception, inertial impaction, thermophoresis, diffusiophoresis, airflow turbulence and electrostatic attraction [START_REF] Pruppacher | Microphysics of Clouds and Precipitation[END_REF][START_REF] Andronache | Estimated variability of below-cloud aerosol removal by rainfall for observed aerosol size distributions[END_REF].

The nucleation scavenging, by essence an in-cloud mechanism, dominates the impaction scavenging which is efficient once big drops are precipitating at cloud base. So the impaction scavenging is often referred to as a below-cloud scavenging (BCS) process but this does not rule out the in-cloud capture of interstitial AP by raindrops.

This study focuses on the role of BCS of AP by the raindrops thus ignoring similar effects by the snowflakes and other ice particles (see the discussion in section 4.2). We refer to the recent work of [START_REF] Feng | A size-resolved model for below-cloud scavenging of aerosols by snowfalls[END_REF] for the scavenging of AP by snow particles of different shapes. This work is part of a larger project, which aims at describing all the processes leading to the wet deposition of polydisperse AP in a three-dimensional (3D) cloud resolving model.

According to [START_REF] Pruppacher | Microphysics of Clouds and Precipitation[END_REF], high values of BCS are found both for ultrafine particles (diameters less than 0.01 µm) and for coarse particles (diameters larger than 2 µm). Interestingly, for the intermediate range of particle size referred as the "Greenfield gap", the BCS efficiency decreases by several orders of magnitude, thus making the particles less depleted and so available for in-cloud nucleation scavenging in the upper levels.

The explicit computation of BCS rates fits well with the present capabilities of cloud resolving/mesoscale models for two reasons at least. First, these models are able to produce 3D rain fields in a consistent way (mixing ratios, and even drop concentrations depending on the microphysical scheme) and second, they resolve in detail the air flow patterns which determine the AP exposure duration of air parcels to BCS removal during rainfall events. This contrasts with estimates of BCS in global or chemistry transport models where crude approximations (see a review by [START_REF] Sportisse | A review of parameterizations for modelling dry deposition and scavenging of radionuclides[END_REF]) are made about the size distribution of the raindrops or the rainfall rate estimates [START_REF] Seinfeld | Atmospheric Chemistry and Physics[END_REF][START_REF] Jacobson | Development of mixed-phase clouds from multiple aerosol size distributions and the effect of the clouds on aerosol removal[END_REF][START_REF] Loosmore | Precipitation scavenging of atmospheric aerosols for emergency response applications: testing an updated model with new real-time data[END_REF]Tost et al., 2006;[START_REF] Henzing | A parametrization of size resolved below cloud scavenging of aerosols by rain[END_REF][START_REF] Feng | A 3-mode parametrization of below-cloud scavenging of aerosols for use in atmospheric dispersion models[END_REF] to cite a few. In return, considering polydisperse AP and raindrops at each gridpoint of 3D domains of simulation implies double integrations of the BCS rates over the two size distributions [START_REF] Mircea | Precipitation scavenging coefficient: influence of measured aerosol and raindrop size distributions[END_REF][START_REF] Andronache | Estimated variability of below-cloud aerosol removal by rainfall for observed aerosol size distributions[END_REF][START_REF] Jung | A moment model for simulating raindrop scavenging of aerosols[END_REF]. For this reason and because of the vanishing part of the distributions, a special attention needs to be drawn upon computation efficiency and accuracy, as stressed by [START_REF] Sportisse | A review of parameterizations for modelling dry deposition and scavenging of radionuclides[END_REF].

This study presents a flexible BCS module and illustrates some applications in the 3D mesoscale/cloud resolving model MesoNH [START_REF] Lafore | The Meso-NH atmospheric simulation system. Part I: adiabatic formulation and control simulations[END_REF]. First, the BCS module has been developed based on the widely used collision efficiency parameterization of [START_REF] Slinn | A reevaluation of the role of thermophoresis as a mechanism of in-and below-cloud scavenging[END_REF]. The original module employs a quadrature method to integrate over the raindrop and the AP size distributions. Then zero-dimensional tests of the module have been done on a set of rain rates and AP concentration data from the COPS experiment (Wulfmeyer and coll., 2008), which took place in the Vosges-Black Forest area in July 2007 (see http://www.cops2007.de/). Finally, the module has been implemented and coupled to the microphysical schemes of MesoNH.

Two simulations have been performed to illustrate the BCS module and its potential interactions with the transport. We studied an idealized case of a shallow tropical rainband in a 2D kinematic framework from the HaRP campaign [START_REF] Cohard | A comprehensive two-moment warm microphysical bulk scheme. II: 2D experiments with a non-hydrostatic model[END_REF] and a squall line case from the COPT81 experiment in West Africa [START_REF] Lafore | A numerical investigation of the organization and interaction of the convective and stratiform regions of tropical squall lines[END_REF]. For both simulations, sensitivity tests are performed on the characteristics of the initial size distributions of AP and on the initial vertical profile of AP spectra.

DEVELOPMENT OF THE BCS MODULE

Total number BCS rate

A BCS module has been developed to describe the wet removal of polydisperse AP by rain precipitation inside the cloud and especially below the cloud base in the three-dimensional model MesoNH.

BCS of AP is related to the AP size distribution but also to the raindrop size distribution through the BCS coefficient. Here a difficulty arises due to the complex mathematical expression of the BCS coefficient which must be integrated over wide distributions, up to infinity for analytical, log-normal or gamma laws. In order to simplify crudely the computation of a BCS rate, Tost et al. (2006) and [START_REF] Loosmore | Precipitation scavenging of atmospheric aerosols for emergency response applications: testing an updated model with new real-time data[END_REF] assume a monodisperse distribution of the raindrops (or equivalently an empirical "representative" raindrop diameter as shown in Fig. 4 of [START_REF] Sportisse | A review of parameterizations for modelling dry deposition and scavenging of radionuclides[END_REF]) to determine the BCS coefficient γ(d p ) of AP of size d p and the falling raindrops. A less empirical method used by [START_REF] Andronache | Estimated variability of below-cloud aerosol removal by rainfall for observed aerosol size distributions[END_REF] and by [START_REF] Feng | A 3-mode parametrization of below-cloud scavenging of aerosols for use in atmospheric dispersion models[END_REF] takes explicitly into account the spectral nature of the raindrop distribution by performing a numerical integration up to a large raindrop diameter. In order to keep an accurate estimate of BCS effects in MesoNH, we choose to compute BCS coefficients by means of a Gauss quadrature method [START_REF] Press | Numerical Recipes in FORTRAN: The Art of Scientific Computing[END_REF]. Infinite integrals for BCS coefficients are computed as a weighted sum for a small number (here ∼ 20) of optimized abscissas (each abscissa is deduced from a particle or drop diameter by a change of variable). Then, the final BCS rate of AP is obtained by integrating the BCS coefficient γ(d p ) over the AP size distribution using another Gauss quadrature, adapted to the AP size distribution.

The AP are represented by the sum of several log-normal distributions n pi (d p )dd p depending on the number of modes:

dN p = l i=1 n pi (d p )dd p = l i=1 N pi √ 2πd p lnσ i e -ln(d p /d pi ) √ 2lnσ i 2 dd p
where l is the number of modes with index i, d p is the AP diameter, N pi , the number concentration, σ i , the geometric standard deviation of the lognormal distribution, and d pi is the median diameter. For the time being, we use a single-moment scheme for the AP in MesoNH (the number concentration N p of each mode is a prognostic variable), so σ i and d pi are kept constant, depending on the AP types 1 . Traditionally three modes (l = 3) are superimposed: the Aitken nuclei mode, the accumulation mode and the coarse mode.

The raindrop size distribution n R (D d ) in MesoNH is modeled by a gen-

eralized gamma distribution n R (D d ) = (α/Γ(ν)) λ αν R D αν-1 d exp(-λ R D d ).
In this study, n R (D d ) is reduced to the classical Marshall-Palmer law by taking α = ν = 1 for the two shape parameters:

n R (D d )dD d = N 0 e -λ R D d dD d (1)
where N 0 = 8 × 10 -3 m -3 mm -1 , λ R is the slope parameter and D d the drop diameter.

The method used to calculate the BCS coefficient follows the concept of the collision efficiency E between an AP and a raindrop [START_REF] Slinn | A reevaluation of the role of thermophoresis as a mechanism of in-and below-cloud scavenging[END_REF][START_REF] Pruppacher | Microphysics of Clouds and Precipitation[END_REF][START_REF] Seinfeld | Atmospheric Chemistry and Physics[END_REF]. studies about BCS are based on this concept [START_REF] Andronache | Estimated variability of below-cloud aerosol removal by rainfall for observed aerosol size distributions[END_REF][START_REF] Loosmore | Precipitation scavenging of atmospheric aerosols for emergency response applications: testing an updated model with new real-time data[END_REF][START_REF] Henzing | A parametrization of size resolved below cloud scavenging of aerosols by rain[END_REF]Tost et al., 2006;[START_REF] Feng | A 3-mode parametrization of below-cloud scavenging of aerosols for use in atmospheric dispersion models[END_REF].

E
There is no available theoretical solution of the Navier-Stokes equation for the prediction of E. Therefore [START_REF] Slinn | A reevaluation of the role of thermophoresis as a mechanism of in-and below-cloud scavenging[END_REF] has formulated a semi-empirical expression of the collision efficiency E, taking into account three of the most understood collection processes: the Brownian diffusion, the interception and the inertial impaction. Some studies (Slinn and Hales, 1971;Wang et al., 1978;Tinsley et al., 2005) describe the relative impact of other processes such as thermo-, diffusiophoresis and electric forces upon E. The inclusion of these processes means to estimate the electric charge of each particle and to determine the local relative humidity along particle trajectories, complicating heavily the computation of E. In general these processes are ignored, as in this study. After non-dimensionalizing the equation of motion for the air, the AP and the raindrops, [START_REF] Slinn | A reevaluation of the role of thermophoresis as a mechanism of in-and below-cloud scavenging[END_REF] found an expression of E that depends on five dimensionless parameters :

Re = D d U t (D d ) ρ a 2 µ a , Sc = µ a ρ a D , St = 2 U t (D d ) τ a D d φ = d p D d , ω = µ w µ a , and St * = 1.2 + (1/12ln(1 + Re)) 1 + ln(1 + Re)
with symbols defined in Table 1. [START_REF] Seinfeld | Atmospheric Chemistry and Physics[END_REF] or [START_REF] Sportisse | A review of parameterizations for modelling dry deposition and scavenging of radionuclides[END_REF] provide expressions for D and τ a which depend on d p . In MesoNH, the terminal velocity of raindrops of diameter D d is approximated by a simple power law dependence in diameter that includes the air density effect expressed by [START_REF] Foote | Terminal velocity of raindrops[END_REF]:

U t (D d ) = a D b d ( ρ a0 ρ a ) 0.4
Finally the analytical Slinn's expression of E that fits experimental data is expanded as:

E(D d , d p ) = 4 Re Sc 1 + 0.4Re 1/2 Sc 1/3 + 0.16Re 1/2 Sc 1/2 +4φ ω -1 + (1 + 2Re 1/2 )φ + St -St * St -St * + 0.667 3/2 ρ p ρ w 1/2
(2)

The analytical expression of the collision efficiency E contains three terms corresponding to distinct physical contributions depending on the AP diameter d p :

1. The first term is the contribution of the Brownian diffusion. It denotes the random motion of small particles under the influence of the surrounding fluid molecules. It dominates for AP with d p < 0.01µm.

2. The second term is the collection by interception. The particles follow the turbulent air motion in the wake of the drops. They are collected when their trajectory is close enough to the drop. This process affects AP with diameter between 0.01 and 2 µm including AP in the "Greenfield gap".

3. The third term refers to the inertial impaction. It occurs when the particles do not follow the mean air flow because of their inertia that takes them off the streamlines. This term is efficient for large particles (d p > 2µm) with a restriction to cases where St > St * in Eq. 2 [START_REF] Feng | A 3-mode parametrization of below-cloud scavenging of aerosols for use in atmospheric dispersion models[END_REF].

Wet deposition by BCS is a 1 st order decay process so that,

∂ψ(d p , t) ∂t = -γ(d p )ψ(d p , t) (3)
where ψ(d p , t) is the AP size distribution in number or mass. γ(d p ) is the BCS coefficient for particles of diameter d p . It results from the integration over the raindrop size distribution of the collision efficiency E times the volume swept by a falling raindrop D d per unit of time,

γ(d p ) = ∞ 0 π 4 D 2 d U t (D d )E(D d , d p )n R (D d )dD d (4)
A Gauss-Laguerre algorithm [START_REF] Press | Numerical Recipes in FORTRAN: The Art of Scientific Computing[END_REF], a variant of the Gauss quadrature technique, is used for the first time to integrate the above expression. The method is based on a n-orthogonal polynomial set associated to the weighting function, here the size distribution of the raindrops (Eq. 1).

The method is a high-order optimal one which is exact when the function to integrate is smoother than a polynomial of degree (2n -1). The χ i symbols denote the weights and x i , the abscissas, then :

∞ 0 x β exp(-x) f (x) dx n i=1 χ i f (x i ).
Finally, the BCS coefficient (Eq. 4) is integrated over the AP diameters d p to obtain the BCS rate with respect to the total number of AP using Eq. 3:

∂N p ∂t BCS = ∞ 0 ∂n p (d p ) ∂t BCS dd p = - ∞ 0 γ(d p ) n p (d p ) dd p . ( 5
)
This time, the integration is performed using a Gauss-Hermite quadrature formula with appropriate weights χ i and abscissas x i :

∞ 0 exp(-x 2 ) f (x ) dx n i=1 χ i f (x i ).
As confirmed by a recent study [START_REF] Andronache | Estimated variability of below-cloud aerosol removal by rainfall for observed aerosol size distributions[END_REF], Fig. 1 shows the BCS coefficient γ(d p ) for increasing rainfall rates R from 0.1 to 100 mm.hr -1 .

The in Fig. 1 that the logarithm of the BCS coefficient is proportional to the logarithm of R, with a good approximation. These results are in agreement with those of previous studies [START_REF] Andronache | Estimated variability of below-cloud aerosol removal by rainfall for observed aerosol size distributions[END_REF][START_REF] Loosmore | Precipitation scavenging of atmospheric aerosols for emergency response applications: testing an updated model with new real-time data[END_REF]Tost et al., 2006;[START_REF] Feng | A 3-mode parametrization of below-cloud scavenging of aerosols for use in atmospheric dispersion models[END_REF].

rainfall rate R = ∞ 0 (π/6)(ρ w /ρ a )D 3 d U t (D d )n R (D d ) d(D d ) is

Mass BCS rate

The computation of the AP scavenging rate allows us to study the evolution of the number concentration for each AP mode (N f ree = N p ) and also to estimate the total mass of AP (m scav ) captured by the drops ; m scav constitutes the depleted part of the AP in mass.

The budget equation of m scav is composed of three terms: a transport tendency (marked "transp"), a mass BCS removal term ("BCS ") and an additional sedimentation term ("sed ") because m scav moves at the terminal fall speed of the drops.

∂m scav ∂t = ∂m scav ∂t transp + ∂m scav ∂t BCS + ∂m scav ∂t sed
The mass of the free AP which are scavenged by rain, is also computed with the Gauss-Hermite quadrature formula, so according to Eq. 3, one gets:

∂m scav ∂t BCS = - ∂m f ree ∂t BCS = π 6 ρ p ∞ 0 d p 3 γ(d p ) n p (d p ) dd p
where both AP sphericity and constant density for each particle mode are assumed.

The mass sedimentation rate of the scavenged particles is taken proportional to the sedimentation rate of the raindrop mixing ratio q R :

1 m scav ∂m scav ∂t sed = 1 q R ∂q R ∂t sed ,
suggesting that statistically, the mass of scavenged AP contained in an individual drop of size D r is proportional to the mass of the drop. One of the goals of COPS was to collect time series of high resolution surface data including as many aerosol, cloud, and precipitation variables as possible, to be used as lower boundary conditions for mesoscale models.

Looking at the BCS process, we are interested by two types of COPS data that are both in situ measurements at ∼ 3 m above ground level :

• Mean rainfall rate calculated each 10 min from precipitation measurements with an Optical Rain Gauge (ORG);

• Time evolution of AP concentrations sampled by a Grimm optical particles counter with a time step of 1 min. The AP concentration is provided for 15 classes of particles: size cut radius (µm) are [0.15; 0.2; 0.25; 0.325; 0.4; 0.5; 0.8; 1; 1.5; 2; 2.5; 3.75; 5; 7.5; 10; 20].

After comparing the two types of COPS data (precipitation intensity and AP concentrations), some sequences have been selected for which the aerosol depletion could be attributed to the wet removal by precipitation:

the assumption is made that other processes could be neglected because of the short time duration of these sequences. Fig. 2 shows the COPS AP concentration data for Julian day 184.38 to 184.8 as well as the precipitation rate attesting the importance of scavenging for large AP.

For the moment, only the first precipitation event (between 184.38 and 184.5) is examined. During this rainfall event, the AP concentration of classes 1 to 8 (r p < 1.5 µm), is not perturbed. In contrast, the AP concentrations of large diameters (classes 9 to 11) vary and are depleted according to the rainfall event. This depletion is partly explained by BCS because, as shown in Fig. 1, the BCS rate is the highest for large diameters (classes 9 to 11)

while small ones (1 to 8) correspond to the "Greenfield gap". Thus, the data suggest the potential effect of BCS upon the evolution of AP in the troposphere. Therefore, the implementation of this process in the aerosol continuity equation is necessary in cloud resolving models such as MesoNH.

In order to simulate the depletion of the AP by BCS using the first pre- For each class, the percentage of AP concentration depleted during the rainfall event, is calculated by comparing the AP concentration at the beginning and at the end of the event. This is done for measured and simulated data. As a result, the computed percentage does not reflect the entire depletion of each class of the COPS measurements. Nevertheless, the influence of the BCS effect is reproduced selectively as a function of the particle class diameter. For class number 3 (accumulation mode) and class number 10 (coarse mode) for instance (see the legend above Fig. 2), the simulated depletion percentages are 0.1% and 87% whereas the measured ones are 10.8%

and 54%, respectively. The discrepancy between the measured and simulated trends can be attributed to other processes such as turbulent transport, microphysics and to the assumption of homogeneous vertical profiles in the precipitating zone.

2D SIMULATIONS WITH MesoNH

The numerical experiments are performed to test and to evaluate the impact of BCS on a multimodal AP population. In all cases, particles are transported by a moist flow defined by warm or mixed-phase microphysical cycles producing rainfalls of varying characteristics along the course of the simulation. The rain mixing ratio

(q R = ∞ 0 (π/6)(ρ w /ρ a )D 3 d n R (D d ) d(D d ))
is used to determine the raindrop size distribution n R (D d ) and to integrate the BCS term (Eqs. 4 and 5) in MesoNH.

Warm shallow convection: the "HaRP" test case

The idealized "HaRP" test case comes from "The Hawaiian Rainband

Project" that took place in 1990 [START_REF] Szumowski | Simple two-dimensional kinematic framework designed to test warm rain microphysical models[END_REF]. A first run is illustrated in Figs 3a-d showing the particle concentrations (grey scale), the rain mixing ratios (contours of 0.1, 0.2, 0.5 and 1.0 g/kg)

It
and the flow circulation with arrows. The simulation is initialized for a single log-normal distribution of AP with a median diameter of 2 µm. The geometric standard deviation σ = 2 and the initial concentration N = 3 cm -3 are representative of marine conditions (see mode "b" in Table 1 of [START_REF] Andronache | Estimated variability of below-cloud aerosol removal by rainfall for observed aerosol size distributions[END_REF]). At the beginning of the simulation, the moist flow converges horizontally in the low levels initiating the formation of a cloud cell. Quickly after reaching the autoconversion threshold, the cloud droplets convert partially into raindrops which start to precipitate and to evaporate below cloud base. At the same time, the AP are entrained in the ascending motion of the flow (see the AP concentration in Fig. 3a taken after 20 min).

After 25 min (Fig. 3b), the AP continue to rise through the rainshaft where they are partially depleted by BCS. Then the rain reaches its maximum intensity at 35 min therefore leading to an efficient depletion of AP by BCS in Fig. 3c. At the end of the simulation (so after 50 min in Fig. 3d), the AP flux diverges at the top of the cloud showing a nearly 50 % reduction of the AP concentration by BCS.

In summary, the 2D simulation reproduces an interesting situation where AP are transported and accumulated by a forced convective flow in a region where raindrops are formed, precipitate and then lead to BCS. A key-aspect of our BCS treatment is that at high resolution, the BCS efficiency in a 2D/3D precipitating system is thoroughly linked to the flow pattern and to the microphysical scheme. This aspect was underevaluated in most previously cited studies where observed or simulated rainfall rates are used to estimate the BCS rate of AP.

A generalization of this "HaRP" case study has been undertaken in order to explore the sensitivity of several AP size distributions to the BCS. In this series of experiment, the median diameter d p varies logarithmically in the interval 10 -3 -10 1 µm while the geometric standard deviation σ is taken in the A similar series of simulation performed with a Gamma (α = 1, ν = 2) raindrop size distribution (Eq. 1) leads to a further slight enhancement of BCS (not shown here). This effect is explained simply by considering the higher concentration of big drops in the "Gamma" case compared to the Marshall-Palmer case.

The "HaRP" results show the difficulties to predict the potential scavenging of multiple AP modes without an accurate BCS module coupled to a microphysical scheme and taking into account the transport of the particles.

The numerous combinations of AP log-normal distribution parameters and rain intensity clearly leads to a large variety of scavenging efficiencies.

Tropical squall line: the "COPT" test case

The "COPT" test case [START_REF] Caniaux | A numerical study of the stratiform region of a fast-moving squall line. Part I: General description and water and heat budgets[END_REF] Although liquid drops and frozen particles can coexist in a squall line, the study focuses on the BCS of AP by rain. In mixed-phase clouds, impaction scavenging by ice is enhanced by phoretic effects [START_REF] Martin | A theoretical determination of the efficiency with which aerosol particles are collected by simple ice crystal plates[END_REF] and probably, by electric forces [START_REF] Miller | A model for the determination of the scavenging rates of submicron aerosols by snow crystals[END_REF][START_REF] Andronache | Diffusion and electric charge contributions to belowcloud wet removal of atmospheric ultra-fine aerosol particles[END_REF] as these clouds become electrified by non-inductive charging mechanisms involving ice-ice shocks [START_REF] Mansell | Electrification and lightning in simulated supercell and non-supercell thunderstorms[END_REF][START_REF] Barthe | Description and first results of an explicit electrical scheme in a 3D cloud resolving model[END_REF]. This last effect clearly makes the in-cloud impaction scavenging of interstitial AP by ice difficult to describe in detail. However as the ice nucleation scavenging should dominate according to [START_REF] Alheit | A theoretical study of the wet removal of atmospheric pollutants. part iv: the uptake and redistribution of aerosol particles through nucleation and impaction scavenging by growing cloud drop and ice particles[END_REF], the impaction scavenging of AP by ice particles is probably of secondary importance. As a result and because we concentrate on (tropical) warm precipitation events, in-and below-cloud scavenging effects of AP by ice and snow are not considered in this study.

The aim of this case experiment is to explore the BCS efficiency of AP layers in a long-lived mixed-phase precipitating system. In this case, we show that the residency time of air parcels (see the discussion about the merits of this aging tracer in Lucas-Picher et al. ( 2008)) is a useful tool to interpret the temporal aspects of the AP removal in rainy areas. Therefore, a series of simulations, without (SCAV 0 ) and with (SCAV 1 ) application of BCS to the dust particles, are performed. In each case, the dust particles are transported by the two-dimensional circulation that develops in the organized convective system. Figs 5a-c show the difference of dust concentrations between the twin simulations SCAV 0 and SCAV 1 taken after 9 hours of simulation. The plots of AP field differences allows to picture the net BCS component of the dust particle budget in a complex precipitating system.

The squall line develops in a series of narrow convective cells on the left followed by a large stratiform area with much less precipitation on the right of the figures. Here, the plotted field difference (SCAV 0 -SCAV 1 ) focuses on the coarse mode concentration of the dust particle distribution. The plotted fraction of the dust particles is always positive because of the reduced fraction of particle concentration when the BCS process is considered. The upper glaciated part of the cloud system and the raining cells are also depicted in the figures. The initial concentration of the coarse mode is 20 cm -3 in all the simulations. A rapid inspection of Figs 5a-c indicates that all the dust layers are affected by the scavenging process which is efficient enough to eliminate much more than 1 cm -3 (scale maxima) of coarse mode particles in the heavy precipitating convective region of the squall line.

As expected the lowest dust layer (Fig. 5a) is the most affected by the Fig. 5c shows that the BCS process still perturbs the field of AP when the dust particles are concentrated in the [4-6] km layer. Surprisingly the low-level rainy area of the squall line where convective cells develop, reveals a lot of scavenged particles even if they were not present initially in this region.

The explanation is that the dust particles are recycled by the density current of the squall line [START_REF] Lafore | A numerical investigation of the organization and interaction of the convective and stratiform regions of tropical squall lines[END_REF]. Clearly speaking, it is the internal dynamics of the squall line (convective cells, mesoscale subsidence and rear-to-front flow below the stratiform region) that transports the air and the dust of the midlevels in a downward branch, and which is responsible for the BCS depletion in the rainy area ahead. This effect is illustrated in Fig. 6 where the time evolution of a selected vertical profile (taken at As expected, BCS is less affecting the upper glaciated region of Fig. 5c in comparison to that of Fig. 5a. The reason is that rain comes from the melting of graupel particles at 4.5 km altitude leading to a very shallow layer where BCS can operate locally during the ascent of the dust particles.

In order to illustrate the importance of transport, here combined to BCS in the squall line (see also [START_REF] Lafore | A numerical investigation of the organization and interaction of the convective and stratiform regions of tropical squall lines[END_REF]), the residency time of air parcels having been in contact with rain, is shown in Fig. 5d after 9 hours of simulation with half-hourly contours starting after 2 hours of simulation.

The residency time (RT) is a 3D tracer which is easy to integrate along the course of the simulation. It includes a transport term followed by a conditional time incrementation step, here characterized by a "rain" mask.

At each time step, the rainy areas are selected by considering the gridpoints passing a mask where the rain mixing ratio is larger than 0.01 g/kg and so where RT is augmented by the value of the timestep.

Without "rain" mask (not shown here), the bulk structure of RT reproduces the mean wind shear with inflow (outflow) conditions below (above) 7.5 km height at the left boundary. More interestingly with the "rain" mask, there is a close resemblance between the fine scale structures of RT in Fig. 5d and the scavenged AP fields in Fig. 5a-c. For instance the wedge pattern of the scavenged AP in the low levels is reflected by a continuous increase of RT at the same place. Similarly, holes of scavenged AP at the melting level (4.5 km) correspond approximately to a reduced RT, meaning that fresher air parcels are less depleted by BCS. Finally, Fig. 5d reveals that a great amount of air in the glaciated part of the squall line has been exposed to the rain during 30 min at least. This demonstrates that the amount of scavenged AP depends on the rainfall rate but that it is also strongly correlated to the residency time of AP under the rainshaft.

The last point to examine is the mass of scavenged AP by wet deposition within the underlying surface of the COPT squall-line. Fig. 7a presents the surface cumulated mass of scavenged AP, assuming a density ρ p = 2600 kg.m -3 , for the 3 layers of AP after 9 hours of simulation. The cumulated rainfall, peaking slightly downwind at 82 mm.hr -1 , is plotted in Fig. 7b for comparison. Fig. 7a shows that BCS is most pronounced for the mid layer ([2-4] km) case where it reaches 0.32 kg.m -2 . Surprisingly, more AP are scavenged from the upper layer ([4-6] km) than from the lower layer ([0-2] km). The explanation comes from the observation that the big drops which are the most efficient to scavenge AP, are present in the top layer while they are partially evaporated when they reach the bottom layer.

In summary, the COPT case offers many opportunities to investigate different aspects of AP scavenging by rain in an organized convective system where resolved flow dynamics and transport of tracers are well analysed [START_REF] Lafore | A numerical investigation of the organization and interaction of the convective and stratiform regions of tropical squall lines[END_REF]. Here the implementation of a performant PPM transport scheme in the model is also crucial to get a deep insight into the complexity of the net BCS aspect of AP at mesoscale.

CONCLUSIONS

A Here the location, with respect to the freezing level, of a layer of dust particles is a key-parameter to analyse the net effect of BCS. Some details of the depleted AP fields can be interpreted with the help of the time residency of air parcels in the rainshaft.

In conclusion, BCS of particles is expected to remove giant condensation nuclei (sea salt particles) and large ice forming nuclei (mineral dust particles)

by precipitation with possible consequences on the nucleating properties of drizzling marine stratus and convective clouds with an ice phase. Therefore the next step of this work will address questions about the nucleation scavenging processes. The combination between BCS and the activation of cloud condensation nuclei is still an intriguing process in the short range evolution of the cloud cover of clean precipitating marine boundary layers (e.g. the formation of open cells in marine stratocumulus [START_REF] Rosenfeld | Switching cloud cover and dynamical regimes from open to closed Benard cells in response to the suppression of precipitation by aerosols[END_REF]

  expresses the number of AP collected by the falling raindrop in the volume of air swept out by the raindrop. A value of E = 1 implies that the raindrop collects all the particles encountered along the path of the fall. Most of the recent 1 The residency time of falling raindrops (H/V f all ∼ 2 -5 × 10 2 s where H ∼ 1000 m is the height of the cloud base and V f all ∼ 2 -5 m/s, the fallspeed of the raindrops) is short compared to the time scale of the BCS (Λ -1 > 10 4 s according to Sportisse (2007)) so, a possible modification of the median diameter d pi of the AP size distribution is clearly a second-order effect compared to the variability of the rainfall rates.

  a function of the slope parameter λ R which defines the Marshall-Palmer distribution of the raindrops (see Eq.1). Three distinct regimes of the BCS coefficient, corresponding to those of the collision efficiency E, i.e. Brownian motion, interception and inertial impaction, can be seen in Fig. 1 with a marked "Greenfield gap" for 0.1 µm< d p < 1 µm. The coarse AP are the most efficiently scavenged, the ultrafine particles are moderately scavenged and the scavenging is low for the intermediate mode. Moreover one can notice

  3. QUALITATIVE APPLICATION TO THE COPS EXPERIMENTCOPS (Convective and Orographically-induced Precipitation Study) was a 3 month international field campaign with part of a program aiming at improving precipitation forecasts(Wulfmeyer and coll., 2008). COPS took place in Summer 2007 in south-western Germany and north-eastern France.

  cipitation event from the COPS data, a mean particle diameter d m has been defined for each class and the corresponding BCS coefficient γ(d m ) calculated. The method consists in initializing the AP concentrations with mea-sured values. Then the mean rainfall intensity data are incorporated in the BCS module for the selected rainy episode. The temporal evolution of the AP concentrations of each class of size N i is only driven by the calculated BCS rates.

  scheme for the microphysics and a highly accurate PPM scheme for the transport of the scalar fields. The computational domain extends over 180×60 gridpoints with a spacing of 50 m in the x and z directions. The time step is 5 sec. The concentrations of the particle modes are continuously replenished through the open lateral boundaries in the lowest 250 m where inflow conditions occur.

  [1-6] range (see the compiled size distribution cases of[START_REF] Andronache | Estimated variability of below-cloud aerosol removal by rainfall for observed aerosol size distributions[END_REF]).The results obtained for a Marshall-Palmer raindrop size distribution are plotted in Fig.4. The figure represents the percentage of removed AP after 50 min of simulation. These scavenged AP concentrations are time and space averaged in the regions where the rain mixing ratio is larger than 0.01 g/kg.Inspection of Fig.4, shows a continuous enhancement of the BCS with σ for d p < 1 µm. This feature results from the fact that increasing the width of the AP log-normal distribution and holding d p fixed, increases the AP size interval over which the scavenging efficiency is integrated (see Fig.1).Consequently and because the weighting effect of the scavenging coefficient of the smallest particles dominates, the BCS increases with σ. The situation appears more complicated in the "Greenfield gap" interval 0.1 -2 µm of d p where the BCS efficiency of AP is reduced. In this region however, the large sensitivity of the AP size distributions is fully explained by the head and the tail of the AP size distribution which involve higher and higher scavenging efficiencies as σ increases. Moving to d p > 2. µm, i.e. to coarse AP mode, one can notice that the BCS is no more dependent on the value of σ. d p is the leading parameter for the BCS in this AP size range. As a result, the three distinct regimes which show up in Fig.4are clearly correlated to the variations of the scavenging coefficients shown in Fig.1.

  is typically a 12 hour simulation of a tropical squall line with kilometer scale-resolved internal circulations, a 2D turbulence scheme and mixed-phase microphysics. The model is initialized with 3 successive layers of aerosol of 2 km depth starting from the ground level. For each layer, the same multimodal population of particles corresponding to the "Dust Layer" case reported by[START_REF] Andronache | Estimated variability of below-cloud aerosol removal by rainfall for observed aerosol size distributions[END_REF], is chosen for the simulation of an African squall line. The domain contains 320×44 gridpoints unevenly spaced in the vertical (∆z = 70 m at ground level and ∆z = 700 m above 12 km). The horizontal resolution is 1.25 km.The model is integrated with a time step of 10 sec. A gravity wave damping layer is inserted between 17 km and the model top at 22.5 km. A constant speed transformation is used to compensate for the motion of the squall line.No fluxes are considered in the surface layer. Convection is initiated by forming a -0.01 K.s -1 artificial cold pool in the low levels of a small domain during 10 min.

  Figs 5a-c differ only by their respective initialization: the initial dust layer was between [0-2] km for Fig. 5a, [2-4] km for Fig. 5b and [4-6] km for Fig. 5c. To summarize, the figures are showing at the same scale, the net scavenging efficiency of each initial layer of coarse dust particle mode (d p1 = 0.55 µm, σ = 2.5).

  scavenging removal because of stronger convergence and circulation flow intensity in the lower levels, and a longer time exposure to the precipitating area. Another interesting feature displayed in Figs5a-cis the fate of the net scavenged particles in the upper glaciated region of the squall line, that is well above the ice freezing level at 4.5 km height. In Fig.5athe relatively large difference of concentration well above the ice melting layer reflects the indirect scavenging impact in the upper glaciated region : the non-scavenged particles in simulation SCAV 0 are transported upwards while they are depleted in simulation SCAV 1 due to efficient BCS in the low levels. Therefore the particle concentration in the squall line upper region is clearly sensitive to BCS of low levels AP sources with potential feedback on the ice nucleation process by selectively modifying the ice nuclei spectra (for instance contact freezing nuclei which efficiency depends on particle size). The same analysis holds for Fig.5bwith a lower scavenging efficiency and less upward transport of the "mid-level" particle layer in the convectively rainy region (leftmost precipitating cells) because of a shorter residency time. In contrast, scavenging is enhanced in the trailing stratiform part of the squall line (right part of the [0-2] km layer in Figs5a-b). This feature is explained by the presence of bigger raindrops (raindrops partially evaporated during their fall) which are more efficient to scavenge AP in the [2-4] km layer.

  x = 225 km) of the [4-6] km AP layer, is plotted with the rain mixing ratio. Three stages are identified in the figure. During the first 3 hours of simulation, the initial vertical distribution of AP is relatively undisturbed. Then, AP are appearing in the lower(upper) levels due to down(up)ward transport by the first convective cells of the squall line. Finally after 6 hours, the AP are continuously filling the atmosphere up to 12 km height. During this period of time, rainy cells are precipitating in the lowest 4.5 km allowing raindrops to scavenge the transported AP.

  BCS module has been implemented in the mesoscale/cloud resolving model MesoNH. The module is based on the collection efficiencies of Slinn (1983) with an accurate and optimized numerical integration over unbounded raindrop and AP size distributions. Then the study explores the sensitivity of BCS to various aspects of the initial AP size distribution and for two precipitating events. Here in contrast with 0D/1D studies and experiments performed at larger scale, the combined effect of transport and BCS acting on AP is stressed and illustrated in the framework of 2D simulations. The HaRP numerical experiment underlines the size dependent depletion of an AP population by BCS as illustrated by changing the median diameter of log-normal distributions: the ultrafine particles are depleted, the intermediate ones remain unaltered and the large particles are scavenged dramatically. The results show also a significant sensitivity to the standard geometric deviation parameter of the AP size distributions. The simulation of the COPT squall line demonstrates the importance of coupling transport and microphysics to track the BCS effects in organized precipitating system.

Figure 1 :Figure 3 :

 13 Figure 1: Scavenging coefficient as a function of particle size and rain rate (calculated with Eq.(4))

Figure 4 :Figure 5 :

 45 Figure 4: Fraction of scavenged AP as a function of the median diameter d p and of the geometric standard deviation σ of the AP log-normal size distribution. The results obtained for a Marshall-Palmer distribution of the raindrops after 50 min of simulation of the "HaRP" case.

Figure 6 :Figure 7 :

 67 Figure 6: Time evolution of vertical profiles taken at x=225 km. The coarse mode particles [cm -3 ] in the [4-6] km layer case are displayed with the rain mixing ratio with solid contours superimposed in log scale from 10 -5 to 10 -3 kg/kg.

Table 1 :

 1 ). Definition and unit of symbols used in this article.

	    	b = 0.8 a = 842	Parameters of U t (D d )	N pi	Number concentration of AP mode i
	d p		AP diameter	n pi (d p )	AP size distribution of mode i
	d pi	Median diameter of mode i	n R (D d ) Raindrop size distribution
	D d	Drop diameter	q R	Rain mixing ratio
	D		AP diffusivity	R	Rainfall rate
	E(D d , d p ) AP-droplets collision efficiency	Re	Raindrop Reynolds number
	f (x i )	Function to integrate	r p	AP radius
	i, l	Index of AP mode	Sc	AP Schmidt number
	m f ree	Mass of free AP	St	AP Stokes number
	m scav	Total mass of scavenged AP	St *	Critical Stokes number
	n		Order of the Gauss quadrature	U t (D d )	Terminal fall velocity of raindrops
	N 0	Marshall-Palmer parameter	x i , x i	Abcissa for Gauss quadrature
	N f ree	Number concentration of free AP
	α		Parameters of n R (D d )	Ψ(d p , t) Mass or number AP size distribution
	β		Parameter of the Gamma function	ρ a0	Reference air density at ground level
	γ(d p )	BCS coefficient	ρ a	Reference density of dry air
	Γ		Gamma function	ρ p	AP density
	λ R	Marshall-Palmer slope parameter	σ i	Standard deviation of mode i
	µ a	Dynamic viscosity of dry air	τ a	Relaxation time of the collected AP
	µ w	Water viscosity	χ i , χ i	Weights for Gauss quadrature
	ν		Parameters of n R (D d )	ω	Viscosity ratio
	φ		Ratio of diameters

Acknowledgements

The authors are deeply grateful to T. Maric, now at University of Washington, Seattle, WA, for the implementation of the PPM transport scheme of

 [START_REF] Colella | The piecewise parabolic method (ppm) for gas-dynamical simulations[END_REF]in MesoNH. We also thank C. Galy-Lacaux of Laboratoire d'Aérologie for useful discussion on the wet deposition of aerosol particles and J. Duron for her precious help in preparing some of the figures.


Wenjuan Jiang
email: wenjuan.jiang@gipsa-lab.grenoble-inp.fr

Alexandre Kruszewski
email: alexandre.kruszewski@ec-lille.fr

Emilia Fridman
email: emilia@eng.tau.ac.il

Jean-Pierre Richard
email: jean-pierre.richard@ec-lille.fr

Delay dependent stability analysis of interval time-delay systems

Keywords: time-varying delay, switching systems, exponential stability, Lyapunov-Krasovskii functional, LMI

This paper considers interval time-varying delay systems. The time-delay interval is divided into several zones and the systems switch among the different zones. Based on Lyapunov-Krasovskii functional methods and linear matrix inequality (LMI) techniques, Exponential stability is exploited for every time-delay zone. The global stability of the switched system is guaranteed if some minimum average dwell time conditions are satisfied. Some numerical examples and comparisons with other works show that the methods greatly enlarge the value of maximum upper-bound of time-delay for the systems.

INTRODUCTION

A variety of stability and control techniques have been developed for general time delay systems [START_REF] Niculescu | Delay effects on stability: a robust control approach[END_REF][START_REF] Richard | Time delay systems: an overview of some recent advances and open problems[END_REF][START_REF] Chiasson | Applications of time delay systems[END_REF]. Of course, when a delay is introduced into a system, the achievable performance (speed, robustness) highly depends on the delay. In the case of a time-varying delay, the guaranteed performance generally depends on the interval of variation of the delay. To give an example, controlling a system through a communication network induces variable time delays and the resulting performance depends on the quality of service (QoS) the network can offer [START_REF] Jiang | A gain scheduling strategy for the control and estimation of a remote robot via internet[END_REF][START_REF] Seuret | Networked control using gps synchronization[END_REF].

In order to reduce the conservatism, some authors used the piecewise analysis method which consist in study the system for smaller delay intervals with a common Lyapunov-Krasovskii Functional (LKF) [START_REF] Yue | A piecewise analysis method to stability analysis of linear continuous/discrete systems with time-varying delay[END_REF]. Another approach is consider the system as a switching system in which each subsystem has a different delay. These delays belong to smaller intervals.

Switched systems are dynamical hybrid systems consisting of a family of continuous-time subsystems and a logical rule that orchestrates the switching between them (Liberzon D., 2003; [START_REF] Hirche | Performance oriented control over networks -switching controllers and switched time delay[END_REF]. In general, the stability of the subsystems themselves are not sufficient for the stability of the overall system. A lot of studies have dealt with the stability analysis and design of switched systems [START_REF] Liberzon | Stability of switched systems: a lie-algebraic condition[END_REF][START_REF] Hespanha | stability of switched systems with average dwell-time[END_REF][START_REF] Xie | Stabilization of switched linear systems with time-delay in detection of switching signal[END_REF][START_REF] Hirche | Towards quality-of-service control of networked control systems: a switched time delay systems approach[END_REF].

Single [START_REF] Sun | Stability of switched systems with time-varying delays: delay-dependent common lyapunov functional approach[END_REF][START_REF] Jiang | A gain scheduling strategy for the control and estimation of a remote robot via internet[END_REF] and multiple LKF analysis methods are the most frequently applied in the stabilization of the whole system. A common LKF can ensure the stability with arbitrary switching but it does not always exist. As an alternative, dwell-time based switching is considered in (Hespanha J.P. and Morse A.S., 1999) but time delay is not considered, while in [START_REF] Yan | Stability analysis of switched time delay systems[END_REF][START_REF] Hirche | Towards quality-of-service control of networked control systems: a switched time delay systems approach[END_REF], the time delay is a constant one.

In this paper, a linear switched time-delay system is considered. Multiple LKF are used in the stability analysis. The time-delay interval is divided into several ones to reduce conservativeness. In this means, the maximum timedelay can be enlarged. The minimum average dwell time is considered to guarantee the global stability of the switched system. The LMI (Linear Matrix Inequalities) conditions are derived following the idea of [START_REF] Park | Stability and robust stability for systems with a time-varying delay[END_REF] (using the convex approach). The first part introduces some sufficient conditions for exponential stability of linear systems with interval time-varying delay. Our main results are presented in the second part , where the linear time-delay system is analyzed as a switched system and where sufficient condition for the global stability are given by using average dwelling time. The last part gives some examples and comparisons with other works.

Notation: Throughout the paper the superscript '𝑇 ' stands for matrix transposition, ℛ 𝑛 denotes the 𝑛 dimensional Euclidean space with vector norm ∥ ⋅ ∥, ℛ 𝑛×𝑚 is the set of all 𝑛 × 𝑚 real matrices, and the notation 𝑃 > 0, for 𝑃 ∈ ℛ 𝑛×𝑛 means that 𝑃 is symmetric and positive definite. The symmetric elements of the symmetric matrix will be denoted by * . 𝐿 2 is the space of square integrable functions 𝑣 : [0, ∞) → 𝑅 𝑛 with the norm

∥𝑣∥ 𝐿2 = [ ∫ ∞ 0 ∥𝑣(𝑡)∥ 2 𝑑𝑡] 1/2 .

EXPONENTIAL STABILITY OF INTERVAL TIME-DELAY SYSTEMS

For the sake of simplicity, the LMI conditions are first derived for exponential stability of a linear time-delay system without switches.

Consider a linear time-delay system:

ẋ(𝑡) = 𝐴𝑥(𝑡) + 𝐴 1 𝑥(𝑡 -𝜏 (𝑡)) (1)
with 𝐴 stable and 𝜏 (𝑡), a fast varying interval delay with the bounds as [ℎ 1 , ℎ 2 ].

In the following, the notation 𝑥 𝑡 stands for 𝑥(𝑡). Applying the following LKFs [START_REF] Fridman | Exponential stability of linear distributed parameter systems with timevarying delays[END_REF]:

𝑉 (𝑡, 𝑥 𝑡 , ẋ𝑡 ) = 𝑥 𝑇 (𝑡)𝑃 𝑥(𝑡) + 1 ∑ 𝑖=0 𝑡-ℎ𝑖 ∫ 𝑡-ℎ𝑖+1 𝑒 2𝛼(𝑠-𝑡) 𝑥 𝑇 (𝑠)𝑆 𝑖 𝑥(𝑠)𝑑𝑠 + 1 ∑ 𝑖=0 (ℎ 𝑖+1 -ℎ 𝑖 ) -ℎ𝑖 ∫ -ℎ𝑖+1 𝑡 ∫ 𝑡+𝜃 𝑒 2𝛼(𝑠-𝑡) ẋ𝑇 (𝑠)𝑅 𝑖 ẋ(𝑠)𝑑𝑠𝑑𝜃 (2)
where ℎ 0 = 0, 𝑃 > 0 and 𝑅 𝑖 , 𝑆 𝑖 ≥ 0. Saying that the system (1) is exponentially stable means that there exist a positive function 𝑊 (𝑡, 𝑥(𝑡), ẋ(𝑡)) such that:

𝑊 (𝑡, 𝑥 𝑡 , ẋ𝑡 ) ≤ 𝑒 -2𝛼(𝑡-𝑡0) 𝑊 (𝑡, 𝑥 𝑡0 , ẋ𝑡0 ) ( 3 )
It can be ensured by using the LKFs (2) and by checking if: V (𝑡, 𝑥 𝑡 , ẋ𝑡 ) + 2𝛼𝑉 (𝑡, 𝑥 𝑡 , ẋ𝑡 ) < 0 ( 4 )

After the differentiation of the function (2), the system is exponentially stable if:

V (𝑡, 𝑥 𝑡 , ẋ𝑡 ) + 2𝛼𝑉 (𝑡, 𝑥 𝑡 , ẋ𝑡 ) ≤ 2𝑥 𝑇 (𝑡)𝑃 ẋ(𝑡) +2𝛼𝑥 𝑇 (𝑡)𝑃 𝑥(𝑡) + ẋ𝑇 (𝑡)[ 1 ∑ 𝑖=0 (ℎ 𝑖+1 -ℎ 𝑖 ) 2 𝑅 𝑖 ] ẋ(𝑡) - 1 ∑ 𝑖=0 (ℎ 𝑖+1 -ℎ 𝑖 )𝑒 -2𝛼ℎ𝑖+1 𝑡-ℎ𝑖 ∫ 𝑡-ℎ𝑖+1 ẋ𝑇 (𝑠)𝑅 𝑖 ẋ(𝑠)𝑑𝑠 + 1 ∑ 𝑖=0 𝑥(𝑡 -ℎ 𝑖 ) 𝑇 𝑒 -2𝛼ℎ𝑖 𝑆 𝑖 𝑥(𝑡 -ℎ 𝑖 ) - 1 ∑ 𝑖=0 𝑒 -2𝛼ℎ𝑖+1 𝑥 𝑇 (𝑡 -ℎ 𝑖+1 )𝑆 𝑖 𝑥(𝑡 -ℎ 𝑖+1 ).
(5)

Noticing that:

𝑡-ℎ0 ∫ 𝑡-ℎ1 ẋ𝑇 (𝑠)𝑅 𝑖 ẋ(𝑠)𝑑𝑠 = 𝑡-𝜏 (𝑡) ∫ 𝑡-ℎ1 ẋ𝑇 (𝑠)𝑅 0 ẋ(𝑠)𝑑𝑠 + 𝑡-ℎ0 ∫ 𝑡-𝜏 (𝑡) ẋ𝑇 (𝑠)𝑅 0 ẋ(𝑠)𝑑𝑠
and applying the Jensen's inequality [START_REF] Gu | Stability of time-delay systems[END_REF])

𝑡-ℎ0 ∫ 𝑡-ℎ1 ẋ𝑇 (𝑠)[(ℎ 2 -ℎ 1 )𝑅 0 ] ẋ(𝑠)𝑑𝑠 ≥ 𝑡-ℎ0 ∫ 𝑡-ℎ1 ẋ𝑇 (𝑠)𝑑𝑠𝑅 0 𝑡-ℎ0 ∫ 𝑡-ℎ1 ẋ(𝑠)𝑑𝑠, 𝑡-ℎ1 ∫ 𝑡-𝜏 (𝑡) ẋ𝑇 (𝑠)(ℎ 2 -ℎ 1 )𝑅 𝑗 ẋ(𝑠)𝑑𝑠 ≥ ℎ 2 -ℎ 1 𝜏 -ℎ 1 𝑡-ℎ1 ∫ 𝑡-𝜏 (𝑡) ẋ𝑇 (𝑠)𝑑𝑠𝑅 1 𝑡-ℎ1 ∫ 𝑡-𝜏 (𝑡) ẋ(𝑠)𝑑𝑠, 𝑡-𝜏 (𝑡) ∫ 𝑡-ℎ2 ẋ𝑇 (𝑠)[(ℎ 2 -ℎ 1 )𝑅 1 ] ẋ(𝑠)𝑑𝑠 ≥ ℎ 2 -ℎ 1 ℎ 2 -𝜏 𝑡-𝜏 (𝑡) ∫ 𝑡-ℎ2 ẋ𝑇 (𝑠)𝑑𝑠𝑅 1 𝑡-𝜏 (𝑡) ∫ 𝑡-ℎ2 ẋ(𝑠)𝑑𝑠. ( 6
)
Here for 𝜏 → ℎ 1 we understand by

1 𝜏 (𝑡) -ℎ 1 𝑡-ℎ1 ∫ 𝑡-𝜏 (𝑡) ẋ(𝑠)𝑑𝑠 = lim 𝜏 →ℎ1 1 𝜏 (𝑡) -ℎ 1 𝑡-ℎ1 ∫ 𝑡-𝜏 (𝑡) ẋ(𝑠)𝑑𝑠 = ẋ(𝑡 -ℎ 1 ). For ℎ 2 -𝜏 (𝑡) → 0 the vector 1 ℎ2-𝜏 (𝑡) ∫ 𝑡-𝜏 (𝑡) 𝑡-ℎ2 ẋ(𝑠)𝑑𝑠 is defined similarly as ẋ(𝑡 -ℎ 2 ).
Denoting

𝑣 1 = 1 𝜏 -ℎ 1 𝑡-ℎ1 ∫ 𝑡-𝜏 (𝑡) ẋ(𝑠)𝑑𝑠, 𝑣 2 = 1 ℎ 2 -𝜏 𝑡-𝜏 (𝑡) ∫ 𝑡-ℎ2 ẋ(𝑠)𝑑𝑠, (7) we obtain V (𝑡, 𝑥 𝑡 , ẋ𝑡 ) + 2𝛼𝑉 (𝑡, 𝑥 𝑡 , ẋ𝑡 ) ≤ 2𝑥 𝑇 (𝑡)𝑃 ẋ(𝑡) +2𝛼𝑥 𝑇 (𝑡)𝑃 𝑥(𝑡) + ẋ𝑇 (𝑡) 1 ∑ 𝑖=0 (ℎ 𝑖+1 -ℎ 𝑖 ) 2 𝑅 𝑖 ) ẋ(𝑡) + 1 ∑ 𝑖=0 𝑒 -2𝛼ℎ𝑖 𝑥 𝑇 (𝑡 -ℎ 𝑖 )𝑆 𝑖 𝑥(𝑡 -ℎ 𝑖 ) - 1 ∑ 𝑖=0 𝑒 -2𝛼ℎ𝑖+1 𝑥 𝑇 (𝑡 -ℎ 𝑖+1 )𝑆 𝑖 𝑥(𝑡 -ℎ 𝑖+1 ) -[𝑥(𝑡) -𝑥(𝑡 -ℎ 1 )] 𝑇 𝑒 -2𝛼ℎ1 𝑅 0 [𝑥(𝑡) -𝑥(𝑡 -ℎ 1 )] -(𝜏 -ℎ 1 )(ℎ 2 -ℎ 1 )𝑣 𝑇 1 𝑒 -2𝛼ℎ2 𝑅 1 𝑣 1 -(ℎ 2 -𝜏 )(ℎ 2 -ℎ 1 )𝑣 𝑇 2 𝑒 -2𝛼ℎ2 𝑅 1 𝑣 2 . ( 8
)
We note that in the latter bound we can substitute for 𝜏 its upper bound ℎ 2 .

We insert free-weighting 𝑛×𝑛-matrices (He Y. et al., 2004[START_REF] He | Delayrange-dependent stability for systems with time-varying delay[END_REF] by adding the following expressions to V :

0 = 2[𝑥 𝑇 (𝑡)𝑌 𝑇 1 + ẋ𝑇 (𝑡)𝑌 𝑇 2 + 𝑥 𝑇 (𝑡 -𝜏 )𝑇 𝑇 1 ][-𝑥(𝑡 -ℎ 1 ) +𝑥(𝑡 -𝜏 ) + (𝜏 -ℎ 1 )𝑣 1 ], 0 = 2[𝑥 𝑇 (𝑡)𝑍 𝑇 1 + ẋ𝑇 (𝑡)𝑍 𝑇 2 ][𝑥(𝑡 -ℎ 2 ) +(ℎ 2 -𝜏 )𝑣 2 -𝑥(𝑡 -𝜏 )]. (9)
We use further the descriptor method [START_REF] Fridman | New lyapunov-krasovskii functionals for stability of linear retarded and neutral type systems[END_REF], where the right-hand side of the expression

0 = 2[𝑥 𝑇 (𝑡)𝑃 𝑇 2 + ẋ𝑇 (𝑡)𝑃 𝑇 3 ][𝐴𝑥(𝑡) + 𝐴 1 𝑥(𝑡 -𝜏 (𝑡)) -ẋ(𝑡)], (10)
with some 𝑛 × 𝑛-matrices 𝑃 2 , 𝑃 3 is added into the righthand side of (8).

Setting 𝜂(𝑡) = 𝑐𝑜𝑙{𝑥(𝑡), ẋ(𝑡), 𝑥(𝑡-ℎ 1 ), 𝑥(𝑡-ℎ 2 ), 𝑣 1 , 𝑣 2 , 𝑥(𝑡-𝜏 )}, we obtain that along the time-delay system (1) V (𝑡, 𝑥 𝑡 , ẋ𝑡 ) + 2𝛼𝑉 (𝑡, 𝑥 𝑡 , ẋ𝑡 ) ≤ 𝜂 𝑇 (𝑡)Φ𝜂(𝑡) < 0, (11)

if the LMI ( 12) is feasible,

where 𝜏 = ℎ 2 and R𝑖 = 𝑒 -2𝛼ℎ𝑖+1 𝑅 𝑖 , S𝑖 = 𝑒 -2𝛼ℎ𝑖+1 𝑆 𝑖 , Ŝ𝑖 = 𝑒 -2𝛼ℎ𝑖 𝑆 𝑖 , Φ 11 = 𝐴 𝑇 𝑃 2 + 𝑃 𝑇 2 𝐴 + 𝑆 0 -𝑒 -2𝛼ℎ1 𝑅 0 + 2𝛼𝑃, Φ 12 = 𝑃 -𝑃 𝑇 2 + 𝐴 𝑇 𝑃 3 , Φ 22 = -𝑃 3 -𝑃 𝑇 3 + 1 ∑ 𝑖=0 (ℎ 𝑖+1 -ℎ 𝑖 ) 2 𝑅 𝑖 , Φ 16 = 𝑃 𝑇 2 𝐴 1 , Φ 26 = 𝑃 𝑇 3 𝐴 1 , 𝜙 3 = -( S0 + R0 -Ŝ1 ), 𝜙 4 = -S1 , 𝜙 5 = -(ℎ 2 -ℎ 1 )(𝜏 -ℎ 1 ) R1 , 𝜙 6 = -(ℎ 2 -ℎ 1 )(ℎ 2 -𝜏 ) R1 . ( 15
)
The latter condition leads for 𝜏 → ℎ 1 and for 𝜏 → ℎ 2 to the following LMIs ( 13) and ( 14), where we deleted the zero column and the zero row: Denote by 𝜂 𝑖 (𝑡) = 𝑐𝑜𝑙{𝑥(𝑡), ẋ(𝑡), 𝑥(𝑡 -ℎ 1 ), 𝑥(𝑡 -ℎ 2 ), 𝑣 𝑖 , 𝑥(𝑡 -𝜏 )}, 𝑖 = 1, 2.

Then ( 13), ( 14) imply (11) because [START_REF] Park | Stability and robust stability for systems with a time-varying delay[END_REF].

ℎ 2 -𝜏 ℎ 2 -ℎ 1 𝜂 𝑇 2 (𝑡)Ψ 1 𝜂 2 (𝑡) + 𝜏 -ℎ 1 ℎ 2 -ℎ 1 𝜂 𝑇 1 (𝑡)Ψ 2 𝜂 1 (𝑡) = 𝜂 𝑇 (𝑡)Ψ𝜂(𝑡) < 0. Thus, Ψ is convex in 𝜏 ∈ [ℎ 1 , ℎ 2 ]
if the four LMIs ( 13) and ( 14) are feasible, thus, the following result is obtained.

Lemma 1. Let there exist 𝑛 × 𝑛-matrices 𝑃 > 0, 𝑅 𝑖 > 0, 𝑆 𝑖 > 0, 𝑖 = 0, 1, 𝑃 2 ,𝑃 3 , 𝑌 1 , 𝑌 2 , 𝑇 1 , 𝑇 2 , 𝑍 1 and 𝑍 2 such that the four LMIs ( 13) and ( 14) with notations given in (15) are feasible. Then the switched system (1) is exponentially stable for delay 𝜏 ∈ [ℎ 1 , ℎ 2 ] with a decay rate 𝛼.

AVERAGE DWELL-TIME FOR SWITCHED TIME-DELAY INTERVALS

In order to reduce the conservatism, one would like to use multiple LKF on smaller delay intervals. To achieves this, the system (1) is rewritten as a switching system:

ẋ(𝑡) = 𝐴𝑥(𝑡) + 𝑁 ∑ 𝑖=0 𝛽 𝑖 𝐴 1 𝑥(𝑡 -𝜏 𝑖 (𝑡)), ( 16
)
where 𝛽 𝑖 ∈ {0, 1} and ∑ 𝑁 𝑖=0 𝛽 𝑖 = 1, 𝑁 corresponds to the number of the time-delay intervals.

As multiple LKF are used for different delay intervals, arbitrary switching cannot be obtained. So a certain criteria of minimum dwell time will be considered to guarantee the global stability [START_REF] Liberzon | Switching in Systems and Control[END_REF].

Here, to simplify the problem, the time-delay interval is equally divided into 𝑁 smaller intervals. Each one is considered for a mode of the switched system. Then, the bounds of the intervals can be written as

[𝜈 𝑖 , 𝜈 𝑖+1 ], 𝑖 ∈ {1, ..., 𝑁 }, with 𝜈 0 = 0, 𝜈 1 = ℎ 1 , 𝜈 𝑁 +1 = ℎ 2 and ∪ 𝑁 𝑖=1 [𝜈 𝑖 , 𝜈 𝑖+1 ] = [ℎ 1 , ℎ 2 ].
The stability of the system (1) equals to the global stability of the one of ( 16). Two methods are usually applied for verifying stability of switched systems: one is to find the common LKF and the other is to use different LKFs but with minimum dwell time for each mode. The common LKF cannot always be found, and in case of exponential stability, only the same value of 𝛼 can be obtained, i.e., the performance in each mode cannot be taken into account and the performance cannot be better in the case of small time-delay. Not using a common LKF reduces the conservatism as each mode can have a completely different behavior. This reduction allows, in most cases, getting better performances or robustness. So, here the latter method is applied.

The stability analysis is performed by using different LKFs for each time-delay interval. In order to make the LKFs comparable, their structure are chosen as follows:

𝑉 𝑖 (𝑡, 𝑥 𝑡 , ẋ𝑡 ) = 𝑥 𝑇 (𝑡)𝑃 𝑖 𝑥(𝑡)

+ 𝑁 ∑ 𝑗=0 𝑡-𝜈𝑗 ∫ 𝑡-𝜈𝑗+1 𝑒 2𝛼𝑖(𝑠-𝑡) 𝑥 𝑇 (𝑠)𝑆 𝑗𝑖 𝑥(𝑠)𝑑𝑠 + 𝑁 ∑ 𝑗=0 (𝜈 𝑗+1 -𝜈 𝑗 ) -𝜈𝑗 ∫ -𝜈𝑗+1 𝑡 ∫ 𝑡+𝜃 𝑒 2𝛼𝑖(𝑠-𝑡) ẋ𝑇 (𝑠)𝑅 𝑗𝑖 ẋ(𝑠)𝑑𝑠𝑑𝜃 (17)
where 𝑖 represents the mode, 𝑉 𝑖 for 𝜏 𝑖 ∈ [𝜈 𝑖 , 𝜈 𝑖 + 𝑟].

These LKFs consider the delay over 𝑁 zones. Note that if set 𝜈 1 = ℎ 1 , 𝑁 = 1, then the previous LKFs (2) are recovered for a time-delay in

𝜏 (𝑡) ∈ [ℎ 1 , ℎ 2 ].
Therefore, to verify the stability of system (1) can be divided into the following two small problems:

Ψ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ Φ11 -(1 -τ )𝑋 Φ12 R0 -𝑌 𝑇 1 𝑍 𝑇 1 (𝜏 -ℎ1)𝑌 𝑇 1 (ℎ2 -𝜏 )𝑍 𝑇 1 𝑌 𝑇 1 -𝑍 𝑇 1 + Φ16 * Φ22 -𝑌 𝑇 2 𝑍 𝑇 2 (𝜏 -ℎ1)𝑌 𝑇 2 (ℎ2 -𝜏 )𝑍 𝑇 2 𝑌 𝑇 2 -𝑍 𝑇 2 + Φ26 * * 𝜙3 0 0 0 -𝑇1 * * * 𝜙4 0 0 0 * * * * 𝜙5 0 ( 𝜏 -ℎ1)𝑇1 * * * * * 𝜙6 0 * * * * * * 𝑇1 + 𝑇 𝑇 1 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ < 0, (12)
Ψ 1 = ⎡ ⎢ ⎢ ⎢ ⎣ Φ11 Φ12 R0 -𝑌 𝑇 1 𝑍 𝑇 1 (ℎ2 -ℎ1)𝑍 𝑇 1 𝑌 𝑇 1 -𝑍 𝑇 1 + Φ16 * Φ22 |𝜏 =ℎ 1 -𝑌 𝑇 2 𝑍 𝑇 2 (ℎ2 -ℎ1)𝑍 𝑇 2 𝑌 𝑇 2 -𝑍 𝑇 2 + Φ26 * * 𝜙3 0 0 -𝑇1 * * * 𝜙4 0 0 * * * * 𝜙6 |𝜏 =ℎ 1 0 * * * * * 𝑇1 + 𝑇 𝑇 1 ⎤ ⎥ ⎥ ⎥ ⎦ < 0, ( 13
) Ψ 2 = ⎡ ⎢ ⎢ ⎢ ⎣ Φ11 Φ12 R0 -𝑌 𝑇 1 𝑍 𝑇 1 (ℎ2 -ℎ1)𝑌 𝑇 1 𝑌 𝑇 1 -𝑍 𝑇 1 + Φ16 * Φ22 |𝜏 =ℎ 2 -𝑌 𝑇 2 𝑍 𝑇 2 (ℎ2 -ℎ1)𝑌 𝑇 2 𝑌 𝑇 2 -𝑍 𝑇 2 + Φ26 * * 𝜙3 0 0 -𝑇1 * * * 𝜙4 0 0 * * * * 𝜙5 |𝜏 =ℎ 2 (ℎ2 -ℎ1)𝑇1 * * * * * 𝑇1 + 𝑇 𝑇 1 ⎤ ⎥ ⎥ ⎥ ⎦ < 0. (14)
(1)Guarantee the stability of the system in each time-delay interval. Apply the Theorem 1 to each time-delay interval, the maximum value of exponential decrease decay, i.e. 𝛼 𝑖 , of each mode can be obtained.

(2)According to the theory of (Hespanha J.P. and Morse A.S., 1999) and [START_REF] Zhai | Stability analysis of switched systems with stable and unstable subsystems: an average dwell time approach[END_REF], an average dwell time has to be calculated to guarantee the global system stability.

Consider during a period of time as 𝑡 𝑎 as showed in the Fig. 1, 𝑁 ′ denotes the number of switching over the time of 𝑡 𝑎 . For generalizing the problem, we consider 𝑁 ′ ≥ 𝑁 .

t

𝑉 𝑉 1 𝑉 2 𝑉 3 𝑉 𝑁 ′ -1 𝑉 𝑁 ′ 𝑡 0 𝑡 1 𝑡 2 𝑡 3 𝑡 𝑁 ′ -1 𝑡 𝑁 ′ 𝑡 𝑎

Fig. 1. Switching time-delay intervals

Decreasing sequence in Fig. 1 implies that there exits 𝜇 𝑖 , 𝑖 ∈ {1, ..., 𝑁 ′ }, the smaller positive real such that:

𝑉 𝑖+1 (𝑡, 𝑥 𝑡 , ẋ𝑡 ) ≤ 𝜇 𝑖 𝑉 𝑖 (𝑡, 𝑥 𝑡 , ẋ𝑡 ). ( 18
)
This is implied by the following conditions:

𝑃 𝑖+1 ≤ 𝜇 𝑖 𝑃 𝑖 , 𝑒 -2𝛼𝑖+1𝜈𝑗 𝑆 𝑗(𝑖+1) ≤ 𝜇 𝑖 𝑒 -2𝛼𝑖𝜈𝑗 𝑆 𝑗𝑖 , 𝑒 -2𝛼𝑖+1𝜈𝑗 𝑅 𝑗(𝑖+1) ≤ 𝜇 𝑖 𝑒 -2𝛼𝑖𝜈𝑗 𝑅 𝑗𝑖 , (19)
with 𝑗 ∈ {1, 2, ..., 𝑁 ′ }.

If the LMIs conditions above are satisfied, one have: 𝑉 𝑖 (𝑡, 𝑥 𝑡 , ṫ𝑡 ) ≤ 𝑒 -2𝛼𝑖(𝑡-𝑡0) 𝑉 𝑖 (𝑡 0 , 𝑥 𝑡0 , ẋ𝑡0 ) for each interval of time-delay. According to these conditions, the smaller positive 𝜇 𝑖 can be calculated.

Notice that:

𝑉 𝑁 ′ (𝑡 𝑎 , 𝑥 𝑡𝑎 , ẋ𝑡𝑎 ) ≤ 𝑒 -2𝛼 𝑁 ′ (𝑡𝑎-𝑡 𝑁 ′ ) 𝑉 𝑁 ′ (𝑡 𝑁 ′ , 𝑥 𝑡 𝑁 ′ , ẋ𝑡 𝑁 ′ ) ≤ 𝜇 𝑁 ′ 𝑒 -2𝛼 𝑁 ′ (𝑡𝑎-𝑡 𝑁 ′ ) 𝑉 𝑁 ′ -1 (𝑡 𝑁 ′ , 𝑥 𝑡 𝑁 ′ , ẋ𝑡 𝑁 ′ ) ≤ ⋅ ⋅ ⋅ ≤ 𝜇 1 𝜇 2 ⋅ ⋅ ⋅ 𝜇 𝑁 ′ 𝑒 -2𝛼1𝑡1 𝑒 -2𝛼2(𝑡2-𝑡1) ⋅ ⋅ ⋅ 𝑒 -2𝛼 𝑁 ′ (𝑡𝑎-𝑡 𝑁 ′ ) 𝑉 1 (𝑡 0 , 𝑥 𝑡0 , ẋ𝑡0 ) (20) 𝜇 1 𝜇 2 ⋅ ⋅ ⋅ 𝜇 𝑁 ′ 𝑒 -2𝛼1𝑡1 𝑒 -2𝛼2(𝑡2-𝑡1) ⋅ ⋅ ⋅ 𝑒 -2𝛼 𝑁 ′ (𝑡𝑎-𝑡 𝑁 ′
) ≤ 1 is sufficient to prove the switching stability for a given period of time 𝑡 𝑎 showed in Fig. 1, so the minimum average dwell time for the system has to satisfy the following condition:

(𝛼 1 -𝛼 2 )𝑡 1 + (𝛼 2 -𝛼 3 )𝑡 2 + ⋅ ⋅ ⋅ + 𝛼 𝑁 ′ 𝑡 𝑎 ≥ 0.5𝑙𝑛(𝜇 1 𝜇 2 ⋅ ⋅ ⋅ 𝜇 𝑁 ′ ). (21)

NUMERICAL EXAMPLES

In this section, we use two numerical examples as the same as in (He Y. et al., 2004) to show our results. To compare with other methods, we here consider the stability conditions, i.e., in lemma (1), 𝛼 = 0.

Example 1. Consider the system (1) with

𝐴 = [ -2 0 0 -0.9 ] , 𝐴 1 = [ -1 0 -1 -1 ] . ( 22
)
As listed in the table Fig. 4, when ℎ 1 = 0, our result greatly enlarge the value of maximum upper-bound of time-delay for the system. While considering for the case of interval variable time-delay, the comparison with the methods of [START_REF] Jiang | On h control for linear systems with interval time-varying delay[END_REF] and [START_REF] He | Delayrange-dependent stability for systems with time-varying delay[END_REF] shows the merits of our methods. Note that when ℎ 1 = 0, the result coincides with [START_REF] Park | Stability and robust stability for systems with a time-varying delay[END_REF].

To show the merits of using several intervals of timedelay, we take the same example with two intervals as suggested in [START_REF] Yue | A piecewise analysis method to stability analysis of linear continuous/discrete systems with time-varying delay[END_REF]. For simplicity, two intervals here are considered, i.e., the total interval is averaged divided into two smaller ones. Then the result in table 4 can be obtained. It is obvious that using multiple LKF reduces conservatism compared with common LKF [START_REF] Yue | A piecewise analysis method to stability analysis of linear continuous/discrete systems with time-varying delay[END_REF] Example 2. Consider the system (1) with

𝐴 = [ 0 1 -1 -2 ] , 𝐴 1 = [ 0 0 -1 1 ] . ( 23
)
In the table Fig. 4, for given lower bounds ℎ 1 of variable time-delay, the corresponding maximum values of upper bounds ℎ 2 are given by considering different methods to guarantee the stability of the system. The comparison shows that our results are much less conservative than the other methods.

As we can find in the two examples, when the upper bound of time-delay becomes bigger, considering an unique interval, it becomes very difficult to find a LKF to guarantee the stability for the system. Even that Lemma 1 has greatly removed the conservation for the result, but for the greater delay, the result is not satisfied. The method of dividing the one interval into several ones enlarges the bounds of time-delays.

CONCLUSION

In this paper, the delay-dependent stability condition is given for interval LTI time-delay system based on the piecewise analysis method and LKF method. The timedelay interval is divided into several smaller ones and the system switches according to the bounds of each interval. Using multiple LKF, exponential stability of each mode is obtain and the global stability is guaranteed by the average dwell time. For sake of simplicity, the delay interval is equally divided, but for some special cases when the timedelay varies much around the bounds of the intervals, an oscillation will be brought to the system. So, some techniques for dividing the time-delay interval can be further developed.

  Allowable upper bound of ℎ 2 with given ℎ 1 is 𝛼 1 = 0.46, 𝛼 2 = 0.01, 𝜇 1 = 𝜇 2 = 15. According to the equation 21, the minimum average dwell time can be calculated for a given period of time.

	Method	ℎ 1	0	1	2	3	4	4.4697
	(Jiang X. and Han Q. L., 2005) ℎ 2 1.01 1.64 2.39 3.20 4.06	-
	(He Y. et al., 2007)	ℎ 2 1.34 1.74 2.43 3.22 4.07	4.47
	(Park P.G. and Ko. J.W., 2007) ℎ 2 1.86	-	-	-	-	-
	Lemma 1	ℎ 2 1.86 2.06 2.61 3.31 4.09	4.47
	Fig. 2. Allowable upper bound of ℎ 2 with given ℎ 1
	Method	𝜈 1	0	0.5	1.0	1.5	2.0
	(Yue D. et al., 2008) 𝜈 3 1.98 2.05 2.16 2.37 2.64
	our result	𝜈 3 2.08 2.34 2.58 2.81 3.05
	Fig. 3. Allowable upper bound with switching delays
	Method	ℎ 1	0	0.3	0.5	0.8	1	2
	(Jiang X. and Han Q. L., 2005) ℎ 2 0.67 0.91 1.07 1.33 1.50 2.39
	(He Y. et al., 2007)	ℎ 2 0.77 0.94 1.09 1.34 1.51 2.40
	Lemma 1	ℎ 2 1.06 1.24 1.38 1.6 1.75 2.58
	Fig. 4.

, but a certain minimum dwell time has to be considered. So, for calculating the minimum dwell time, we divide the interval into two smaller ones. As showed in the equation 21, there is a compromise between the value of 𝛼 and 𝜇. For this example, we consider the time-delay interval as [ℎ1, ℎ2] = [0, 2.0], the best solution


G Zheng

J-P Barbot

D Boutat

T Floquet

J.-P Richard

Causal observability of nonlinear time-delay systems with unknown inputs

This paper investigates the problem of causal observability of the states and unknown inputs of nonlinear time-delay systems. Using the theory of non-commutative rings, the nonlinear time-delay system is transformed into a suitable canonical form to solve the problem. A necessary and sufficient condition is given to guarantee the existence of a change of coordinate leading to such a form. Then, algorithms are provided to check the possibility of obtaining causal estimations of the states and unknown inputs for the studied systems.

I. INTRODUCTION

Observation and estimation are important issues in control theory. For linear time-delay systems, various aspects of the observability problem have been studied in the literature, using different methods such as the functional analytic approach ( [START_REF] Bhat | Modal characterizations of controllability and observability in time delay systems[END_REF]) or the algebraic approach ( [START_REF] Brewer | Linear systems over commutative rings[END_REF], [START_REF] Sontag | Linear systems over commutative rings: A survey[END_REF], [START_REF] Fliess | Controllability and observability of linear delay systems: an algebraic approach[END_REF]). The theory of non-commutative rings has been applied to analyze nonlinear time-delay systems firstly in [START_REF] Moog | The disturbance decoupling problem for time-delay nonlinear systems[END_REF] for the disturbance decoupling problem of nonlinear time-delay systems, and for observability of nonlinear time-delay systems with known inputs in [START_REF] Xia | Analysis of nonlinear time-delay systems using modules over non-commutative rings[END_REF].

In this algebraic framework, the left Ore ring of noncommutative polynomials defined over the field of meromorphic functions is used for the analysis of nonlinear timedelay systems, since the rank of a module over this ring is well defined and can be used to characterize controllability, observability and identifiability of nonlinear time-delay systems.

Based on the algebraic framework proposed in [START_REF] Xia | Analysis of nonlinear time-delay systems using modules over non-commutative rings[END_REF], the problem of causal and non-causal observability of the states and unknown inputs for nonlinear time-delay systems has been studied in [START_REF] Zheng | On finite time observation of nonlinear time-delays systems[END_REF] assuming that the systems can be transformed into a suitable canonical form. This condition is relaxed in this paper by introducing additional information depending on the available measurements. A necessary and sufficient condition is given to guarantee the existence of those new variables.

II. ALGEBRAIC FRAMEWORK

Consider the following nonlinear time-delay system:

   ẋ = f (x(t -iτ )) + s j=0 g j (x(t -iτ ))u(t -jτ ) y = h(x(t -iτ )) = [h 1 (x(t -iτ )), . . . , h p (x(t -iτ ))] T x(t) = ψ(t), u(t) = ϕ(t), t ∈ [-sτ, 0]
(1) where x ∈ W ⊂ R n denotes the state variables, u = [u 1 , . . . , u m ] T ∈ R m is the unknown admissible input and y ∈ R p is the measurable output. Without loss of generality, it can be assumes that p ≥ m. τ is the basic time delay, and it is assumed that all times delays are multiple times of τ . i ∈ S -= {0, 1, . . . , s}, a finite set of constant time-delays. f , g j and h are meromorphic functions 1 where

f (x(t -iτ )) = f (x, x(t -τ ), . . . , x(t -sτ )). ψ : [-sτ, 0] → R n and ϕ : [-sτ, 0] → R m denote
unknown continuous functions of initial conditions. In this work, it is assumed that the system (1) with u = 0 is locally observable, and admits a unique solution for initial conditions ψ and ϕ.

Based on the algebraic framework introduced in [START_REF] Xia | Analysis of nonlinear time-delay systems using modules over non-commutative rings[END_REF], let K be the field of meromorphic functions of a finite number of the variables from {x j (t -iτ ), j ∈ [1, n], i ∈ S -}. With the standard differential operator d, define the vector space E over K:

E = span K {dξ : ξ ∈ K}
which is the set of linear combinations of a finite number of one-forms from dx j (t -iτ ) with row vector coefficients in K. Introduce the backward time-shift operator δ defined by

δ i ξ(t) = ξ(t -iτ ), ξ(t) ∈ K, for i ∈ Z + (2)
and

δ i (a(t)dξ(t)) = δ i a(t)δ i dξ(t) = a(t -iτ )dξ(t -iτ ) (3)
for a(t)dξ(t) ∈ E, and i ∈ Z + . Let K(δ] denote the set of polynomials of the form

a(δ] = a 0 (t) + a 1 (t)δ + • • • + a ra (t)δ ra (4)
where a i (t) ∈ K. The addition in K(δ] is defined as usual, but the multiplication is given by

a(δ]b(δ] = ra+r b k=0 i≤ra,j≤r b i+j=k a i (t)b j (t -iτ )δ k (5)
Note that K(δ] satisfies the associative law and is a noncommutative ring (see [START_REF] Xia | Analysis of nonlinear time-delay systems using modules over non-commutative rings[END_REF]). However, it is proved that the ring K(δ] is a left Ore ring [START_REF] Ježek | Rings of skew polynomials in algebraical approach to control theory[END_REF], [START_REF] Xia | Analysis of nonlinear time-delay systems using modules over non-commutative rings[END_REF], which enables to define the rank of a module over this ring. Let M denote the left module over K(δ]: M = span K(δ] {dξ, ξ ∈ K}, where K(δ] acts on dξ according to (2) and [START_REF] Brewer | Linear systems over commutative rings[END_REF].

With the definition of K(δ], (1) can be rewritten in a more compact form as follows:

   ẋ = f (x, δ) + m i=1 G i u i (t) y = h(x, δ) x(t) = ψ(t), u(t) = ϕ(t), t ∈ [-sτ, 0] (6)
where f (x, δ) = f (x(t -iτ )) and h(x, δ) = h(x(t -iτ )) with entries belonging to K, G i = s j=0 g j i δ j with entries belonging to K(δ]. It is assumed that rank K(δ]

∂h ∂x = p, which implies that [h 1 , . . . , h p ]
T are independent functions of x and its backward time-shifts.

III. NOTATIONS AND DEFINITIONS

Observability definitions for nonlinear systems without delays have been given given in [START_REF] Hermann | Nonlinear controllability and observability[END_REF] and in [START_REF] Diop | Nonlinear observability, identifiability and persistent trajectories[END_REF]. A definition of observability for time-delay systems is given in [START_REF] Marquez-Martinez | Observability and observers for nonlinear systems with time delays[END_REF]. A more generic definition is stated here as follows:

Definition 1: System (1) is locally observable if the state x(t) ∈ W ⊂ R n can be expressed as a function of the output and its derivatives with their backward and forward shifts as follows:

x(t) = α(y(t -jτ ), . . . , y (k) (t -jτ )) (7)
for j ∈ Z and k ∈ Z + . It is locally causally observable if [START_REF] Floquet | Super twisting algorithm based stepby-step sliding mode observers for nonlinear systems with unknown inputs[END_REF] is satisfied for j ∈ Z + and k ∈ Z + , and locally non-causally observer if [START_REF] Floquet | Super twisting algorithm based stepby-step sliding mode observers for nonlinear systems with unknown inputs[END_REF] is satisfied for j ∈ Z and k ∈ Z + . Following the same principle of Definition 1, the following definition for the reconstruction of the unknown inputs is given.

Definition 2: The unknown input u(t) can be estimated if it can be written as a function of the derivatives of the output and its derivatives with backward and forward shifts, i.e.

u(t) = β(y(t -jτ ), . . . , y (k) (t -jτ )) (8)
for j ∈ Z and k ∈ Z + . It can be causally estimated if ( 8) is satisfied for j ∈ Z + and k ∈ Z + , and non-causally estimated if ( 8) is satisfied for j ∈ Z and k ∈ Z + . Definition 3:

(Unimodular matrix) [14] A matrix A ∈ K n×n (δ] is said to be unimodular over K(δ] if it has a left inverse A -1 ∈ K n×n (δ], such that A -1 A = I n×n . Definition 4: (Change of coordinate) [14] z = φ(δ, x) ∈ K n×1 is a causal change of coordinate over K for system (1) if there exist locally a function φ -1 ∈ K n×1 and some constants c 1 , • • • , c n ∈ N such that diag{δ ci }x = φ -1 (δ, z).

The change of coordinate is bicausal over

K if max{c i } = 0, i.e. x = φ -1 (δ, z).
Note that the relative degree for nonlinear systems without delays is well defined via the Lie derivative (see [START_REF] Isidori | Nonlinear control systems (3rd edition)[END_REF]). Then many efforts have been done to extend the classical Lie derivative for nonlinear time-delay systems. Several researchers tried to extend the Lie derivative to nonlinear time-delay systems (see [START_REF] Germani | An asymptotic state observer for a class of nonlinear delay systems[END_REF], [START_REF]A new approach to state observation of nonlinear systems with delayed output[END_REF], [START_REF] Oguchi | Sliding-mode control of retarded nonlinear systems via finite spectrum assignment approach[END_REF], [START_REF] Oguchi | Input-output linearization of retarded non-linear systems by using an extension of lie derivative[END_REF]) in the framework of commutative rings. In what follows the relative degree and the observability indices for nonlinear time-delay systems are defined by extending the Lie derivative in the algebraic framework of [START_REF] Xia | Analysis of nonlinear time-delay systems using modules over non-commutative rings[END_REF], from the non-commutative rings point of view.

Let f (x(t -jτ )) and h(x(t -jτ )) for 0 ≤ j ≤ s be an n and p dimensional vectors, respectively, with entries

f r ∈ K for 1 ≤ r ≤ n and h i ∈ K for 1 ≤ i ≤ p. Let ∂h i ∂x = ∂h i ∂x 1 , • • • , ∂h i ∂x n ∈ K 1×n (δ] (9)
where for 1 ≤ r ≤ n:

∂h i ∂x r = s j=0 ∂h i ∂x r (t -jτ ) δ j ∈ K(δ]
then the Lie derivative for nonlinear systems without delays can be extended to nonlinear time-delay systems in the framework of [START_REF] Xia | Analysis of nonlinear time-delay systems using modules over non-commutative rings[END_REF] as follows

L f h i = ∂h i ∂x (f ) = n r=1 s j=0 ∂h i ∂x r (t -jτ ) δ j (f r ) (10)
For j = 0, [START_REF] Hermann | Nonlinear controllability and observability[END_REF] is the classical definition of the Lie derivative of h along f . For h i ∈ K, define

L Gi h i = ∂h i ∂x (G i )
Using the above definition of Lie derivative, the relative degree can be defined in the following way:

Definition 5: (Relative degree) System (6) has relative degree

(ν 1 , • • • , ν p ) in an open set W ⊆ R n if, for 1 ≤ i ≤ p, the following conditions are satisfied : 1) for all x ∈ W , L Gj L r f h i (x, δ) = 0, for all 1 ≤ j ≤ m and 0 ≤ r < ν i -1; 2) there exists x ∈ W such that ∃j ∈ [1, m], L Gj L νi-1 f h i (x, δ) = 0. If for 1 ≤ i ≤ p, (1) is satisfied for all r ≥ 0, then we set ν i = ∞.
Since ( 6) is assumed to be locally observable when u = 0, then one can define, for (6), the so-called observability indices introduced in [START_REF] Krener | Ad f,g ), (ad f,g ) and locally (ad f,g ) invariant and controllability distributions[END_REF]. Let

F k := span K(δ] dh, dL f h, • • • , dL k-1 f h for 1 ≤ k ≤ n. It can be shown that the filtration of K(δ]- module satisfies F 1 ⊂ F 2 ⊂ • • • ⊂ F n .
Then define 6) is locally observable with u = 0. Reorder, if necessary, the output components of ( 6), such that

d 1 = rank K(δ] F 1 and d k = rank K(δ] F k -rank K(δ] F k-1 for 2 ≤ k ≤ n. Let k i = card {d k ≥ i, 1 ≤ k ≤ n} . Then (k 1 , • • • , k p )
rank K(δ] { ∂h1 ∂x , • • • , ∂L k 1 -1 f h1 ∂x , • • • , ∂hp ∂x , • • • , ∂L kp-1 f hp ∂x } = k 1 + • • • + k p = n Since rank K(δ] ∂h ∂x = p, the observability indices (k 1 , • • • , k p ) for (h 1 , • • • , h p )
are well defined, but the order may be not unique.

IV. CANONICAL FORM AND CAUSAL OBSERVABILITY

Hereafter are recalled some results given in [START_REF] Zheng | On finite time observation of nonlinear time-delays systems[END_REF]. Proposition 1: [START_REF] Zheng | On finite time observation of nonlinear time-delays systems[END_REF] For 1 ≤ i ≤ p, denote k i the observability indices and ν i the relative degree index for y i of (6), and note ρ i = min {ν i , k i }. Then there exists a change of coordinate φ(x, δ) ∈ K n×1 , such that (6) can be transformed into the following form:

żi,j = z i,j+1 (11)
żi,ρi = V i = L ρi f h i (x, δ) + m j=1 L Gj L ρi-1 f h i (x, δ)u j ( 12
)
y i = C i z i = z i,1 (13) ξ = α(z, ξ, δ) + β(z, ξ, δ)u ( 14
)
where

z i = (z i,1 , • • • , z i,ρi ) T = h i , • • • , L ρi-1 f h i T ∈ K ρi×1 α ∈ K µ×1 , β ∈ K µ×1 (δ] with µ = n - p j=1 ρ j C i = (1, 0, • • • , 0) ∈ R 1×ρi Moreover if k i < ν i , one has V i = L ρi f h i = L ki f h i . For (11), note H(x, δ) = Ψ(x, δ) + Γ(x, δ)u (15) with H(x, δ) = (V 1 , • • • , V p ) T Ψ(x, δ) = L ρ1 f h 1 , • • • , L ρp f h p T and Γ(x, δ) =    L G1 ρ1-1 f h 1 • • • L Gm L ρ1-1 f h 1 . . . . . . . . . L G1 L ρp-1 f h p • • • L Gm L ρp-1 f h p    (16)
where

H(x, δ) ∈ K p×1 , Ψ(x, δ) ∈ K p×1 and Γ(x, δ) ∈ K p×m (δ]. Since Γ ∈ K p×m (δ] with m ≤ p, if rank K(δ] Γ = m,
there exists a matrix Ξ ∈ K p×p (δ] such that

ΞΓ = Γ 0
where Γ ∈ K m×m (δ] has full rank m. The set Φ is introduced as follows:

Φ = {dh 1 , • • • , dL ρ1-1 f h 1 , • • • , dh p , • • • , dL ρp-1 f h p } (17)
Then the following theorem can be stated. Theorem 1: [START_REF] Zheng | On finite time observation of nonlinear time-delays systems[END_REF] Consider the system (6) with outputs (y 1 , . . . , y p ) and the corresponding (ρ 1 , . . . , ρ p ) with ρ i = min{k i , ν i } where k i and ν i are the observability indices and the relative degree indexes, respectively. If

rank K(δ] Φ = n
where Φ is defined in [START_REF] Oguchi | Input-output linearization of retarded non-linear systems by using an extension of lie derivative[END_REF], then there exists a change of coordinate φ(x, δ) such that ( 6) can be transformed into [START_REF] Isidori | Nonlinear control systems (3rd edition)[END_REF][START_REF] Ježek | Rings of skew polynomials in algebraical approach to control theory[END_REF][START_REF] Krener | Ad f,g ), (ad f,g ) and locally (ad f,g ) invariant and controllability distributions[END_REF][START_REF] Marquez-Martinez | Observability and observers for nonlinear systems with time delays[END_REF] with dim ξ = 0.

Moreover, if the change of coordinate is locally bicausal over K, then the state x(t) of ( 6) is locally causally observable, and if Γ ∈ K m×m (δ] is also unimodular over K(δ], then the unknown input u(t) of ( 6) can be causally estimated as well. Due to the particular triangular structure of the system ( 11)-( 12), the whole state of the system, as well as the unknown input, can be estimated in finite time using higher order sliding mode observers [START_REF] Floquet | Super twisting algorithm based stepby-step sliding mode observers for nonlinear systems with unknown inputs[END_REF].

When rank K(δ] Φ < n, the problem is to estimate the state ξ. It may be still possible provided that some additional conditions are satisfied. This is detailed in the next Section.

V. EXTENDED CASE

In [START_REF] Barbot | An observation algorithm for nonlinear systems with unknown inputs[END_REF], a constructive algorithm to solve the above mentioned problem for nonlinear systems without delays has been studied. The result of this section can be seen as an extension of the work [START_REF] Barbot | An observation algorithm for nonlinear systems with unknown inputs[END_REF] to treat the observation problem for time-delay systems with unknown inputs of the form [START_REF] Fliess | Controllability and observability of linear delay systems: an algebraic approach[END_REF]. The objective is to generate additional variables from the available measurement and unaffected by the unknown input such that an extended canonical form similar to (11)-( 12) can be obtained for the estimation of the remaining state ξ.

For this, consider Φ defined in [START_REF] Oguchi | Input-output linearization of retarded non-linear systems by using an extension of lie derivative[END_REF]. If rank K(δ] Φ = j, one can select j linearly independent vectors over K(δ] from Φ, denoted as

Φ = {dz 1 , • • • , dz j } . Note £ = span R[δ] {z 1 , • • • , z j }
where R[δ] is the commutative ring of polynomials of δ with coefficients belonging to the field R and let £(δ] be the set of polynomials of δ with coefficients over £. Define the module spanned by the elements of Φ over £(δ] as follows

Ω = span £(δ] {ξ, ξ ∈ Φ} . ( 18
)
Define also

G = span R[δ] {G 1 , . . . , G m }
and its left annihilator

G ⊥ = span R[δ] {ω ∈ Ω | ωg = 0, ∀g ∈ G}.
Based on the above definitions, let us state the main result of the paper.

Theorem 2: Consider the system (6) with outputs y = (y 1 , • • • , y p ) T and the corresponding (ρ 1 , . . . , ρ p ) with ρ i = min{k i , ν i } where k i and ν i are the observability indices and the relative degree indexes, respectively. Suppose rank K(δ] Φ < n where Φ is defined in [START_REF] Oguchi | Input-output linearization of retarded non-linear systems by using an extension of lie derivative[END_REF]. There exist l new independent outputs over K suitable to the causal estimation problem if and only if

rank K H = l where H = span R[δ] {ω ∈ G ⊥ ∩ Ω | ωf / ∈ £}. ( 19
)
Moreover, the l additional outputs, denoted ȳi , 1 ≤ i ≤ l, are given by: ȳi = ω i f mod £

where ω i ∈ H. Proof:

Denote Q i = q i 1 , • • • , q i p
as a 1×p vector with q i j ∈ K(δ] for 1 ≤ j ≤ p. One has

Q i Γ = Q i    L G1 L ρ1-1 f h 1 • • • L Gm L ρ1-1 f h 1 . . . . . . . . . L G1 L ρp-1 f h p • • • L Gm L ρp-1 f h p    =      Q i      ∂L ρ 1 -1 f h1 ∂x . . . ∂L ρp-1 f hp ∂x           [G 1 , • • • , G m ]
because of the associativity law over K (δ]. Then according to the definition (9), one gets

Q i Γ = ω i [G 1 , • • • , G m ] = ω i G
where

ω i = n c=1 p j=1 q i j ∂L ρ j -1 f
hj ∂xc dx c . Moreover, one can check that

ω i f =      Q i      ∂L ρ 1 -1 f h1 ∂x . . . ∂L ρp-1 f hp ∂x           f = Q i           ∂L ρ 1 -1 f h1 ∂x . . . ∂L ρp-1 f hp ∂x      f      = Q i    L ρ1 f h 1 . . . L ρp f h p    = Q i Ψ.
According to [START_REF] Moog | The disturbance decoupling problem for time-delay nonlinear systems[END_REF], one has

Q i H = Q i (Ψ + Γu) = ω i f + ω i Gu ( 20
)
where

H = (V 1 , • • • , V p )
T is a vector which can be estimated in finite time but which is affected by the unknown input. Thus, one can generate l additional information suitable to solve the estimation problem if and only if one can find

l independent {Q 1 , • • • , Q l } vectors over K(δ] such that for each Q i = {q i 1 , • • • , q i p } satisfying q i j ∈ £(δ], one has Q i Γ = 0 and Q i H /
∈ £, Thus, one has to prove that the following conditions are equivalent:

1) there exist l row vectors [START_REF] Xia | Analysis of nonlinear time-delay systems using modules over non-commutative rings[END_REF]. Necessity: Suppose item 1) is satisfied, then according to [START_REF] Zheng | On finite time observation of nonlinear time-delays systems[END_REF], one has

Q i = q i 1 , • • • , q i p , with q i j ∈ £(δ], such that rank K(δ] {Q 1 , • • • , Q l } = l, Q i Γ = 0 and Q i H / ∈ £. 2) rank K H = l with H defined in
Q i Γ = ω i G = 0 ⇒ ω i ∈ G ⊥ and Q i H = ω i f / ∈ £. Thus rank K(δ] {Q 1 , • • • , Q l } = l implies rank K {ω 1 , • • • , ω l } = l. Since L ρi-1 f h i ∈ £ and q i ∈ £(δ], then ω i = n c=1 p j=1 q i j ∂L ρ j -1 f
hj ∂xc dx c ∈ Ω. Thus, ω i ∈ H defined in [START_REF] Xia | Analysis of nonlinear time-delay systems using modules over non-commutative rings[END_REF].

Sufficiency: Suppose rank K H = l. Then one can find l independent 1-forms over

K: {ω 1 , • • • , ω l } with ω i ∈ G ⊥ ∩ Ω which implies there exist l independent vectors over K(δ]: {Q 1 , • • • , Q l } with entries belonging to £(δ] such that rank K(δ] {Q 1 , • • • , Q l } = l, since for each Q i one has Q i Γ = ω i G = 0 and Q i H = ω i f / ∈ £. The variable ȳi = Q i H = ω i f mod £
can be used as an additional output since 1) H can be estimated in finite time;

2) Q i has entries in £(δ];

3) ȳi do not belong to the current set £ of measured variables.

Theorem 2 means that if rank K H = l, one can define an extended output y T , ȳ T ∈ R p+l with observability indices ki and relative degree indexes νi . Then, similarly to [START_REF] Oguchi | Input-output linearization of retarded non-linear systems by using an extension of lie derivative[END_REF], a new Φ can be defined. If rank K(δ] Φ = n, the state and the unknown inputs can be estimated in finite time. Otherwise, if rank K(δ] Φ < n and if Theorem 2 is still valid, the procedure is iterated until Theorem 1 is true, resulting in the causal observation of the states and unknown inputs for nonlinear time-delay systems of the form [START_REF] Fliess | Controllability and observability of linear delay systems: an algebraic approach[END_REF].

In the tables Algorithm 1 and 2, two algorithms are proposed in order to check the possible causal estimation of the states and unknown inputs. Algorithm 2 is used to determine whether the studied system can be transformed into the proposed canonical form, by using Theorem 2. If it is, Algorithm 1 then checks whether the states and unknown inputs of the system are causally observable.

Let us note that it is the locally bicausal change of coordinate which makes the state of system locally causally observable. And it is the unimodular characteristic of Γ over K(δ] which guarantees the causal reconstruction of unknown inputs. The non-causal case was studied in [START_REF] Zheng | On finite time observation of nonlinear time-delays systems[END_REF].

Algorithm 1 Check the causal estimation of the states and unknown inputs 1: Input: TDS with x ∈ R n , y ∈ R p and u ∈ R m . 2: Output: "Failed" or ("OK", x, û).

ν i ← RelaDeg(y i );
Get relative degree for y i 7:

k i ← ObsInd(y i ); Get observability indices for y i 8:

ρ i ← max{ν i , k i }; 9: end for 10: Φ ← {dh 1 , • • • , dL ρ1 f h 1 , • • • , dh p , • • • , dL ρp f h p }; 11: ρ ← {ρ i , • • • , ρ p }; 12: j ← rank K(δ
Φ ← {dz 1 , • • • , dz j };
Select j linearly independent vectors over K 7:

£ ← span R[δ] {z 1 , • • • , z j }; 8: Ω ← span £(δ] {ξ, ξ ∈ Φ}; 9: H ← span R[δ] {ω ∈ G ⊥ ∩ Ω | ωf / ∈ £}; 10: l ← rank K H; 11: if l > 0 then Theorem 2 is satisfied 12: ∃l 1-forms, s.t. H = span R[δ] {ω 1 , • • • ω l }; 13: Ȳ ← {ω i f mod £, 1 ≤ i ≤ l}; New outputs 14:
for each ȳi ∈ Ȳ do 15:

ν i ← RelaDeg(y i );

16:

k i ← ObsInd(y i ); 17: ρ i ← max{ν i , k i }; 18: ρ ← ρ ∪ {ρ i }; 19: Y ← Y ∪ {ȳ i }; 20: Φ ← Φ ∪ {dȳ i , • • • , dL
              ẋ1 = -δx 1 + δx 4 u 1 ẋ2 = -δx 3 + x 4 ẋ3 = x 2 -δx 4 u 1 ẋ4 = u 2 y 1 = x 1 y 2 = δx 1 + x 3 (21)
One can check that

ρ 1 = ρ 2 = 1 and Φ = {dx 1 , δdx 1 + dx 3 }whichgivesrank K(δ] Φ = 2 < n. Set G =span R[δ] {G 1 , . . . , G m }, one has G ⊥ = span R[δ] {dx 1 + dx 3 , dx 2 } Since rank K(δ] Φ = 2, £ = span R[δ] {x 1 , δx 1 + x 3 } and Ω = span £(δ] {ξ, ξ ∈ Φ} = span £(δ] {dx 1 , dx 3 } . One obtains Ω ∩ G ⊥ = span £(δ] {dx 1 , dx 3 } ∩span R[δ] {dx 1 + dx 3 , dx 2 } = span £(δ] {dx 1 + dx 3 }
From the definition of H in [START_REF] Xia | Analysis of nonlinear time-delay systems using modules over non-commutative rings[END_REF], one can check that rank K H = 1, which gives the one-form ω = dx 1 + dx 3 , satisfying ω ∈ Ω ∩ G ⊥ and ωf = -δx 1 + x 2 / ∈ £. Thus, a new output ȳ1 = h 3 is given by In this paper, causal observation of nonlinear time-delay systems with unknown inputs has been studied. The system has to be transformed in a canonical form suitable for the design of a finite time observer. A necessary and sufficient condition for the existence of a change of coordinate has been given. A constructive algorithm has been proposed to check the possibility to obtain causal estimations of the states and unknown inputs for the studied systems.

G.

  Zheng is with Project ALIEN, INRIA Lille-Nord Europe, 40, avenue Halley, 59650 Villeneuve d'Ascq, France gang.zheng@inria.fr J.-P. Barbot is with ECS ENSEA, 6 Avenue du Ponceau, 95014 Cergy-Pontoise, and Project ALIEN, INRIA, France barbot@ensea.fr D. Boutat is with ENSI de Bourges PRISME, 10 Boulevard de Lahitolle, 18020 Bourges, France driss.boutat@ensi-bourges.fr T. Floquet and J.-P. Richard are with LAGIS FRE CNRS 3303, Ecole Centrale de Lille, BP 48, 59651 Villeneuve d'Ascq, and Project ALIEN, INRIA, France {thierry.floquet, jean-pierre.richard}@ec-lille.fr

  are the observability indices and p i=1 k i = n since (

ȳ1 = h 3 = 1 Tx 1 = y 1 x 2 = y 3 x 3 = y 2 -δy 1 x 4 = -δ 2 y 1 + 1 0

 3113411 ωf mod £ = x 2 = δy 1 + (1 -δ) ẏ1 + ẏ2(22)which yields ν 3 = 2 and k 3 = 2, thus ρ 3 = 2. Finally one has the new Φ as follows:Φ = {dx 1 , δdx 1 + dx 3 , dx 2 , -δdx 3 + dx 4 }It can be checked that rank K(δ] Φ = 4, and the new £ is£ = span R[δ] {x 1 , δx 1 + x 3 , x 2 , -δx 3 + x 4 }This gives the following change of coordinatez = φ(x, δ) = (x 1 , δx 1 + x 3 , x 2 , -δx 3 + x 4 )T it is easy to check that it is bicausal over K(δ], sincex = φ -1 = z 1 , z 3 , z 2 -δz 1 , z 4 + δz 2 -δ 2 z δy 2 + ẏ3where the new output y 3 is defined in (22). Moreover, the matrix Γ with the new output y 3 is ΞΓ = I 2×2Consequently, according to Theorem 1, u 1 and u 2 are causally observable. A straightforward computation yieldsu 1 = ẏ1 -δ 3 y1+δ 2 y2+δ ẏ3 u 2 = -δ ẏ1 + ÿ3 VII. CONCLUSION

3 :

 3 Main: Check observability of states and inputs 4: Y ← {y 1 , • • • , y p } 5: for each y i ∈ Y do

	6:

  ] Φ; 13: result ← Search Diffeo(Φ, Y, ρ, j); 14: if result=("Failed") then Search the change of coordinate to transform TDS into the canonical form1: Function Search Diffeo(Φ, Y, ρ, j) 2: if j = n thenThe rank of Φ over K(δ] is equal to the dimension of the system

			Algorithm 2 3: return (Φ, Y, ρ);
			4: end if
			5: while j < n do
			6:
	15:	return ("Failed");
	16: end if
	17: (Φ, Y, ρ) ← result;
	18: z = φ(x, δ) ← (Φ, Y, ρ);	Get the change of
		coordinate
	19: if φ is bicausal over K then
	20:	x = φ -1 (δ, z) ∈ K n×1 ;	Causally observable
	21: else
	22:	x ← Null;	Non-causal observation
	23: end if
	24: Calculate Γ(x, δ)	According to Eq. (16)
	25: if rank K(δ] Γ = m then 26: ∃Ξ s.t. ΞΓ = Γ 0 27: if Γ is unimodular over K(δ] then with rank K(δ] Γ = m;
	28:	û = β(δ, y, • • • , y (n) )	Causal estimation
	29:	else
	30:	û ← Null;	Non-causal estimation
	31:	end if
	32: else
	33:	û ← Null;	Might not be estimated
	34: end if
	35: return ("OK", x, û);

i.e. quotients of convergent power series with real coefficients[START_REF] Conte | Nonlinear control systems: An algebraic setting[END_REF],[START_REF] Xia | Analysis of nonlinear time-delay systems using modules over non-commutative rings[END_REF].


Christian Beauger

Philippe Grosseau

Bernard Guilhot

Denis Huguenin

Philibert Iacconi

Low temperature synthesis of ZnGa2O4 : Mn

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