Hafiz Ahmed

Rosane Ushirobira

Denis Efimov

Damien Tran

Mohamedou Sow

Jean-Charles Massabuau

Automatic spawning detection in oysters: a fault detection approach

Keywords: Fault detection, signal processing, oysters, spawning detection

Using measurements of valve activity in a population of bivalves under natural environmental condition (16 oysters in the Bay of Arcachon, France), an algorithm for the automatic detection of spawning period of oysters is proposed. The algorithm is based on the fault detection approach and it works through the estimation of velocity of valves movement activity, which can be obtained by calculating the time derivative of the valves distance. A summarized description on the method used for the derivative estimation is provided, followed by the associated signal processing and decision making algorithm to determine spawning from the velocity signal. A protection from false spawning detection is also considered by analyzing the synchronicity in spawning. Through this study, it is shown that spawning in a population of oysters living in their natural habitat (i.e. in the sea) can be automatically detected without any human expertise, like visual analysis.

I. INTRODUCTION

Global awareness regarding our aquatic environment increases daily. This consciousness combined with the pressure from the general public, scientists and environmental protectionists and various other reasons lead to more and more strict rules and regulations for the protection of environment. In order to better protect the aquatic environment, a detailed analysis of the behavior of marine animals in their natural habitat could be essential. From this necessity, various monitoring approaches have already been proposed to observe the behavior of marines organisms. For example, Mussel Watch [START_REF] Goldberg | The mussel watch-a first step in global marine monitoring[END_REF], which is considered the first step towards a global marine monitoring. These monitoring methods made available a lot of data collected from the field regarding marine animals. Specifically, it opened up a new horizon for a more detailed biological and mathematical analysis of marine animals behavior that was not possible before.

Hafiz Ahmed, Rosane Ushirobira and Denis Efimov are with the Non-A team, Inria, Parc Scientifique de la Haute Borne, 40 avenue Halley, 59650 Villeneuve d'Ascq, France, {hafiz.ahmed, rosane.ushirobira, denis.efimov@inria.fr}. Damien Tran, Mohamedou Sow and Jean-Charles Massabuau are with Université Bordeaux 1, Centre National de la Recherche Scientifique (CNRS), Unité Mixte de Recherche (UMR) 5805 Environnement et Paléoenvironnement Océanique et Continental (EPOC), Arcachon, France and the EA team at UMR 5805 EPOC-OASU, Bordeaux, France, {d.tran, m.sow, jc.massabuau@epoc.u-bordeaux1.fr}. Denis Efimov is also with CRIStAL UMR 9189, Ecole Centrale de Lille, Avenue Paul Langevin, 59651 Villeneuve d'Ascq, France.

Hafiz Ahmed is partially supported by the regional council of Nord-Pas de Calais, France D. Efimov Bivalves, like oysters and mussels, are widespread throughout the world and represent a key element in the ecology and economy of coastal ecosystems. A large scale monitoring of bivalves has been taking place for quite a long time [START_REF] Scanes | Oyster watch: Monitoring trace metal and organochlorine concentrations in sydney's coastal waters[END_REF]. One such project is "MolluSCAN eye", a remarkable monitoring solution that records the valve activity of bivalves in various geographical location in Europe, Arctic and Tropics1 . The data obtained from this project makes possible the analysis of behavior of various bivalves from field measurements. These datasets lead to the successful study on various aspects of oysters behavior like the existence of rhythms [START_REF] Mat | Evidence for a plastic dual circadian rhythm in the oyster crassostrea gigas[END_REF], [START_REF] Mat | Looking for the clock mechanism responsible for circatidal behavior in the oyster crassostrea gigas[END_REF], [START_REF] Tran | Field chronobiology of a molluscan bivalve: How the moon and sun cycles interact to drive oyster activity rhythms[END_REF], growth rates [START_REF] Schwartzmann | In situ giant clam growth rate behavior in relation to temperature: A one-year coupled study of high-frequency noninvasive valvometry and sclerochronology[END_REF], behavioral response to harmful algae exposure [START_REF] Tran | Behavioral responses of crassostrea gigas exposed to the harmful algae alexandrium minutum[END_REF], [START_REF] Mat | Genetic and genotoxic impacts in the oyster crassostrea gigas exposed to the harmful alga alexandrium minutum[END_REF], growth rate recording [START_REF] Tran | Influence of sex and spawning status on oxygen consumption and blood oxygenation status in oysters crassostrea gigas cultured in a mediterranean lagoon (thau, france)[END_REF], ecological monitoring [START_REF] Ahmed | Dynamical model identification of population of oysters for water quality monitoring[END_REF] and so on.

In this work, we were interested in a particular behavior of oysters, i.e. their behavior during reproduction also known as spawning. Spawning is characterized by rapid and rhythmic contraction and relaxation of the valve to expel eggs in the water. Spawning observation is important in domains like aquaculture, ecology etc. In [START_REF] Ahmed | Dynamical model identification of population of oysters for water quality monitoring[END_REF], it was shown that the deviation of valve activity from normal behavior (i.e. slow and non rhythmic contraction and relaxation of valve) can be used as an automatic tool to suspect pollution in surrounding water. Spawning behavior is a deviation from normal behavior but not caused by pollution, so it is necessary to distinguish this particular behavior from any other behavior. The motivation of the current work comes from this necessity. The basis of our work is the valve activity dataset of the MolluSCAN eye project recorded in the Bay of Arcachon, France in 2007.

The goal of this paper is to develop an algorithm that can automatically detect spawning in oysters. In [START_REF] Galtsoff | Physiology of reproduction of ostrea virginica i. spawning reactions of the female and male[END_REF], P. S. Galtsoff did this detection by visual inspection under laboratory conditions. Moreover, due to various technical limitations (for example, unavailability of appropriate sensors in 1938), the recording of valve movements were limited in order to measure the data. This type of approach required constant attention, labor and experience to analyze and detect some spawning behavior as it occurs once or twice a year. An increasing number of online measurements for different populations of oysters makes this manual detection approach unpractical due to a big amount of data to be permanently analyzed by an expert. Therefore a solution is an automatic detection of the spawning. In this perspective, spawning can be then considered as a fault [START_REF] Zolghadri | Fault Diagnosis and Fault-Tolerant Control and Guidance for Aerospace Vehicles: From theory to application[END_REF]. This paper proposes the algorithm to detect that fault.

The outline of the paper is as follows: a brief description of the measurements scheme and experiments is given in section II (a more detailed information can be found in [START_REF] Sow | Water quality assessment by means of hfni valvometry and high-frequency data modeling[END_REF], [START_REF] Schmitt | Statistical study of bivalve high frequency microclosing behavior: Scaling properties and shot noise analysis[END_REF], [START_REF] Tran | Estimation of potential and limits of bivalve closure response to detect contaminants: application to cadmium[END_REF]). The processing of obtained data through section II is discussed in section III. Spawning behaviors of oysters are summarized in section IV while the detection of spawning is discussed in section V. Section VI contains the results and the conclusion of this work can be found in section VII.

II. MEASUREMENT SYSTEM DESCRIPTION

The monitoring site is situated in the Bay of Arcachon, France, at the Eyrac pier (Latitude: 44 • 40 N, Longitude: 1 • 10 W). Sixteen Pacific oysters, Crassostrea gigas, measuring from 8 cm to 10 cm in length were permanently installed on this site. These oysters were all from the same age group (1.5 years old) and came from the same local supplier. They also all grew in the Bay of Arcachon. They were permanently submerged on the sea bottom (at 3 m to 7 m deep in the water, depending on the tide activity).

The basic principles of the measurement process has been first described in [START_REF] Tran | Estimation of potential and limits of bivalve closure response to detect contaminants: application to cadmium[END_REF]. It was further modified later on in [START_REF] Chambon | Influence of the parasite worm polydora sp. on the behaviour of the oyster crassostrea gigas: A study of the respiratory impact and associated oxidative stress[END_REF]. A significant advantage of the developed monitoring system at sea lies in the fact that it is completely autonomous, without in situ human interference for at one full year. Each animal is equipped with two lightweight coils (sensors), ≈ 100 mg each, attached on the edge of each valve. These coils measure 2.5 × 2.5 × 2 mm and were coated with a resin sealing before attaching them on the valves. One of the coils sends a highfrequency sinusoidal signal, that is received by another coil. Measurements are performed every 0.1 sec successively (with the frequency 10 Hz) for one of the sixteen animals. This means that the behavior of a particular oyster is measured every 1.6 sec. Every day, 54000 triplets (1 distance, 1 stamped time value, 1 animal number) are collected for each oyster. The strength of the electric field produced between the two coils is proportional to the inverse of distance between the point of measurement and the center of the transmitting coil, which leads to estimation of the distance between coils.

On the shore, a second electronic unit takes care of the data acquisition and transmission. This unit is equipped with a GSM/GPRS modem and uses Linux operating system for driving the first control unit submerged in the water, managing the data storage with a time stamp, accessing the Internet, and transferring the data. The whole process is done by using an original self-developed software module and hardware that run with mobile phone technology. At the end of a working day, the collected data is transmitted to a central workstation server located in the Marine Station at Arcachon, France. The valve activity data is stored in a central database and the public access to this database is available via a website (under some restrictions on the amount of data) 2 . A schematic description of the monitoring system is presented in Fig. 1.

III. DATA PREPARATION

To summarize the description given in the previous section, after a specified duration of the experiment, the electronic system provides a matrix of values of opening of valves for each oyster P i,j ∈ [0, 12 mm] for 1 ≤ i ≤ N, 1 ≤ j ≤ n, where n = 16 is the number of oysters in the experiment and N is the number of measured points, N = 54000×N days where N days > 0 is the number of days of measurement. In our case, the dataset of 2007 collected at the Eyrac pier in the Bay of Arcachon has been used. So, N days = 365. For the population, i.e. for 16 oysters, the measurement system gave us total 315.36 million data points. The valve opening activity of three different oysters is illustrated in Fig. 2 where we can see that the valve distance activity is not uniform among the oysters. For example, for oyster N°7 the valve distance is moving around 2 -6 mm most of the time while for oyster N°12, it varies between 5-10mm. This diversity is originated by the difference of shell anatomy between animals and variable placement of sensors on valves. To compensate for this variability we normalized the data between certain bounds. Here, the normalization was done considering a 6 days window and with a bound [0, 1], where 0 represents complete closing of the valve while 1 represents complete opening of the valve. The formula to calculate the normalized value is:

p i,j = P i,j -min(P k,j ) max(P k,j ) -min(P k,j )
where k = i -N t , • • • , i with N t the total number of data points in last 6 days, P i,j is the actual valve distance and p i,j is the normalized valve distance. For further work, we have considered this normalized/scaled data instead of the original one [START_REF] Ahmed | Dynamical model identification of population of oysters for water quality monitoring[END_REF].

IV. SPAWNING OF OYSTERS

For oysters, the experimental study of spawning dates back to 1938, when American biologist Paul S. Galtsoff published his seminal work on the physiology of reproduction of oysters [START_REF] Galtsoff | Physiology of reproduction of ostrea virginica i. spawning reactions of the female and male[END_REF]. Based upon laboratory observation, spawning is a specific/particular type of shell/valve activity of female oysters [START_REF] His | Contribution à l'étude biologique de l'huitre dans le bassin d'Arcachon. Activité valvaire de Crassostrea angulata et de Crasssotrea gigas ; application à la reproduction de l'huitre japonaise[END_REF], [START_REF] Bernard | Ecologie de la reproduction de l'huitre creuse, Crassostrea gigas, sur les cotes atlantiques francaises: vers une explication de la variabilite du captage[END_REF]. In [START_REF] Galtsoff | Physiology of reproduction of ostrea virginica i. spawning reactions of the female and male[END_REF], an ostreograph was used to measure the valve activity and that data was later used to study spawning. Under normal environmental conditions shell/valve movements are characterized by long relaxation periods which may vary from a few minutes to hours and are often interrupted by secondary contractions. While during spawning (see Fig. 3, oyster 1 and 3), it can be seen that a series of rapid contractions and relaxations are occurring following one after another with remarkable regularity and continuing for about 30-40 minutes. Consistency in the amplitude of the relaxation curve, especially during the first half of the reaction and the remarkable rhythmicity of the contractions are the most distinctive features of the sexual reaction of a female. This phenomenon doesn't occur under any other circumstance. Burst of valve activity can be seen in other cases as well like under the influence of some external excitation (for example pollution or chemical injection) but (1) their frequency is never so regular, (2) will last for shorter period of time and (3) will have long relaxation period. It was also known that spawning propagates from one to another and eventually over a large fraction of the oyster community [START_REF] Galtsoff | Physiology of reproduction of ostrea virginica i. spawning reactions of the female and male[END_REF].

So, any rhythmic behavior to be considered as spawning should have certain characteristics. They are, for example, regularity in rhythm, consistency in amplitude, happening for 30 -40 minutes, short relaxation period, synchronization of spawning in the population, etc. However, spawning can happen with mild characteristics also (see Fig. 3, oysters 7 and 15). For example, instead of 30 -40 minutes duration, it can last 10 -20 minutes. In this work, we will focus on detecting any spawning behavior with strong characteristics or clearly distinguishable. So, the spawning behavior can be considered as a deviation from normal behavior. In Fault Detection literature, this is known as fault (i.e. deviation from normal behavior). Hence, the detection of this fault is equivalent to the detection of spawning.

In our experiment, we have collected the data of 16 oysters. From that dataset, we also found similar behavior as reported in [START_REF] Galtsoff | Physiology of reproduction of ostrea virginica i. spawning reactions of the female and male[END_REF]. In our case, the oysters happened to spawn on the 15th July, 2007. The data of 4 oysters including both spawning and non-spawning oysters can be seen in Fig. 3. A close look on oysters N°3 and N°15, can be seen in Fig. 4. In this figure, we can clearly see that oyster N°3 fulfills all the criteria to call its rhythmic behavior as female spawning according to [START_REF] Galtsoff | Physiology of reproduction of ostrea virginica i. spawning reactions of the female and male[END_REF] (regular rhythmicity, consistency in the amplitude during rhythmic behavior, happening for about 30 -40 minutes, contraction followed by relaxation or vice-versa with short relaxation period, etc.). However, if we look at oyster N°15, it is evident from its behavior that it is not spawning with distinguishable characteristics (like 30 -40 minutes duration, very regular contraction and relaxation, etc.) although some rhythmicity can be seen. So, oyster N°15 is not showing the type of spawning that we are interested in detecting in this work. Similarly, we can say that oyster N°1 is spawning and N°7 is not spawning in Fig. 3.

V. AUTOMATIC DETECTION OF SPAWNING Detection of spawning is totally done by visual analysis until now. In one of our previous works, we have tried to establish a relation between water quality and abnormal valve activity [START_REF] Ahmed | Dynamical model identification of population of oysters for water quality monitoring[END_REF]. There, we have showed that the deviation of valve activity from normal behavior, if it occurs in the whole animal group, can be used as an indicator for change in water quality. Since spawning behavior is a total deviation from normal behavior, according to [START_REF] Ahmed | Dynamical model identification of population of oysters for water quality monitoring[END_REF], it could also be considered as an indicator for change in water quality. However, in reality this is a totally normal behavior having little to do with the occurrence of poor water quality. So, by automatic spawning detection, we will be able to differentiate spawning behavior from numerous other abnormal motion. Moreover, it will save time and labor of visually analyzing the data to find the period of spawning. So, automatic detection of spawning can be useful in the above mentioned cases.

A typical pattern of spawning behavior (see Fig. 3) is a series of contraction and relaxation of valves. In other words, the velocity of valve movement is fluctuating in a regular rhythmic manner. Since the kinetic energy is related to the velocity, it will also fluctuate very rapidly and regularly. By passing this high frequency kinetic energy like a signal through a specially designed low-pass filter, we can expect to extract some features related to spawning. They can then be used for the detection of spawning.

The main idea of our spawning detection algorithm is first to calculate the velocity of valve movement. Velocity is nothing but the time derivative of the measured distance. Once we estimate the velocity, this information can be used to calculate the energy (square of velocity). Then by passing this energy through a low-pass filter and comparing the value with some pre-defined threshold, we can detect spawning which is a typical approach used for fault detection [START_REF] Zolghadri | Fault Diagnosis and Fault-Tolerant Control and Guidance for Aerospace Vehicles: From theory to application[END_REF]. We have chosen energy signal as the marker because it is easier to obtain a pattern from the square of a velocity signal (having higher amplitude during spawning) than the velocity itself because of the amplitude. The details of the spawning detection algorithm are given below.

A. Velocity Estimation

Since velocity is the time derivative of valve distance activity, we need to differentiate the valve distance to get the velocity. There exist various effective methods like [START_REF] Mboup | Numerical differentiation with annihilators in noisy environment[END_REF], [START_REF] Perruquetti | Finite-time observers: application to secure communication[END_REF], [START_REF] Efimov | A hybrid robust non-homogeneous finitetime differentiator[END_REF]. Out of various methods available, we have chosen the technique proposed in [START_REF] Mboup | Numerical differentiation with annihilators in noisy environment[END_REF] based on algebraic method because of its simplicity, effectiveness and noise compensation. A more detailed description of the approach can be found in [START_REF] Mboup | Numerical differentiation with annihilators in noisy environment[END_REF], [START_REF] Ushirobira | Algebraic parameter estimation of a multi-sinusoidal waveform signal from noisy data[END_REF]. According to this method for a real-valued signal y(t), analytic on some time interval, the first-order time derivative estimate for discrete time case can be written as:

ŷk = 6 M 2 T s M l=0 1 -2 l M y k-l (1)
where t k-l = (k -l)T s , y k-l = y(t k-l ), T s is the sampling time and M is the number of summation steps i.e. T = M T s .

B. Filtering of energy signal and spawning detection

From the normalized valve distance data p i,j , using equation (1) we can easily estimate the velocity v i,j . From the estimated velocity, the kinetic energy like signal can be calculated by taking square of the velocity signal as:

E i,j = v 2 i,j ,
where E i,j is the energy and v i,j is the velocity. This signal will be passed through the following low-pass like filter: where f i,j is filtered signal, T s = 1.6 sec. is the sampling period and µ, γ are the parameters of the filters. The filtered signal obtained through (2) will then be passed through another typical low-pass filter to eliminate the remaining high frequency fluctuation of the energy signal to differentiate the spawning and non spawning case. In some cases, the pattern of spawning and non-spawning signals are very similar (Oyster 8 and 12 in Fig. 8). Double low-pass filtering helps to differentiate in that case. The transfer function of a typical first order low pass filter has the following form:

f i+1,j = f i,j + T s (min(γ, |E i,j |) -µf i,j ) (2)
H(s) = 1 1 + s ωc (3)
where ω c is the cut-off frequency of the filter. This doubly filtered energy signal will then be used for the detection of the spawning. If the signal exceeds some threshold, it will be considered as spawning otherwise no spawning:

Spawning = 1, fi,j ≥ β 0, otherwise (4)
where fi,j is the doubly filtered energy signal and β is the threshold for spawning.

The flowchart of the process can be seen in Fig. 5. The parameters of the algorithm were tuned on a trial and error basis since establishing some tuning is very difficult because of the wide variations of behavior among different animals. The values of the parameters are omitted for space limitations.

VI. RESULTS AND DISCUSSIONS

As mentioned in Section IV, the spawning oysters can be identified by visual inspection of the data. After such an expert evaluation, the oysters N°1, 3, 4, 10, 11, 12, 13 and 16 were classified as clearly visible spawning oysters with spawning period of about 30 -40 minutes and having very regular rhythmic contraction and relaxation of valves. By this we will find the spawning oysters. Next the proposed algorithm has been applied. For example, the data of oysters N°1 and 5 can be seen in Fig. 6. According to the criteria mentioned in section IV, oyster N°1 is classified as a typical spawner while N°5 appears to exhibit a much lesser typical burst of contraction. Note that to be certain that 5 was really spawning or not would require the ability to check if it produces eggs or not, which was clearly out of the scope of the present work. The final output of our algorithm in this case can be seen in Fig. 7. From Fig. 6, it is clear that our algorithm can successfully detect the start and end of spawning. The result for the total population can be seen in Fig. 8. From Fig. 8, we can see that our algorithm can detect the spawning for all the clear cut spawning oysters. There is no false alarm and zero misdetection. This proves the effectiveness of our algorithm. However, in some cases we can see a little delayed detection. The delay appears from filtering and threshold selection mainly which is intentional to eliminate the chance of false detection. From Fig. 8, some abnormal contractions in the male populations can be seen also which is very similar to spawning. It can be said in this case that the male oysters are showing some kind of reaction to their female counterpart.

As mentioned in Section IV, according to [START_REF] Galtsoff | Physiology of reproduction of ostrea virginica i. spawning reactions of the female and male[END_REF], spawning propagates from one to another and then induces an increased burst activity that eventually spread to most individuals. This is a very important criteria to check if the oysters are really spawning, which may tolerate the effect of wrong detection of the algorithm for certain oysters. If they are really spawning, we would see the spawning propagating, otherwise not. One way to check this criteria is to calculate the number of oysters synchronously spawning at a given time. The number of synchronously spawning oysters in our case can be seen in Fig. 9 where we remark that the spawning spread from one oyster to another and then eventually spread to the entire community. At the peak of spawning, all female oysters were spawning. This is a very strong proof of spawning and also tells the effectiveness of our spawning detection algorithm. Moreover, we can see that it took very short time (approximately 7 minutes in this particular case) to propagate spawning from the first spawning oyster to the second one and much less for the rest of the oysters.

VII. CONCLUSIONS This paper presents an algorithm to automatically detect the start and end of the spawning period of a population of oysters, that is, to detect individual spawning events and the day of spawning. The developed algorithm is based on the estimation of velocity of valve movement. The algorithm was then tested on a population to validate its effectiveness. The obtained results are very promising and open up the scope of real-time spawning detection of oyster population in marine environment.

In future works, the algorithm can be tested on new measurement sites and different bivalves species. This work considered only strongly distinguishable spawning cases. Less evident cases should also be detected in future works. Moreover, delay improvement can also be done since there exists a gap in the detection for some oysters. Very noisy/faulty measurement detection system can also be included into the algorithm to prevent any false/early/late detection. Once a spawning day has been automatically identified, it is thus a solid working base for an expert to screen the records and classify the atypical bursts of contractions.

Figure 1 .

 1 Figure 1. [14] Synoptic representation of the system, from field to laboratory: (1) Sixteen oysters equipped with a couple of electrodes each; cables from oysters to 1st level electronic card; 1st level electronic card in a waterproof box (this whole set is permanently immersed); (2) electrical connection between the first and the second electronic cards; (3) second electronic card (emerged); (4) GPRS antennae; (5) GPRS and Internet connection; (6) Marine Station of Arcachon; (7) daily update on Internet

Figure 2 .

 2 Figure 2. Actual valve distance of three different oysters during 10 days

Figure 3 .

 3 Figure 3. Spawning and non-spawning behavior of oysters

Figure 4 .

 4 Figure 4. Behavior of oyster N°3 and 15 (zoomed)

Figure 5 .

 5 Figure 5. Flow chart of spawning detection process

Figure 6 .Figure 7 .

 67 Figure 6. Valve activity of oyster 1 and 5 (x-axis is in hour and y-axis is normalized valve distance)

Figure 8 .

 8 Figure 8. Spawning detection for the population

http://molluscan-eye.epoc.u-bordeaux1.fr/

http://molluscan-eye.epoc.u-bordeaux1.fr/


Hafiz Ahmed

Rosane Ushirobira

Denis Efimov

Wilfrid Perruquetti

On conditions of robust synchronization for multistable systems

This paper deals with a robust synchronization problem for multistable systems evolving on manifolds in the context of Input-to-State Stability (ISS) framework. Based on a recent generalization of classical ISS theory to multistable systems (may have unstable equilibrium), a robust synchronization protocol is designed with respect to a compact invariant set of the unperturbed system. The invariant set is assumed to admit a decomposition without cycles (basically no homoclinic nor heteroclinic orbits may exist). Numerical simulation examples illustrate our theoretical results.

I. INTRODUCTION

The first scientific observation regarding a synchronization event can be dated back to seventeenth century by a Dutch scientist C. Huygens, who observed a synchronous behavior in coupled pendulums [START_REF] Kuznetsov | Synchronization of two metronomes[END_REF]. This type of behavior can be seen in a variety of systems in various domains, for example biological, chemical, mechanical, etc. [START_REF] Blekhman | Synchronization in science and technology[END_REF], [START_REF] Osipov | Synchronization in Oscillatory Networks[END_REF], [START_REF] Pikovsky | Synchronization: A Universal Concept in Nonlinear Sciences[END_REF], [START_REF] Strogatz | Sync: How Order Emerges from Chaos in the Universe, Nature, and Daily Life[END_REF]. Synchronization of complex dynamical systems and/or network of systems has attracted a lot of attention from multidisciplinary research communities over the last decades because of their pervasive presence in nature, technology and human society. A collective behavior arises from the interconnection of dynamical systems and it has various potential application domains. For example, transient stability in power network [START_REF] Dorfler | Synchronization and transient stability in power networks and nonuniform kuramoto oscillators[END_REF], cooperative multitasking and formation control [START_REF] Wei | Distributed Consensus in Multi-vehicle Cooperative Control[END_REF]. The core of synchronization is the collective objective of agents in a network to reach a consensus about certain variables of interest.

The existing literature on the synchronization problem is very vast and covers many areas. In [START_REF] Gazi | Swarm Stability and Optimization[END_REF], the problem of formation control is investigated in swarms within the framework of output regulation in nonlinear systems. Detailed study regarding the control and synchronization of chaos can be found in [START_REF] González-Miranda | Synchronization and Control of Chaos: An Introduction for Scientists and Engineers[END_REF]. The paper [START_REF] Lewis | Cooperative Control of Multi-Agent Systems[END_REF] extends optimal control and adaptive control design methods to multi-agent nonlinear systems on communication graphs. Recent advances in various aspects of cooperative control of multi-agent systems can be found in The authors are with the Non-A team @ Inria, Parc Scientifique de la Haute Borne, 40 Hafiz Ahmed is partly supported by the regional council of Nord-Pas de Calais, France.

This work was supported in part by Région Nord-Pas de Calais, France, by the Government of Russian Federation (Grant 074-U01) and the Ministry of Education and Science of Russian Federation (Project 14.Z50.31.0031). [START_REF] Shamma | Cooperative Control of Distributed Multi-Agent Systems[END_REF]. The theoretical framework for design and analysis of distributed flocking algorithms can be found in [START_REF] Olfati-Saber | Flocking for multi-agent dynamic systems: algorithms and theory[END_REF].

In this work, we consider the synchronization problem for multistable systems based on the framework of ISS. This is a very well established technique for the study of stability and robustness of nonlinear systems. The ISS property provides a natural framework of stability analysis with respect to input perturbations (see [START_REF] Dashkovskiy | Input to state stability and allied system properties[END_REF] and references therein). The classical definition allows stability properties with respect to arbitrary compact invariant sets (and not simply equilibria) to be formulated and characterized. However, the implicit requirement that these sets should be simultaneously Lyapunov stable and globally attractive, makes the basic theory not applicable for a global analysis of many dynamical behaviors of interest, having multistability [START_REF] Angeli | Detection of multistability, bifurcations and hysteresis in a large class of biological positive-feedback systems[END_REF], [START_REF] Angeli | Multi-stability in monotone input/output systems[END_REF], [START_REF] Gelig | Stability of nonlinear systems with non unique equilibrium[END_REF], periodic oscillations [START_REF] Stan | Analysis of interconnected oscillators by dissipativity theory[END_REF] or hidden attractors [START_REF] Leonov | Hidden attractors in dynamical systems. from hidden oscillations in hilbert-kolmogorov, aizerman, and kalman problems to hidden chaotic attractor in chua circuits[END_REF], just to name a few, and only local analysis remains possible [START_REF] Chaves | Bistable biological systems: A characterization through local compact input-to-state stability[END_REF]. Attempts were made to overcome such limitations by introducing the notions of almost global stability [START_REF] Rantzer | A dual to Lyapunov's stability theorem[END_REF] and almost input-to-state stability [START_REF] Angeli | An almost global notion of input-to-state stability[END_REF], etc.

Recently, the authors in [START_REF] Angeli | On input-to-state stability with respect to decomposable invariant sets[END_REF] have proposed that the most natural way of relaxing ISS condition for systems with multiple invariant sets is equivalent to relax the Lyapunov stability requirement [START_REF] Efimov | Global lyapunov analysis of multistable nonlinear systems[END_REF] (rather than the global nature of the attractivity property). Using this relatively mild condition, the authors [START_REF] Angeli | On input-to-state stability with respect to decomposable invariant sets[END_REF] have generalized the ISS theory as well as, the related literature on time invariant autonomous dynamical systems on compact spaces [START_REF] Nitecki | Filtrations, decompositions, and explosions[END_REF] for multistable systems. In our current work, the results presented in [START_REF] Angeli | On input-to-state stability with respect to decomposable invariant sets[END_REF] are extended to provide sufficient conditions for the existence of robust synchronization for multistable systems.

The rest of the paper is organized as follows. Section II introduces some preliminaries about decomposable sets and the notions of robustness with respect to ISS. Our main results and the family of nonlinear systems being considered can be found in section III. Numerical simulation examples are given to illustrate these results in section IV. The paper closes with some concluding remarks in section V.

II. PRELIMINARIES

For an n-dimensional C 2 connected and orientable Riemannian manifold M without a boundary, let the map f (x, d) : M × R m → T x M be of class C 1 , and consider a nonlinear system of the following form:

ẋ(t) = f (x(t), d(t)) (1)
where the state x ∈ M and d(t) ∈ R m (the input d(•) is a locally essentially bounded and measurable signal) for t ≥ 0.

We denote by X(t, x; d(•)) the uniquely defined solution of (1) at time t fulfilling X(0, x; d(•)) = x. Together with (1) we will analyze its unperturbed version:

ẋ(t) = f (x(t), 0). (2)
A set S ⊂ M is invariant for the unperturbed system (2) if X(t, x; 0) ∈ S for all t ∈ R and for all x ∈ S. For a set S ⊂ M define the distance to the set |x| S = min a∈S δ(x, a) from a point x ∈ M , where the symbol δ(x 1 , x 2 ) denotes the Riemannian distance between x 1 and x 2 in M , |x| = |x| {0} for x ∈ M or the usual Euclidean norm of a vector x ∈ R n . For a signal d : R → R m the essential supremum norm is defined as

d ∞ = ess sup t≥0 |d(t)|. A function α : R + → R + is said to belong to class K, i.e., α ∈ K, if it is continuous, strictly increasing and α(0) = 0. Furthermore, α ∈ K ∞ if α ∈ K and it is unbounded, i.e. lim s→∞ α(s) = ∞.
For any points x, y ∈ M , the α and ωlimit sets for (2) can be defined as follows:

α(x) := {y | y = lim n→-∞ X(x, t n ) where t n -∞} ω(x) := {y | y = lim n→∞ X(x, t n ) where t n ∞} .

A. Decomposable sets

Let Λ ⊂ M be a compact invariant set for (2).

Definition 1. [24]

A decomposition of Λ is a finite and disjoint family of compact invariant sets Λ 1 , . . . , Λ k such that

Λ = k i=1 Λ i .
For an invariant set Λ, its attracting and repulsing subsets are defined as follows:

W s (Λ) = {x ∈ M : |X(t, x, 0)| Λ → 0 as t → +∞}, W u (Λ) = {x ∈ M : |X(t, x, 0)| Λ → 0 as t → -∞}.

Define a relation on the family of invariants sets of M

: for W ⊂ M and D ⊂ M , we write W ≺ D if W s (W) ∩ W u (D) = . Definition 2. [24] Let Λ 1 , . . . , Λ k be a decomposition of Λ, then 1. An r-cycle (r ≥ 2) is an ordered r-tuple of distinct indices i 1 , . . . , i r such that Λ i1 ≺ . . . ≺ Λ ir ≺ Λ i1 . 2. A 1-cycle is an index i such that [W u (Λ i ) ∩ W s (Λ i )] - Λ i = .
3. A filtration ordering is a numbering of the Λ i so that

Λ i ≺ Λ j ⇒ i ≤ j.
As we can conclude from Definition 2, the existence of an rcycle with r ≥ 2 is equivalent to the existence of a heteroclinic cycle for (2) [START_REF] Guckenheimer | Structurally stable heteroclinic cycles[END_REF]. And the existence of a 1-cycle implies the existence of a homoclinic cycle for (2) [START_REF] Guckenheimer | Structurally stable heteroclinic cycles[END_REF].

Let a compact set W ⊂ M be containing all α and ω limit sets of (2). Definition 3. The set W is called decomposable if it admits a finite decomposition without cycles, W = k i=1 W i , for some non-empty disjoint compact sets W i , which form a filtration ordering of W, as detailed in definitions 1 and 2.

B. Robustness notions

The following robustness notions for systems in (1) have been introduced in [START_REF] Angeli | On input-to-state stability with respect to decomposable invariant sets[END_REF]. Definition 4. We say that the system (1) has the practical asymptotic gain (pAG) property if there exist η ∈ K ∞ and a non-negative real q such that for all x ∈ M and all measurable essentially bounded inputs d(•) the solutions are defined for all t ≥ 0 and the following holds:

lim sup t→+∞ |X(t, x; d)| W ≤ η( d ∞ ) + q.
(3)

If q = 0, then we say that the asymptotic gain (AG) property holds.

Definition 5. We say that the system (1) has the limit property (LIM) with respect to W if there exists µ ∈ K ∞ such that for all x ∈ M and all measurable essentially bounded inputs d(•) the solutions are defined for all t ≥ 0 and the following holds:

inf t≥0 |X(t, x; d)| W ≤ µ( d ∞ ).
Definition 6. We say that the system (1) has the practical global stability (pGS) property with respect to W if there exist β ∈ K ∞ and q ≥ 0 such that for all x ∈ M and all measurable essentially bounded inputs d(•) the following holds for all t ≥ 0:

|X(t, x; d)| W ≤ q + β(max{|x| W , d ∞ }).
It has been shown in [START_REF] Angeli | On input-to-state stability with respect to decomposable invariant sets[END_REF] that to characterize (3) in terms of Lyapunov functions the following notion is appropriate: Definition 7. We say that a C 1 function V : M → R is a practical ISS-Lyapunov function for (1) if there exists K ∞ functions α 1 , α 2 , α and γ, and scalars q ≥ 0 and c ≥ 0 such that

α 1 (|x| W ) ≤ V (x) ≤ α 2 (|x| W + c),
the function V is constant on each W i and the dissipation inequality below holds:

DV (x)f (x, d) ≤ -α(|x| W ) + γ(|d|) + q.
If the dissipation inequality holds for q = 0, then V is said to be an ISS-Lyapunov function.

Notice that existence of α 2 and c follows (without any additional assumptions) by standard continuity arguments.

The main result of [START_REF] Angeli | On input-to-state stability with respect to decomposable invariant sets[END_REF] connecting these robust stability properties is states below: Theorem 8. Consider a nonlinear system as in (1) and let a compact invariant set containing all α and ω limit sets of (2) W be decomposable (in the sense of Definition 3). Then the following facts are equivalent.

1. The system admits an ISS Lyapunov function.

2.

The system enjoys the AG property.

3.

The system admits a practical ISS Lyapunov function. 4. The system enjoys the pAG property. 5. The system enjoys the LIM property and the pGS.

A system (1) for which this list of equivalent properties is satisfied, is called ISS with respect to the set W [START_REF] Angeli | On input-to-state stability with respect to decomposable invariant sets[END_REF].

III. SYNCHRONIZATION OF MULTISTABLE SYSTEMS

In this section we will consider the following family of nonlinear systems:

ẋi (t) = f i (x i (t), u i (t), d i (t)) , i = 1, . . . , N, N > 1, (4)
where the state x i ∈ M i and M i is an n i -dimensional C 2 connected and orientable Riemannian manifold without a boundary, the control u i (t) ∈ R mi and the external disturbance d i (t) ∈ R pi (u i (•) and d i (•) are locally essentially bounded and measurable signals) for t ≥ 0, the map f i (x i , u i , d i ) :

M i × R mi × R pi → T x M i is C 1 , f i (0, 0, 0) = 0. Denote the common state vector of (4) as x = [x T 1 , . . . , x T N ] T ∈ M = N i=1 M i (M is the corresponding Riemannian manifold of dimension n = N i=1 n i where the family (4) behaves), d = [d T 1 , . . . , d T N ] T ∈ R p with p = N i=1
p i is the total exogenous input. Assumption 1. For all i = 1, . . . , N , each system in (4) has a compact invariant set W i containing all α and ω limit sets of ẋi (t) = f i (x i (t), 0, 0), W i is decomposable in the sense of Definition 3, and the system enjoys the AG property with respect to inputs u i and d i as in Definition 4.

In this case, from Theorem 8, there exist

C 1 ISS-Lyapunov functions V i : M i → R with K ∞ functions α 1i , α 2i , α 3i , γ ui and γ di such that α 1i (|x i | Wi ) ≤ V i (x i ) ≤ α 2i (|x i | Wi + c i ), c i ≥ 0, (5) DV i (x i )f i (x i , u i , d i ) ≤ -α 3i (|x i | Wi ) + γ ui (|u i |) + γ di (|d i |)
for all i = 1, . . . , N . For a future reference, define also the invariant set of disconnected and unperturbed (

u i = d i = 0) family W = N i=1 W i ⊂ M (0 ∈ W). Then, by definition, there exist functions ν 1 , ν 2 ∈ K ∞ such that ν 1 (|x| W ) ≤ N i=1 |x i | Wi ≤ ν 2 (|x| W ) (6)
for all x ∈ M . Since the set W is compact, then there are functions ν 3 , ν 4 ∈ K ∞ and a scalar c 0 ≥ 0 such that for all x ∈ M ,

|x| ≤ ν 3 (|x| W ) + c 0 , |x| W ≤ ν 4 (|x|). (7)
Thus, in this work we will consider the family (4) under Assumption 1, i.e. a family of robustly stable nonlinear systems. In general, the sets W i include equilibrium (at the origin, for instance) and limit cycles of agents in (4). There exist many works devoted to synchronization and design of consensus protocols for such a family or oscillatory network [START_REF] Li | Consensus of multiagent systems and synchronization of complex networks: A unified viewpoint[END_REF], [START_REF] Zhao | Global bounded synchronization of general dynamical networks with nonidentical nodes[END_REF], [START_REF] Pogromsky | A partial synchronization theorem[END_REF]. The goal of this study is to find a condition under which the existence of a global synchronization/consensus protocol for d = 0 implies robust synchronization in (4) for a bounded d = 0.

Let a C 1 function y(x) : M → R q , y(0) = 0 be a synchronization measure for (4). We say that the family (4) is synchronized (or reached the consensus) if y(x(t)) ≡ 0 for all t ≥ 0 on the solutions of the network under properly designed control actions

u i (t) = ϕ i [y(x(t))] (8)
(ϕ i : R q → R mi is a C 1 function, ϕ i (0) = 0) for d(t) ≡ 0, t ≥ 0. In this case the set A = {x ∈ W : y(x) = 0} contains the synchronous solutions of the unperturbed family in ( 4) and the problem of synchronization of "natural" trajectories is considered since A ⊂ W. In fact, due to the condition ϕ i (0) = 0 the convergence of y (synchronization/consensus) implies that the solutions of the interconnection belong to W, the conditions of convergence of the synchronizing/consensus output y can be found in [START_REF] Li | Consensus of multiagent systems and synchronization of complex networks: A unified viewpoint[END_REF], [START_REF] Zhao | Global bounded synchronization of general dynamical networks with nonidentical nodes[END_REF], [START_REF] Pogromsky | A partial synchronization theorem[END_REF]. Therefore, it is assumed that the controls ϕ i (y) ensure the network global synchronization, while decomposability in general follows from Assumption 1. We are going to show that in the setup as above, by selecting the shapes of ϕ i , it is possible to guarantee robust synchronization of (4) for any measurable and essentially bounded input d.

To this end, by continuity arguments there exist functions η 1 , η 2 , µ i ∈ K ∞ with a scalar η 0 ≥ 0 such that for all x ∈ M :

|y(x)| ≤ η 0 + η 1 (|x| W ), |y(x)| ≤ η 2 (|x|), (9)
|ϕ i (y)| ≤ µ i (|y|)
(note that the first two inequalities are related through [START_REF] Wei | Distributed Consensus in Multi-vehicle Cooperative Control[END_REF]). Then the following intermediate result can be established under Assumption 1 for (4), [START_REF] Gazi | Swarm Stability and Optimization[END_REF].

In this paper proofs are omitted due to space limitation.

Proposition 9. Let Assumption 1 be satisfied for (4). Then there exist ϕ i , i = 1, . . . , N in ( 8) such that the interconnection ( 4), ( 8) has pGS property with respect to the set W.

To proof the above proposition, we need the following restriction to satisfy:

µ i (s) ≤ γ -1 ui N -1 α 4 • η -1 1 (0.5s) . (10)
Assumption 2. For µ i satisfying [START_REF] Lewis | Cooperative Control of Multi-Agent Systems[END_REF] and all ϕ i : R

+ → R + that satisfy |ϕ i (y)| ≤ µ i (|y|) , ∀y ∈ R p ,
the compact set A contains all α and ω limit sets of ( 4), ( 8) for d = 0, and it is decomposable.

Note that by definition of the set A, |x(t)| W ≤ |x(t)| A ≤ |x(t)| W + z for a scalar z ≥ 0 for all x ∈ M , then the pGS property with respect to the set A has also been proven.

Therefore, in the setup used in this work the boundedness of trajectories (boundedness of |x(t)| W implies the same property for |x(t)| according to [START_REF] Wei | Distributed Consensus in Multi-vehicle Cooperative Control[END_REF]) follows by a proper selection of the interconnection gain in [START_REF] Gazi | Swarm Stability and Optimization[END_REF], i.e. by decreasing the control gain a certain robustness of (4), ( 8) is inherited after individual systems as it is stated in Assumption 1.

Theorem 10. Let assumptions 1 and 2 be satisfied for (4), [START_REF] Gazi | Swarm Stability and Optimization[END_REF]. Then there exist ϕ i , i = 1, . . . , N in (8) such that the interconnection (4), [START_REF] Gazi | Swarm Stability and Optimization[END_REF] has AG property with respect to the set A.

Roughly speaking this qualitative result says that if the synchronized output y is related with |x| W as in [START_REF] González-Miranda | Synchronization and Control of Chaos: An Introduction for Scientists and Engineers[END_REF] and each system in the network is robustly stable as in Assumption 1, then the system can be robustly synchronized by a sufficiently small feedback proportional to y.

IV. EXAMPLES AND SIMULATIONS

A. Application to nonlinear pendulums with friction

This example is taken from [START_REF] Zhao | Global bounded synchronization of general dynamical networks with nonidentical nodes[END_REF]. Consider a network of nonlinear non-identical pendulums for i = 1, . . . , N , N > 1:

ẋ1i = x 2i , ẋ2i = -Ω 2 i sin(x 1i ) -κx 2i + d i , (11)
where the state x i = [x 1i , x 2i ] takes values on the cylinder M i := S × R, the exogeneous disturbance d i (t) ∈ R, κ is a constant parameter and Ω 2 i is the angular frequency of individual pendulums.The unperturbed system has a Hamiltonian H(x i ) = 0.5x 2 2i + Ω 2 i (1 -cos(x 1i )) and Ḣ = x 2i d i -κx 2 2i . Each unperturbed system has two equilibrium [0, 0] and [π, 0] (the former is attractive and the later one is a saddle-point), thus W i = {[0, 0] ∪ [π, 0]} is a compact set and containing all αand ω-limit sets of (11) for d i = 0. In addition, it is straightforward to check that W i is decomposable in the sense of Definition 3. Lemma 11. [START_REF] Efimov | Iss of multistable systems with delays: application to droop-controlled inverter-based microgrids[END_REF] For each i = 1, . . . , N , the systems in ( 11) is ISS with respect to the set W i .

Consequently, Assumption 1 is satisfied for (11) (since admitting an ISS Lyapunov function is equivalent to enjoying AG property according to Theorem 8) and we may select the synchronization measure y for the network. Since, in [START_REF] Zhao | Global bounded synchronization of general dynamical networks with nonidentical nodes[END_REF], authors have considered the first coordinate as synchronization measure, we are going to follow the same idea. The synchronization measure in this work is:

y = A sin(x 1 ),
where x 1 = [x 11 , . . . , x 1N ] T and A ∈ R N ×N is a Metzler matrix as in the first example.

Since the global boundedness of trajectories of (11) for bounded inputs is proven in Lemma 11, then a local analysis around equilibria is sufficient to show the synchronization measure convergence. It is straightforward to check that linearized around equilibria dynamics has y = 0 as a stable and attractive manifold. By this, the convergence of y is guaranteed locally. Then by taking,

ϕ i (y) = βy i , β > 0,
we may suppose that Assumption 2 is satisfied for some sufficiently small β. The results of simulations confirm this conclusion which can be seen in Fig. 1

V. CONCLUSIONS

Based on an extension of the ISS framework to systems with multiple invariant sets and evolving on (non-compact) a manifold, sufficient conditions for robust synchronization were derived. The condition imposed on the controller (ϕ i (0) = 0) made sure that the convergence of the synchronization measure implies that the interconnection (system-controller) belongs to the decomposable set W. Practical global stability analysis of the interconnection was done with respect to W. The asymptotic gain property of the interconnection with respect to the set of synchronous solutions A (A ⊂ W) was also proved.

Numerical simulations proved the effectiveness of our method to network of both identical and nonidentical nodes. However, our results are applicable only to systems that allow decomposition without cycles.

2 Figure 1 . 2 Figure 2 .

 2122 Figure1. The result of simulation for[START_REF] Shamma | Cooperative Control of Distributed Multi-Agent Systems[END_REF]

  avenue Halley, 59650 Villeneuve d'Ascq, France. Denis Efimov and Wilfrid Perruquetti are with CRIStAL UMR 9189, Ecole Centrale de Lille, Avenue Paul Langevin, 59651 Villeneuve d'Ascq, France. E-mail: {hafiz.ahmed, rosane.ushirobira, denis.efimov, wilfrid.perruquetti @ inria.fr} D. Efimov is with Department of Control Systems and Informatics, Saint Petersburg State University of Information Technologies Mechanics and Optics (ITMO), 49 av. Kronverkskiy, 197101 Saint Petersburg, Russia.


A Note on Continuous Delayed Sliding Mode Control

Keywords: Sliding mode, Time delay, Input-to-state stability

A new continuous sliding mode control approach is introduced with the dedicated mathematical tools. A timedelay modification/approximation of sign function is proposed, and it is shown that by substituting this new "sign" realization in the conventional sliding mode algorithms the main advantages of the sliding mode tools are preserved (like time of convergence and uniformity of stability with respect to matched disturbances), while the chattering is seriously reduced. These results are illustrated and confirmed by numerical simulations for the first order sliding mode control and the super-twisting algorithm.

I. INTRODUCTION

The sliding mode control and estimation algorithms became very popular nowadays due to their strong advantages: finite time of convergence and compensation of matched disturbances [START_REF] Perruquetti | Sliding Mode Control in Engineering[END_REF], [START_REF] Fridman | Sliding Mode Enforcement after[END_REF], [START_REF] Boiko | Discontinuous Control Systems: Frequency-Domain Analysis and Design[END_REF]. Despite of that the sliding mode solutions have also several drawbacks mainly originated by the impossibility of a perfect practical realization of sliding motion. Among these shortages it is necessary to mention the chattering phenomenon, which is a high frequency oscillation of control signal when trajectories stay around the sliding surface. The appearance of chattering may physically destroy actuator and/or degrade the performance of transients [START_REF] Utkin | Sliding Modes in Control Optimization[END_REF], [START_REF] Levant | Chattering analysis[END_REF]. There are several approaches for chattering reduction [START_REF] Fridman | Sliding Mode Enforcement after[END_REF]. One of the most popular deals with the high order sliding mode algorithms [START_REF] Levant | Sliding order and sliding accuracy in sliding mode control[END_REF], [START_REF] Edwards | Continuous higher order sliding mode control based on adaptive disturbance compensation[END_REF], where the discontinuity, for example sign[y(t)] = y(t) |y(t)| , y(t) = 0 0, y(t) = 0 , with y(t) ∈ R and |•| is the absolute value, appears not in the control/estimated variable itself but in its derivative of higher order [START_REF] Perruquetti | Sliding Mode Control in Engineering[END_REF], [START_REF] Fridman | Sliding Mode Enforcement after[END_REF]. However, since a discontinuity is present in the system closed loop, it may negatively influence the transients through the corresponding derivative (if this derivative has a physical meaning in the system). That is why frequently for practical realization of sliding mode algorithms different approximations of sign functions are used [START_REF] Ambrosino | Robust model tracking control for a class of nonlinear plant[END_REF], [START_REF] Burton | Continuous approximation of variable structure control[END_REF], like for instance ŝign[y(t)] = y(t)

+ |y(t)|
with some sufficiently small > 0, or see also [START_REF] Canudas-De Wit | Smoothing strategies for highgain control[END_REF] for an interesting general sigmoid approximation. Such a control based on approximated sign functions received the name of continuous sliding mode [START_REF] Khalil | Nonlinear Systems[END_REF], [START_REF] Oza | Continuous second order sliding mode based robust finite time tracking of a fully actuated biped robot[END_REF], [START_REF] Shtessel | An asymptotic second-order smooth sliding mode control[END_REF]. The main drawback of existent approximations is that a chattering reduction is achieved by a price of quality loss (a practical convergence in the presence of matched disturbances with an exponential rate [START_REF] Ambrosino | Robust model tracking control for a class of nonlinear plant[END_REF], [START_REF] Esfandiari | Stability analysis of a continuous implementation of variable structure control[END_REF]). In the present work a development of the continuous sliding mode control is presented that is based on a sign approximation using the time-delay framework, which guarantees the quality preservation (there is no static error in the presence of a matched disturbance and the speed of convergence is faster than exponential). The dynamical systems subjected by a time-delay, whose models are described by functional differential equations, find their applications in different areas of science and technology [START_REF]Applications of Time Delay Systems[END_REF]. Analysis of delay influence on the system stability is critical for many natural and human-developed systems [START_REF] Gu | Stability of Time-Delay Systems[END_REF], [START_REF] Kolmanovsky | Stability of functional differential equations[END_REF], [START_REF] Richard | Time-delay systems: an overview of some recent advances and open problems[END_REF]. The problem of stability investigation in timedelay systems is much more complicated than for ordinary differential equations since design of a Lyapunov-Krasovskii functional or a Lyapunov-Razumikhin function is a complex issue.

Usually appearance of a delay in the system leads to the performance degradation and augmentation of analysis complicacy. However in this work we are going to introduce delay in a proper way into the system, in order to make an approximation of the sign function, proving certain performance in the system under chattering reduction. The obtained approximation is very simple and can be easily implemented in control/estimation systems based on the sliding mode algorithms: In this work the new continuous delayed sliding mode control framework will be substantiated, it is based on such an approximation and admits the following advantages:

sign τ [y(t)] = y(
• the control algorithms demonstrate faster than exponential convergence to the origin; • exact cancellation of matched disturbances is preserved despite of approximation; • chattering reduction. The outline of this work is as follows. The preliminary definitions for time-delay systems are given in Section 2.

The motivating example of the first order sliding mode control algorithm is considered in Section 3. An extension to high order sliding mode algorithms is presented in Section 4. Several examples are considered for illustration of the obtained results. II. PRELIMINARIES Consider a functional differential equation of retarded type [START_REF] Kolmanovsky | Stability of functional differential equations[END_REF]:

dx(t)/dt = f (x t , d(t)), t ≥ 0, (1)
where

x(t) ∈ R n and x t ∈ C [-τ,0] is the state function, x t (s) = x(t + s), -τ ≤ s ≤ 0 (we denote by C [-τ,0] , 0 < τ < +∞ the Banach space of continuous func- tions φ : [-τ, 0] → R n with the uniform norm ||φ|| = sup -τ ≤ς≤0 |φ(ς)|, where | • | is the standard Euclidean norm) ; d(t) ∈ R m is an essentially bounded measurable input, i.e. ||d|| ∞ = sup t≥0 |d(t)| < +∞; f : C [-τ,0] × R m → R n is a locally bounded functional, f (0, 0) = 0.
The representation (1) includes pointwise or distributed time-delay systems with either constant or variable time delay τ (t) ∈ [0, τ ]. We assume that solutions of the system (1) satisfy the initial functional condition

x 0 ∈ C [-τ,0] .

A. Discontinuous functional differential equations

It is known from the theory of functional differential equations [START_REF] Kolmanovsky | Stability of functional differential equations[END_REF] that under the above assumptions the system (1) with a locally Lipschitz f has a unique solution x(t, x 0 , d) satisfying the initial condition x 0 for the input d(t), which is defined on some finite time interval [-τ, T ) (we will use the notation

x(t) to reference x(t, x 0 , d) if the origin of x 0 and d is clear from the context). If function f (φ, d) = f 0 (φ(0), φ, d)
, where f 0 : R 2n+m → R n , and it is discontinuous with respect to φ(0) and on a zero measure set N ⊂ R n only, following [START_REF] Filipov | Differential Equations with Discontinuous Righthand Sides[END_REF], [START_REF] Heemels | Input-to-state stability and interconnections of dicontinuous dynamical systems[END_REF] we will consider its multivalued extension (define

B ε (x) = {y ∈ R n : |x -y| ≤ ε} as a closed ball of radius ε > 0 around x ∈ R n ): F (φ, d) = ε>0 conv[f 0 {B ε (φ(0)) \ N , φ, d}],
which is non-empty, compact, convex and upper semicontinuous [START_REF] Heemels | Input-to-state stability and interconnections of dicontinuous dynamical systems[END_REF] for any d ∈ R m . In particular, the multivalued extension of sign τ [y(t)] can be defined as follows:

sign τ [y(t)] = y(t) max{|y(t)|,|y(t-τ )|} , y(t) = 0 [-1, 1], y(t) = 0 .
In this case, instead of (1), the solutions of the following functional differential inclusion will be considered:

dx(t)/dt ∈ F (x t , d(t)), t ≥ 0, (2)
and for any initial condition x 0 ∈ C [-τ,0] the set of corresponding solutions of (2) initiated at x 0 for the input d can be denoted as S(x 0 , d).

For a locally Lipschitz continuous function V : R n → R + the upper directional Dini derivative for (2) is defined as follows:

D + V [x t (0)]F (x t , d) = sup v∈F (xt,d) lim h→0 + sup V [x t (0) + hv] -V [x t (0)] h .
A continuous function σ : R + → R + belongs to class K if it is strictly increasing and σ (0) = 0; it belongs to class K ∞ if it is also radially unbounded. A continuous function r,•) is a strictly decreasing to zero for any fixed r ∈ R + .

β : R + × R + → R + belongs to class KL if β(•, r) ∈ K and β(

B. Stability definitions

Let Ω be an open neighborhood of the origin in C [-τ,0] , and denote D set of all essentially bounded measurable input d : R + → R m with ||d|| ∞ < D for some 0 < D < +∞. Definition 1. [START_REF] Moulay | Finitetime stability and stabilization of time-delay systems[END_REF] The system (2) for d = 0 is said to be (a) stable if there is σ ∈ K such that for any x 0 ∈ Ω and all x(t, x 0 , 0) ∈ S(x 0 , 0), |x(t, x 0 , 0)| ≤ σ(||x 0 ||) for all t ≥ 0;

(b) asymptotically stable if it is stable and lim t→+∞ |x(t, x 0 , 0)| = 0 for any x 0 ∈ Ω and all x(t, x 0 , 0) ∈ S(x 0 , 0); (c) finite-time stable if it is stable and for any x 0 ∈ Ω there exists 0 ≤ T x0 < +∞ such that x(t, x 0 , 0) = 0 for all t ≥ T x0 and all x(t, x 0 , 0) ∈ S(x 0 , 0). The functional T 0 (x 0 ) = inf{T x0 ≥ 0 : x(t, x 0 , 0) = 0 ∀t ≥ T x0 } is called the settling time of the system [START_REF] Perruquetti | Sliding Mode Control in Engineering[END_REF].

If Ω = C [-τ,0] , then the corresponding properties are called global stability/asymptotic stability/finite-time stability.

If the above properties hold for all x(t, x 0 , d) ∈ S(x 0 , d) with any x 0 ∈ Ω and any d ∈ D, then the corresponding properties are called uniform stability/asymptotic stability/finite-time stability.

For the forthcoming analysis we will need Lyapunov-Razumikhin theorem, which is given below (the variants for ordinary differential equations are given in [START_REF] Kolmanovsky | Stability of functional differential equations[END_REF], [START_REF] Gu | Stability of Time-Delay Systems[END_REF]).

Theorem 1. Let α 1 , α 2 ∈ K and η : R + → R + be a continuous nondecreasing function. If there exists a Lipschitz continuous function V : R n → R such that

α 1 (|x|) ≤ V (x) ≤ α 2 (|x|), ∀x ∈ R n ,
and the derivative of V along a solution x(t) of (2) satisfies

D + V [x(t)]F (x t , 0) ≤ -η(|x(t)|) if V [x(t + θ)] ≤ V [x(t)] (3)
for all θ ∈ [-τ, 0], then (2) is stable. If, in addition, η ∈ K and there exists a continuous nondecreasing function p(s) > s for s > 0 such that the condition (3) is strengthened to

D + V [x(t)]F (x t , 0) ≤ -η(|x(t)|) if V [x(t+θ)] ≤ p{V [x(t)]} for all θ ∈ [-τ, 0], then (2) is asymptotically stable. If in addition α 1 ∈ K ∞ , then (2
) is globally asymptotically stable. If the above properties are satisfied not only for d = 0 but for any d ∈ D, then (2) is uniformly stable/asymptotically stable.

Since the conditions of Theorem 1 are formulated for a Lyapunov function V , then the proof is the same for a differential equation or inclusion (it is omitted in this note).

Similarly, in order to analyze robustness of functional differential inclusions we will use an extension of the inputto-state stability (ISS) framework proposed in [START_REF] Teel | Connections between Razumikhin-type theorems and the ISS nonlinear small gain theorem[END_REF].

Definition 2. The system (2) is called ISS, if there exist β ∈ KL and γ ∈ K such that for all t ≥ 0 it holds

|x(t, x 0 , d)| ≤ β(||x 0 ||, t) + γ(||d|| ∞ )
for all x 0 ∈ C [-τ,0] , all essentially bounded measurable inputs d and all x(t, x 0 , d) ∈ S(x 0 , d). Definition 3. A locally Lipschitz function V : R n → R + is called ISS Lyapunov-Razumikhin function for (2), if there exist α 1 , α 2 ∈ K ∞ and χ τ , χ d , α ∈ K such that the following conditions hold:

α 1 (|x|) ≤ V (x) ≤ α 2 (|x|) ∀x ∈ R n ; V (x(t)) ≥ χ τ [ sup θ∈[-τ,0] V {x(t + θ)}] + χ d (|d(t)|) ⇒ D + V [x(t)]F (x t , d) ≤ -α(V [x(t)]),
for all solutions of (2) x(t, x 0 , d) ∈ S(x 0 , d) for all x 0 ∈ C [-τ,0] and all essentially bounded measurable inputs d. Theorem 2. [START_REF] Teel | Connections between Razumikhin-type theorems and the ISS nonlinear small gain theorem[END_REF] If there exists an ISS Lyapunov-Razumikhin function V for (2) and χ τ (s) < s for all s > 0, then the functional differential inclusion (2) is ISS from d to x with the gain γ(s

) = α -1 1 • χ d (s).

III. APPROXIMATION OF THE FIRST ORDER SLIDING MODE ALGORITHM

A. The rate of convergence

First of all, let us consider the system

ẋ(t) = -ksign τ [x(t)] (4)
for x 0 ∈ C [-τ,0] with τ > 0 and k > 0, in order to evaluate the rate of convergence in a system with the introduced approximation of the sign function, in this case

F (x t , 0) = -ksign τ [x(t)].
First of all note that for (4) the origin is the equilibrium: once x(t ) = 0 for some t ≥ 0, then x(t) = 0 for all t ≥ t . For any x 0 ∈ C [-τ,0] consider

D + V [x(t)]F (x t , 0) = -k x 2 (t) max{|x(t)|, |x(t -τ )|} = -kα(t)|x(t)| ≤ 0, α(t) = 1, |x(t)| ≥ |x(t -τ )| |x(t)| |x(t-τ )| , |x(t)| < |x(t -τ )| . Therefore, V (t) = V [x(t)]
is not increasing and the system is stable. Using this fact and introducing the sequence of

a i = |x(iτ )| for i ≥ 1, it is possible to show that a i+1 ≤ a 2 i a i + kτ , i ≥ 1. (5)
Moreover, the rate of convergence of this sequence is polynomial, and for small amplitudes of x(t), the state of ( 4) converges to zero with a hyper-exponential speed 1 . To confirm this conclusion, the results of this sequence calculation for different initial values of a 1 and for kτ = 1 are shown in Fig. 1 in logarithmic scale.

The following result can be proven (all proofs are excluded due to space limitations).

Lemma 1. The system (4) is globally asymptotically stable, and it has a hyper-exponential rate of convergence. Consequently, the proven rate of convergence in (4) is not finite-time, and by this it is slower than in the original system

ẋ(t) = -ksign[x(t)], (6)
but anyway it is much faster than an exponential rate, which can be obtained using standard approximations. Another consequence of this result is that in (4) there is no chattering of the right-hand side (that explains the paper title). For an illustration, the results of simulation of ( 4) and ( 6) with the discretization step h = 10 -4 are given in Fig. 2 for x(s) = 1, s ∈ [-τ, 0] and k = τ = 1. The solid red lines represent the variables x(t) and ẋ(t) generated by the system (4), the dash blue lines correspond to the variables of the system (6).

B. Uniform and robust stability/convergence

Now let us consider the perturbed version of (4):

ẋ(t) = -ksign τ [x(t) + d 1 (t)] + d 2 (t), (7)
1 I.e. the rate of convergence to the origin is much faster than exponential [START_REF] Polyakov | Implicit lyapunov-krasovski functionals for time delay systems[END_REF]. 4) and [START_REF] Levant | Sliding order and sliding accuracy in sliding mode control[END_REF] where d 1 (t) ∈ R and d 2 (t) ∈ R represent the measurement noise and the additive disturbance respectively. As usual in the sliding mode theory, we will be interested in the case when ||d 2 || ∞ < D and k > D, in this case

ẋ(t) ∈ F (x t , d 1 ), (8)
F (x t , d 1 ) = -ksign τ [x(t) + d 1 (t)] + [-D, D]. Consider for (8) a Lyapunov function candidate V [x(t)] = 1 2 x 2 (t): D + V [x(t)]F (x t , d 1 ) ≤ -kx(t)sign τ [x(t) + d 1 (t)] + |x(t)|D. Let √ 2b-1 √ 2bp -1 +1 k ≥ D +
for some > 0, where 0 < p < 1 and b > 0.5, and the following restriction is satisfied on trajectories of (8):

V [x(t)] ≥ b|d 1 (t)| 2 + b|d 1 (t -τ )| 2 + p 2 V [x(t -τ )] (9) ≥ max{b|d 1 (t)| 2 , b|d 1 (t -τ )| 2 , p 2 V [x(t -τ )]},
then the following relation can be proven under (9):

D + V [x(t)]F (x t , d 1 ) ≤ -2V [x(t)].
Thus V is an ISS Lyapunov-Razumikhin function for [START_REF] Ambrosino | Robust model tracking control for a class of nonlinear plant[END_REF] with respect to an extended input [d 1 (t) d 1 (t-τ )] T and according to Theorem 2 the functional differential inclusion [START_REF] Ambrosino | Robust model tracking control for a class of nonlinear plant[END_REF]

is ISS from [d 1 (t) d 1 (t -τ )] T to x(t): Lemma 2. Let √ 2b-1 √ 2bp -1 +1 k ≥ D + for some > 0, 0 < p < 1 and b > 0.5. Then (8) is ISS from d 1 (t) to x(t) with the gain γ(s) = 2 √ 2bs.
The results of simulation with the discretization step h = 10 -4 in the Euler method for [START_REF] Edwards | Continuous higher order sliding mode control based on adaptive disturbance compensation[END_REF] and the conventional sliding mode system

ẋ(t) = -ksign[x(t) + d 1 (t)] + d 2 (t) (10)
are shown in figures 3 and 4 for k = 1 and τ = 1 (in this case b = 10 3 , p = 0.96 and = 10 -2 for the given D = 0.9 and k = 1). In the top of Fig. 3 the variable x(t) is shown for both (7) (red line) and (10) (blue line), for d 1 (t) = 0 and d 2 (t) = 0.9 sin(5t). At the bottom of the figure the variable u(t) = sign τ [x(t) + d 1 (t)] for (7) (red line) and the signal -d 2 (t) (green dash line) are plotted (the variable 7) and [START_REF] Canudas-De Wit | Smoothing strategies for highgain control[END_REF] with d 1 (t) = 0.1 sin(10t) u(t) represents a kind of "control" approximation in a sliding mode control system) with its zoomed version for the last second of simulation. The variable u(t) = sign[x(t) + d 1 (t)] for [START_REF] Canudas-De Wit | Smoothing strategies for highgain control[END_REF] is not shown since it represents a pure chattering after approximately 1 second of time. As we can conclude, the system (7) has much less chattering and its signal u(t) estimates d 2 (t) directly without a need of filtering, which is usual for conventional sliding mode control and estimation systems. In Fig. 4 the results of simulation are shown for d 1 (t) = 0.1 sin(10t) and d 2 (t) = 0.9 sin(5t), the variables x(t) and u(t) are shown for (7) (red lines) and (10) (blue lines). In this case, again, the variable x(t) has a similar behavior for ( 7) and ( 10), but u(t) contains less chattering for [START_REF] Edwards | Continuous higher order sliding mode control based on adaptive disturbance compensation[END_REF].

IV. APPROXIMATION OF HIGH ORDER SLIDING MODE

ALGORITHMS

The main idea of the proof of Lemma 2 has the following interpretation: if for sign 0 [x(t) + d 1 (t)] the sliding mode system has an ISS Lyapunov function V , and since sign τ [x(t) + d 1 (t)] can be reduced to sign 0 [x(t) + d 1 (t)] under the Lyapunov-Razumikhin arguments, then V is an ISS Lyapunov-Razumikhin function for the system with sign τ [x(t) + d 1 (t)]. Thus a simple extension for high order sliding mode systems can be obtained.

Consider a conventional sliding mode system:

ẋ(t) ∈ G x(t), sign{s[x(t)]}, d(t) , (11)
where x(t) ∈ R n is the state, d(t) ∈ R m is the disturbance, G : R n+m+k → R n and s : R n → R k are continuous functions, s(t) = s[x(t)] ∈ R k is the sliding surface variable (application of the sign function to a vector argument is understood elementwise).

Assumption 1. Let the system (11) be ISS from d to x and there exist an ISS Lyapunov function V : R n → R n such that

α 1 (|x|) ≤ V (x) ≤ α 2 (|x|) ∀x ∈ R n ; V (x(t)) ≥ χ d (|d(t)|) ⇒ D + V [x(t)]G x(t), sign{s[x(t)]}, d ≤ -α(V [x(t)])
for some α 1 , α 2 ∈ K ∞ and χ d , α ∈ K.

Let us introduce another time-delayed modification of the sign function for the high order case:

sign τ [x(t)] =                s(t) max{|s(t)|,|s(t-τ )|} , V [x(t)] < pV [x(t -τ )]
and |s(t)| = 0

s(t) |s(t)| , V [x(t)] ≥ pV [x(t -τ )] and |s(t)| = 0 [-1, 1] k , |s(t)| = 0
for some constant 0 < p < 1, where V has been introduced in Assumption 1. As before, sign 0 [x(t)] = sign[x(t)].

Theorem 3. Let Assumption 1 be satisfied for the system [START_REF] Khalil | Nonlinear Systems[END_REF]. Then for any 0 < τ < +∞ the system

ẋ(t) ∈ G x(t), sign τ [x(t)], d(t) , (12)
is also ISS from d to x.

For an illustration of the result of this theorem consider the super-twisting differentiator:

ẋ1 (t) = -k 1 |x 1 (t)|sign[x 1 (t)] + x 2 (t), (13) ẋ2
(t) = -k 2 sign[x 1 (t)] + d(t), where 0 < k 1 < k 2 are parameters, G (x, s, d) = -k 1 |x 1 |sign[x 1 ] + x 2 -k 2 s + d .
As a Lyapunov function candidate take one from [START_REF] Moreno | Strict lyapunov functions for the supertwisting algorithm[END_REF] (the Lyapunov function proposed in [START_REF] Polyakov | Reaching time estimation for "supertwisting" second order sliding mode controller via lyapunov function designing[END_REF] also can be used): the blue dash lines correspond to the variables of ( 13), the variable u(t) in both cases demonstrate the function sign(•) behavior and its approximation. As we can conclude from this numerical experiment, the higher order sliding mode algorithm behavior is quite nice and smooth (see the variables x 1 (t) and x 2 (t)), but even in this case a chattering reduction is achieved by a slight decrease of the convergence rate using the new continuous delayed sliding mode algorithm.

V (x) = ξ T P ξ, where ξ = [ |x 1 (t)|sign[x 1 (t)]

V. CONCLUSIONS

In this note the Continuous Delayed Sliding Mode control has been presented, which demonstrates how by a delay introduction in the closed-loop of a nonlinear system it is possible to improve the quality of transients. An approximation of the sign function has been given based on delayed values of the processed variable, which ensures chattering reduction for the sliding mode control and estimation algorithms. The delay has been introduced in a way that under the Lyapunov-Razumikhin conditions the system is reduced to its original non-delayed form, then its stability can be easily concluded. Efficiency of the proposed approach has been demonstrated on examples: the chattering can be avoided (see Fig. 2) or tolerated (see figures 3-5) by a slight decreasing of the convergence rate, in addition, the matched disturbances are compensated.

Figure 1 .

 1 Figure 1. The results of the sequence (5) calculation

Figure 2 .

 2 Figure 2. Results of simulation for (4) and[START_REF] Levant | Sliding order and sliding accuracy in sliding mode control[END_REF]

Figure 3 .

 3 Figure 3. Trajectories of (7) and (10) for d 1 = 0

Figure 4 .

 4 Figure 4. Trajectories of (7) and (10) with d 1 (t) = 0.1 sin(10t)

x 2 ]Figure 5 .

 25 Figure 5. Results of simulation for super-twisting algorithm

  Efimov, A. Polyakov, W. Perruquetti and J.-P. Richard are with Non-A team @ Inria, Parc Scientifique de la Haute Borne, 40 av. Halley, 59650 Villeneuve d'Ascq, France and CRIStAL (UMR-CNRS 9189), Ecole Centrale de Lille, BP 48, Cité Scientifique, 59651 Villeneuve-d'Ascq, France. L. Fridman is with Departamento de Ingeniería de Control y Robótica, División de Ingeniería Eléctrica, Facultad de Ingeniería UNAM, Edificio T, Segundo piso Ciudad Universitaria D.F., México. D. Efimov is with Department of Control Systems and Informatics, Saint Petersburg State University of Information Technologies Mechanics and Optics (ITMO), 49 Kronverkskiy av., 197101 Saint Petersburg, Russia. This work was supported in part by collaboration program CNRS/CONACyT(209731), CONACyT grant 132125, PAPIIT-UNAM grant IN 113613, by the Government of Russian Federation (Grant 074-U01) and the Ministry of Education and Science of Russian Federation (Project 14.Z50.31.0031).

D.