Johan Sarrazin

Emmanuel Promayon

Michael Baumann

Jocelyne Troccaz

Hand-Eye Calibration of a Robot -UltraSound Probe System without any 3D Localizers

3D UltraSound (US) probes are used in clinical applications for their ease of use and ability to obtain intraoperative volumes. In surgical navigation applications a calibration step is needed to localize the probe in a general coordinate system. This paper presents a new hand-eye calibration method using directly the kinematic model of a robot and US volume registration data that does not require any 3D localizers. First results show a targeting error of 2.34 mm on an experimental setup using manual segmentation of five beads in ten US volumes.

I. INTRODUCTION

Today, 3D US volumes are increasingly used in clinical applications. The trade-off between the ease of use, the convenience, the image quality and the acquisition time of 3D volumes facilitate the development of new US guidance techniques in medicine. During medical procedures intraoperative US volumes can be used to register other intraoperative US volumes or to guide the surgeon during interventions. These US volumes can also be registered with preoperative MR or CT images where the definition of targets can bring valuable information.

However, the US probe alone does not provide any way of positioning each US volume in a single reference frame. In a majority of medical applications, the calibration of the probe is needed to register all data in the same mathematical space. To overcome this shortcoming the literature mainly presents two types of calibration techniques using external 3D localizers [START_REF] Bergmeir | Comparing calibration approaches for 3D ultrasound probes[END_REF]. In these techniques, a tracking reference frame is fixed on the probe and US volumes are acquired from a tracked phantom immersed in a water bath, as it explained in [START_REF] Bergmeir | Comparing calibration approaches for 3D ultrasound probes[END_REF], [START_REF] Baumann | 3-D Ultrasound probe calibration for computer-guided diagnosis and therapy. Proceedings of Computer Vision Approaches to Medical Image Analysis[END_REF] and [START_REF] Prager | Rapid Calibration for 3-D Freehand Ultrasound[END_REF], or by directly registering US volumes [START_REF] Schlosser | Automatic 3D ultrasound calibration for image guided therapy using intramodality image registration[END_REF].

During ultrasound guided prostate biopsy intervention, a well-distributed tissue sampling has to be performed using TransRectal UltraSound volumes (TRUS) in order to assess the presence of prostate cancer. A robot described in [START_REF] Poquet | An endorectal ultrasound probe comanipulator with hybrid actuation combining brakes and motors[END_REF] that holds the 3D US probe may be used to efficiently sample areas in the prostate. The calibration of this system could be an effective way to collect additional data to enable registrations. The calibration methods describe above have the drawbacks of requiring an external equipment such as tracking device and calibration rig. The reference frame fixed on the probe has to be detected by the tracking device which generally comes with a loss of movements' freedom. The use of large water basin or any other bulky device may also be problematic for day to day clinical use or involve additional electrical risks.

From these observations, a new method to calibrate a robot -US probe holder system by using the robot kinematic model and a 3D US volume registration algorithm was developed. Probe positions and orientations can be calibrated to the US volumes. In a previous work [START_REF] Sarrazin | 3D Ultrasound probe calibration using robotic arm and image registration[END_REF], a feasibility study of this method was performed. Preliminary results were obtained by simulating the robot via an optical Polaris localizer (NDI-Northern Digital Inc., Canada) and a tracker reference frame rigidly fixed on the 3D US probe. This initial study gave us a RMS targeting error of 1.83mm. It also allowed us to obtain a test bench to evaluate the accuracy of the robot compared to the existing tracking methods using optical localizers. Another interesting finding was that the manual segmentation of glass beads, serving as fiducials, in US volumes of average quality may induce a significant noise in the calibration result.

In this paper, we present the first results obtained by this hand-eye calibration method using the actual robot on a real experiment setup with a phantom. We evaluated the robustness of the method by performing a sensitivity study based on the impact of the segmentation noise on one hand, and the impact of the noise induced by both the robot sensors and by the registration method on the other hand. To describe this work, the method is presented in Section II. In the third section, the targeting error resulting during the experimental setup using the robot and the sensitivity study based on simulated data are detailed and discussed. Finally, concluding remarks are provided in the last section.

II. METHOD

The method presented in this paper adapts a standard handeye calibration algorithm [START_REF] Tsai | A new technique for fully autonomous and efficient 3D robotics hand/eye calibration[END_REF], [START_REF] Zhao | Hand-eye calibration using convex optimization[END_REF] to estimate the rigid transformation X from US volume frames to the robot coordinate system. A phantom is used to acquire all the necessary US volumes and robot positions. As it is represented in Fig. 1 we respectively use an US volume registration algorithm and the kinematic model of the robot to obtain the A i and B i matrices. A i are the transformations from the current 3D US volume to a reference US volume acquired at the beginning of the experiment. B i are the transformations from the current 3D US probe position to the robot coordinate system. The calibration matrix X is then obtained by solving the problem AX = XB with a minimum number of three measures with different rotations and translations.

Robot Base

US Reference volume

i th (i+1) th X X A i A i+1 B i+1 B i

A. Replacing the 3D localizers systems

In the standard calibration methods found in the literature a tracking system is used to obtain the position of the US probe in the localizer reference frame. Rigid trackers (usually called rigid bodies) are fixed on the 3D US probe. The rigid bodies provide the B i transformations (Fig. 1). Some methods track the phantom using another rigid body to obtain the transformation between the localizer and the phantom, i.e the A i transformation (Fig. 1). If the phantom is not tracked itself, then a registration method is used to register the US volumes between each other. Our method does not require any additional 3D localizer by using:

• a 3D US registration method, • a phantom specifically designed for the calibration,

• the kinematic model of a robot. The 3D US volume registration was used to compute the A i transformations. This registration method is provided by the Urostation (Koelis, La Tronche, France). It is the same algorithm used during TRUS biopsy session. The registration accuracy of this application is 0.8mm. A phantom made up of PolyVinyl Chloride (PVC) [START_REF] Hungr | A realistic deformable prostate phantom for multimodal imaging and needle-insertion procedures[END_REF] was built with a prostate casting from a segmented MRI, see Fig. 2 a). Its size is 16cm x 12cm x 9cm. The phantom was built with a target Young Modulus of 150kP a. For this experiments the same phantom is used for the calibration and the acquisitions. As a result, the speed of sound in PVC may not be taken into account for this evaluation. A series of clay beads with a diameter smaller than 2mm were added in the phantom in order to compute a targeting error after the calibration process (section II-B).

The kinematic model of a robot, developed by [START_REF] Poquet | An endorectal ultrasound probe comanipulator with hybrid actuation combining brakes and motors[END_REF], was used to obtain the B i transformations. The use of the robot instead of a 3D localizer is possible because the 3D US probe is rigidly fixed to the robot arm. The rotations and the translations of the probe are therefore taken into account in the robot's kinematic model. A reference position (orientation and translation) of the probe can be computed at any time in the robot coordinate system, which then can be used to solve the hand-eye problem.

B. Measuring the calibration error

To measure an error in the calibration process, the properties of our PVC phantom are used. The beads are manually segmented in the US volumes, Fig. 2 b), and then projected in the robot coordinate system thanks to the X matrix obtained by the hand-eye calibration and the B i matrices measured during the acquisitions. The Root Mean Square (RMS) targeting error is then computed as the distance of each projected bead to the mean position in the robot frame:

RM S = n i=1 m j=1 (Q j R -P j iR ) 2 n.m where Q j R = 1 n n i=1
P j iR is the mean position of the j th bead in the robot frame and P j iR is the estimated position of the j th bead in the i th US volumes transformed to the robot coordinate system. Note that P j iR = B i .X -1 . P j iU S , where P j iU S is the position of the bead j manually segmented in the US volume frame i.

C. Generation of virtual bench tests

To perform the sensitivity study, we developed a specific tool in the CamiTK framework [START_REF] Promayon | Using CamiTK for rapid prototyping of interactive Computer Assisted Medical Intervention applications[END_REF] to generate synthetic data sets and add some noise to them. Virtual reference frames symbolizing respectively the robot and the reference US volume are randomly generated. Then, a virtual tool composed of two reference frames, where the generated X transformation between these two coordinate systems is known, is displaced randomly in the scene. Each frame of the tool is computed in the virtual robot frame and in the virtual reference volume frame to obtain the A i and the B i matrices. A combination of Gaussian or Uniform noise can be added randomly to these virtual acquisitions in rotations and/or translations. Finally we solved the hand-eye problem to obtain the calibration matrix X and observe the impact of the noise on this matrix. This method mimics the experiments with the actual robot in terms of rotations and translations. To study the impact of the manual segmentation step, nine points are generated in the reference US volumes and computed in each new US volume. To have realistic data set, these points represent a cube and its centroid of a size similar to a classical TRUS volume. A Gaussian noise is added to these points to study the impact of the segmentation noise.

III. EXPERIMENTS

On the basis of the previous feasibility study, we decided to evaluate the calibration method with both synthetic and real data sets to respectively observe the impact of the noise in the calibration process and the ability of the robot to precisely track the probe movements. Each subsection presents the experimental setup, the results and the related discussion.

A. Impact of the segmentation

To study the impact of the segmentation error on RM S 12 virtual acquisitions were generated. Gaussian noise is added progressively to the synthetic point coordinates from σ = 0 to σ = 1 by step of 0.1. Fig. 3 presents the variation of RM S in function of the additional noise. The graph clearly shows that the quality of segmentation impacts the final calibration error. In our previous study, we observed that the error due to the manual segmentation was ranged from 1 to 1.5 voxels, where voxel size was 0.57mm x 0.57mm x 0.57mm. The maximum deviation error of a manually segmented bead was consequently between 0.99mm and 1.48mm. Fig. 3 shows that such a segmentation error is equivalent to a centered gaussian noise ranged from 0.5mm to 0.9mm STD in the final RM S .

Error RMS Min Error

If we crosscheck these results with our previous calibration results, the error induced by the segmentation can stand for up to 50% of the total RM S . It also means that with a perfect calibration the accuracy of the clinical application is highly dependent on the US volumes and segmentation quality.

B. Impact of the kinematic model and the image registration

To study the impact of the kinematic model and the image registration on RM S 150 set of 20 virtual acquisitions were generated with a translation motion range of [-5; 5]cm and a rotation range of [0; 60] o . Gaussian noise is added progressively to the 150 sets from σ = 0 to σ = 1 by step of 0.1, generating in total 33,000 acquisitions, each including a couple of noisy A i and B i . The noise mimics the volume registration noise in the A i transformations and the noise induced by the robot sensors in the B i transformations. For each of the 1,650 sets, RM S was computed depending on the number of acquisitions used in the hand-eye algorithm, using 4, 8, 12, 16 or 20 acquisitions. Fig. 4 presents the root mean square of the RM S errors in function of the additional noise level and the number of acquisition used for the calibration in the case of an exact segmentation. The graph shows that the noise in A i and B i transformations impacts the calibration results and that RM S also depends on the number of acquisitions used to solve the AX = XB system. Although using only 4 acquisitions for the calibration step would be faster in a medical application, the resulting error would be high. While using 8 acquisitions seems better, it shows more noise sensitivity with some local peaks. Using 12 or more acquisitions does not change significantly the resulting error. The difference range in RM S for the calibration performed with 12, 16 and 20 acquisitions is less than 0.12mm.

Twelve acquisitions seem to be a good compromise to obtain a robust and accurate result while minimizing the acquisition time required by this step.

C. Robot Experiments

For the first experiments using the robot probe holder and the dedicated phantom a series of 14 acquisitions were made, that includes 14 kinematic states of the robot and 14 US volumes. Acquisitions 1 to 12 were used to solve the AX = XB problem. Five beads were then manually segmented in eight US volumes used to solve the hand-eye algorithm (i.e., v1 to v8) in order to evaluate the accuracy of the process. Two other US volumes independent of the calibration process were also used to assess the robustness of the calibration (i.e., v9 and v10). Fig. 5 shows RM S for each segmented beads for volumes v 1 to v 10 . 5 shows that the maximum targeting error is given by the US volume v 5 . The errors can come from the registration step, the segmentation step or can be dependent to the experimental setup (for instance it can be due to the physical constraints applied to the end-effector of the robot). If we consider volume v 5 as an outlier, RM S can notably be improved to 2.12mm.

Compared to the literature, our results of 2.34mm with the targeting point error are very encouraging. [START_REF] Bergmeir | Comparing calibration approaches for 3D ultrasound probes[END_REF], [START_REF] Huang | Image registration based 3D TEE-EM calibration[END_REF] and [START_REF] Schlosser | Automatic 3D ultrasound calibration for image guided therapy using intramodality image registration[END_REF] respectively found a RMS targeting point error of 3.5mm by using 12 2D images, 2.37mm with 56 US volumes and 2.0mm with 36 US volumes. These hand-eye methods use accurate 3D localizers to obtain B i matrices and registration tools to obtain A i matrices. The comparison of these results with our sensitivity study shows that the RM S can rise sharply with a rotation and a translation noise (lesser than respectively σ = 0.5 o and σ = 0.5mm) induced by the robot sensors and by the registration tool. Compared to our previous work, the use of a robot instead of a 3D localizer generates an error of approximately 0.5mm. This error might mainly be induced by the use of the robot sensors but it is acceptable compared to the calibration results using a tracker device found in the literature.

On the basis of these results, other paths can be explored to improve our method such as the design of a better and more echogenic phantom to have the possibility to segment automatically shapes and points and to improve the registration. During medical interventions using a 3D US probe, the number of acquisitions is limited. As a consequence, the hand-eye calibration problem cannot be solved with a high number of positions to have a robust X matrix. Rejection criteria algorithms could so be implemented to eliminate acquisitions considered as outliers due to a bad registration [START_REF] Schlosser | Automatic 3D ultrasound calibration for image guided therapy using intramodality image registration[END_REF] or induced by a noise in the robot.

IV. CONCLUSIONS

We have proposed the first results of a new method to calibrate an ultrasound probe held by a robot by using the robot kinematic model and the image registrations as inputs of a hand-eye calibration problem. Unlike most of ultrasound probe calibration methods found in the literature, our method does not require any additional 3D localizer. Despite the use of less accurate systems compared to tracking devices, the resulting precision is good. Generated and real experiments have been performed to evaluate the performance of our method. Numerically, the calibration with 12 acquisitions seems to be enough as shown in our experiments. Finally, the trade-off between the ease of use of the method and the resulting RMS targeting error of 2.34mm is very encouraging.

Fig. 1 .

 1 Fig. 1. Application of the hand-eye calibration problem with a robot probe holder and a 3D US probe.

Fig. 2 .

 2 Fig. 2. (a) The specifically designed phantom contains glass beads and a look-alike prostate. (b) Beads and prostate are visible in the US images and can be manually segmented.

Fig. 3 .

 3 Fig. 3. RM S in function of the segmentation additional noise.

Fig. 4 .

 4 Fig.[START_REF] Schlosser | Automatic 3D ultrasound calibration for image guided therapy using intramodality image registration[END_REF]. Evolution of the calibration error when Gaussian noise is gradually added to the A i and the B i transformations.

Fig. 5 .

 5 Fig. 5. RM S for each segmented beads for volumes v 1 to v 10 .

TABLE I

 I RM S FOR FIVE SEGMENTED BEADS.

		RMS	Min	Max	Std
	All US volumes	2.34	0.87	3.82	0.73
	Without v5 US volume	2.12	0.84	3.31	0.61

RM S for this experiment is 2.34mm with a standard deviation of 0.73mm, see Table

I

. Fig

ACKNOWLEDGMENTS

The authors would like to thank M.A. Vitrani, P. Mozer, and D. Reversat for their technical and logistic support to use the robot as well as S. Voros and A. Agustinos for their help on the hand-eye calibration.

*This work is supported by Koelis, ANRT and ANR grants (PROSBOT project ANR-11-TECS-0017 and CAMI labex ANR-11-LBX-0004)


S Bejan

S Azou

P Morel

C Diouf

M Telescu

N Tanguy

A Sharaiha

A joint Linearization/Companding Approach for Improving a CO-OFDM Transmitter

Keywords: Coherent Optical OFDM, Semiconductor optical amplifier (SOA), Digital Predistortion, PAPR, Linearization

The joint use of peak-to-average power ratio (PAPR) reduction and linearization via digital predistortion is investigated in this letter, with the view to improve the performances of coherent optical OFDM (CO-OFDM) systems employing a semiconductor optical amplifier (SOA). PAPR reduction is performed via Wang's nonlinear companding transform (WNCT), which has been recently pointed out as a pertinent choice for optical communications, and a Filter Lookup Table (FLUT) scheme is considered for linearizing the transmitter. Experimental results prove the effectiveness of the proposed scheme, as a lower EVM is achieved with respect to system implementations using only PAPR reduction or linearization.

I. INTRODUCTION

OFDM is considered as a strong candidate for future Wavelength-Division-Multiplexing (WDM) systems due to its many advantages, such as high spectral efficiency, simple compensation of linear channel impairments, dynamic bandwidth allocation capability in a multiuser context, and powerful digital signal processing (DSP)-based implementation [START_REF] Cvijetic | OFDM for Next-Generation Optical Access Networks[END_REF]. It has been recently demonstrated in [START_REF] Khaleghi | Experimental Validation of Numerical Simulations and Performance Analysis of a Coherent Optical-OFDM Transmission System Employing a Semiconductor Optical Amplifier[END_REF] that SOA may be a pertinent choice for amplifying Coherent Optical-OFDM (CO-OFDM), for its interesting features such as large optical bandwidth, small size and possibility of integration at limited cost. Another study [START_REF] Schmuck | Demonstration of an SOA-assisted open metroaccess infrastructure for heterogeneous services[END_REF], conducted in an industrial context, demonstrates the practical feasibility of a SOA-based metro-access infrastructure, with the simultaneous operation of multiple services using different coherent modulation formats. Note also that SOA might be suitable for some low cost coherent systems [START_REF] Agmon | Bi-directional ultra-dense polarization-diverse OFDM/WDM PON with laserless colorless 1Gb/s ONUs based on Si PICs and <417 MHz mixed-signal ICs[END_REF] [START_REF] Dong | Monolithic silicon photonic integrated circuits for compact 100+ Gb/s coherent optical receivers and transmitters[END_REF]. However, the authors in [START_REF] Khaleghi | Experimental Validation of Numerical Simulations and Performance Analysis of a Coherent Optical-OFDM Transmission System Employing a Semiconductor Optical Amplifier[END_REF] report that SOA-based systems may suffer from nonlinear effects such as four wave mixing (FWM) and phase-amplitude coupling, as a result of large PAPR that is typically associated with OFDM signals. In a recent paper [START_REF] Azou | Performance Improvement of a SOA-based Coherent Optical-OFDM Transmission System via Nonlinear Companding Transforms[END_REF], we considered various companding schemes for improving a SOA-based CO-OFDM system performances. Thanks to an optimization of Wang's predistortion function [START_REF] Wang | Nonlinear Companding Transform for Reduction of Peak-to-Average Power Ratio in OFDM Systems[END_REF], we obtained a significant Error Vector Magnitude (EVM) reduction and a better power efficiency. This change in the dynamics of OFDM signal translates into less sensitivity to nonlinear channel impairments but the inherent nonlinear behavior of the SOA still significantly affects the performances.

In this letter we investigate the benefits of digital predistortion for linearizing the transmitter characteristics. Such techniques have been widely considered for radio amplifiers [START_REF] Guan | Green Communications: Digital Predistortion for Wideband RF Power Amplifiers[END_REF] but very few results have been reported for optical communication systems so far [START_REF] Vieira | Behavioral Modeling of Radio-over-Fiber Links Using Memory Polynomials[END_REF], [START_REF] Liu | Digital predistortion for RSOAs as external modulators in radio over fiber systems[END_REF], [START_REF] Hekkala | Predistortion of Radio Over Fiber Links: Algorithms, Implementation, and Measurements[END_REF]. An experimental validation of WNCT is addressed here for the first time. Moreover, we propose a joint use of companding and linearization, the latter being performed via a FLUT with a comparison to standard Lookup Table (LUT) implementation; to the knowledge of the authors, such a combination is completely novel.

II. SYSTEM MODEL

We first proceeded to the simulation of the CO-OFDM system with the objective of designing and comparing various candidate algorithms for PAPR reduction and linearization. Our implementation relies on a self-developed SOA model [START_REF] Khaleghi | Experimental Validation of Numerical Simulations and Performance Analysis of a Coherent Optical-OFDM Transmission System Employing a Semiconductor Optical Amplifier[END_REF], which has been fitted to simulate a commercially available bulk 750 µm long SOA (INPHENIX-IPSAD1501), so that it yields a very good matching between simulated results and experimental results. The overall structure of the system is illustrated in Fig. 1; it is a classical layout featuring a PAPR reduction block and a predistorter designed to reduce the nonlinear effects occurring between D/A input and SOA output. We recently investigated the issue of PAPR reduction in [START_REF] Azou | Performance Improvement of a SOA-based Coherent Optical-OFDM Transmission System via Nonlinear Companding Transforms[END_REF], with a special focus on nonlinear companding techniques. With such methods, the original signal samples x n have their amplitude modified according to a particular nonlinear function h, the resulting signal y n = h(x n ) is then converted into analog waveform. At the receiver side, the noisy signal r n = y n + v n is then transformed by the de-companding function in order to recover original signal (plus noise):

x n = h -1 (r n ) x n + h -1 (v n ).
In the present work, we consider that the crest factor reduction is achieved via WNCT [START_REF] Wang | Nonlinear Companding Transform for Reduction of Peak-to-Average Power Ratio in OFDM Systems[END_REF] as it has proved to be very effective in our previous study (balanced system performances considering criterions such as EVM, complexity and flexibility). Throughout the paper we consider the case of a 5 Gb/s QPSK/CO-OFDM transmission spread over 128 subcarriers, with a SOA biased at 200mA providing a gain of 19 dB (with a noise factor of 7 dB). The signal is sampled at 24 GS/s (oversampling factor of 8).

III. LINEARIZATION VIA DIGITAL PREDISTORTION A. General principle

Linearization methods aim at counteracting the nonlinear effects of a device under test (DUT; essentially the SOA in our case) in various ways so that the overall system appears to be linear [START_REF] Guan | Green Communications: Digital Predistortion for Wideband RF Power Amplifiers[END_REF]. A digital baseband predistortion (DPD) is typically adopted as it can compensate irregular characteristics and also efficiently operate over a wide dynamic range. Following this approach, we consider that the DUT is preceded by a DPD block (Fig. 1), denoted as the predistorter, aiming at reproducing the inverse DUT nonlinear behavior. Numerous DPD techniques have been proposed until now, with various performances and complexity. To this day there is no consensus as to the best possible scheme, even though a series of comparisons in terms of performance and complexity have been published, like the one in [START_REF] Tehrani | A comparative analysis of the complexity/accuracy tradeoff in power amplifier behavioral models[END_REF]. Choosing a particular approach is actually closely dependent upon the component to linearize and the application constraints.

LUT is a basic behavioral model for dealing with memoryless nonlinearity (NL). Following this approach a number of correction factors are pre-computed from curve fitting of the amplitude/amplitude (AM/AM) and amplitude/phase (AM/PM) measured characteristics of the DUT. The complex gain associated to the device is then stored in two LUTs, the output of the predistorter being expressed as

x out (n) = G(|x in (n)|)x in (n) where G(|x in (n)|)
is the instantaneous complex gain compensating the nonlinearity and x in , x out denoting the input and output signals of the predistorter. This compact implementation may be suitable for real-time applications, but due to quantization errors, the effectiveness is not always guaranteed. Moreover, the LUT scheme exhibits no memory effects, hence it cannot compensate the dynamical nonlinearities of the DUT.

B. The FLUT approach

In this work, we considered the FLUT scheme initially investigated in [START_REF] Jardin | Filter look up table method for power amplifiers linearization[END_REF] which consists in combining a LUT, for the memoryless NL effects, with a codebook of FIR filters (Fig. 2). Hence, the FLUT structure can be considered as an extension of a Hammerstein model in which a static NL is followed by a single linear filter. The LUT is filled with N 1 complex gain values, a particular gain G i(n) at index i(n) ∈ {0, 1, ..., N 1 -1} being selected according to the input signal magnitude |z(n)|, which is uniformly quantized. Hence, the output of the LUT block can be expressed as z l (n) = G i(n) z(n). The index i(n) further determines which filter in the codebook will be applied to the LUT distorted signal z l (n). For a given input magnitude index, the input-output relation characterizing the FLUT system is

z lf (n) = L-1 k=0 h i(n) (k)G i(n-k) z(n -k) (1)
where h j is the jth filter impulse response of length L.

Gain LUT Filter Codebook

Quantiz. In the FLUT system, learning is based on the combination of a direct approach the adaptation of the gains of the LUT and of an indirect approach for the adaptation of the filter impulse responses (see [START_REF] Jardin | Filter look up table method for power amplifiers linearization[END_REF] for details). We investigated the linearization via FLUT both in simulations and experimentally. As for every DPD scheme, the key parameters of the FLUT system (namely LUT size N 1 , filter codebook size N 2 , filter length L) need to be tuned so as to get a good performance improvement. The co-simulation scheme we used, based on an accurate SOA model, largely facilitated this tuning process. Note that the behavior of the linearization block alone is considered in this section, in order to clarify the influence of its key parameters, whereas joint tuning of linearization and PAPR reduction is investigated in the next section (practical experiments). Once the parameters {N 1 , N 2 , L} have been chosen, we extracted the FLUT coefficients according to the following three main steps: (1) At the first step, the DPD block is initialized as transparent at the transmitter side. The LUT is loaded with values equal to one, and the impulse responses of the FIR filters are unit impulses. With this initial setting, the output of the predistorter block is equal to its input. At the receiver side the coefficients of the LUT are computed via direct learning (aiming at directly minimizing the error between z(n) and the amplified signal, after some scaling by a reference gain G 0 ) and the performances of the system without linearization are evaluated; (2) For the second step, on the transmitter side, the LUT has the coefficients previously identified and the filter codebook still has unit impulse responses. On the receiver side, filters impulse responses are computed via indirect learning (in this case a post-distorter is first designed at the DUT output with no predistortion, and then it is used as a predistorter after removing the postdistorter) and a LMS algorithm, and the performance of the system using the LUT block only is evaluated. The adaptation step µ of the LMS algorithm should be correctly tuned, as it is well known that it can significantly affect the result (in this study we identified µ = 0.1 as a pertinent choice). ( 3) At the third step, the LUT is kept unchanged and the filter codebook is filled with the coefficients determined previously. The performances of FLUT can then be evaluated. Once the FLUT has been identified, it can be used permanently in the transmission chain until a major change occurs; then a new learning cycle needs to be performed. It should be emphasized that to ensure a proper operation of the inverse model (the FLUT system here), a one-to-one mapping between the input and the output must be guaranteed. Otherwise, the inversion of the nonlinear effects may not be possible. In case the SOA operates above its saturation threshold, the high peaks of the output signal will be clipped, which may result in inverse computation issues and predistorter divergence. This problem can be partly alleviated by appropriately choosing FLUT parameters, but, as it will subsequently be shown, combining FLUT with a PAPR reduction technique (WNCT here) is instrumental in successfully dealing with it. Another point to be mentioned is that the comparison of the performance of the linearized system with those obtained by the SOA alone should be conducted for the same average output power. In our study, we chose to operate with a reference gain G 0 = P out max /P in max with the maximum SOA output power P out max being very close to the saturation limit, P in max being the corresponding power at SOA input; consequently, the average input power of the linearized system must be larger than that of the SOA alone in order to obtain the same average output power.

The influence of the LUT size N 1 has been investigated for a codebook of 128 filters of length L = 5. Two different input power values were considered, -26 dBm (close to compression point, with an OSNR of 22.5 dB) and -17 dBm (in the saturation region, with an OSNR of 31.9 dB). The results (EVM against N 1 ) are depicted in Fig. 3.a for three different implementations, conventional CO-OFDM system (with no PAPR reduction), linearization via LUT only and linearization via FLUT. For the latter, the cases of a learning sequence and a test sequence are clearly distinguished for a fair evaluation. It can be clearly seen that the FLUT scheme offers better performance with respect to the LUT scheme for low N 1 values, best results being obtained for 32 bins. Above this threshold some convergence difficulties may occur. Around this optimum LUT size, similar EVM values are noticed. Then, we studied the influence of the filter codebook size N 2 , for N 1 = 32 (Fig. 3.b). As can be observed, for P in = -17 dBm, a local minimum in EVM is reached for N 2 = 16. For the unsaturated case, any change in N 2 has almost no impact on EVM, the three implementations giving close performances.

The required length L of the FIR filters is related to the memory effects intrinsic to the SOA. For identifying its value, we evaluated EVM for L ∈ [START_REF] Cvijetic | OFDM for Next-Generation Optical Access Networks[END_REF][START_REF] Vieira | Behavioral Modeling of Radio-over-Fiber Links Using Memory Polynomials[END_REF]. The results are illustrated in Fig. 3.c; a good value of the filter length appears to be L = 3, for both power values. A slight degradation in performance can be noticed above this value for the saturated case.

We also examined the impact of the signal bandwidth (BW) on the linearization performance (Fig. 4) for two operating points (P out = -1.3 dBm, P out = 3.1 dBm), by varying the OFDM symbol duration while keeping 128 subcarriers with an oversampling factor of 8. It can be seen that the SOA behavior remains the same up to 100 MHz, FWM being the main nonlinear effect in this range. Between 100 MHz and 1 GHz, the degradation in FLUT performances denotes an increased memory effect, and cross-phase/crossgain modulations (XPM/XGM) tend to appear. A BW above 1 GHz leads to lower amplitude-phase coupling and FWM effect.

IV. EXPERIMENTAL SETUP AND RESULTS

Figure 5 shows the experimental setup we considered. As described in the previous section, a preliminary task was the tuning of the FLUT predistorter, together with the WNCT by simulations (Fig. 1). The settings were chosen as follows: N 1 = 32, N 2 = 16, L = 3 and µ = 0.1 for FLUT; c = 0.43 and k 2 = -0.25 for WNCT. For the WNCT/FLUT combining, the setting was changed by taking N 2 = 512 and L = 2 (these values being identified in simulations in the same way as in section 3), the other coefficients remaining unchanged. The test sequence was a QPSK/OFDM signal with 128 subcarriers, 3GHz bandwidth, oversampled by a factor of 4, with a data rate of 5Gb/s. Due to memory limit, the total number of OFDM symbols transmitted for one run was limited to 19 (the chosen bandwidth also results from the same constraint), the first symbol being used as a preamble for acquisition; 11 subcarriers were set to zero to establish a guard band on the spectrum sides and 4 subcarriers near DC component were null subcarriers. Moreover, 8 subcarriers were used as pilots, the remaining 104 subcarriers being used to transmit data. The signals were generated in Matlab and uploaded to an Agilent M8190A AWG operating at 12 GSa/s; the I/Q components then passed through two RF amplifiers with 16 GHz bandwidth before entering the optical IQ modulator, with a half-wave voltage of about 5.1 V and a 1540 nm CW light applied at its input (ECL laser with 100 kHz linewidth). The transmission link includes only one SOA as a booster. At the receiver side, an Agilent N4391A Optical Modulation Analyzer (OMA) was used as a local oscillator for coherent detection. In order to validate the predistortion techniques, raw signal was recorded at 40 GSa/s and sent to a receiver implemented under Matlab; the main steps prior to FLUT learning and decompanding steps were timing recovery, phase correction and equalization. EVM as a function of SOA output power was analysed for different implementations (WNCT, LUT, FLUT, WNCT/LUT, WNCT/FLUT) in comparison with conventional SOA-based CO-OFDM system not using PAPR reduction nor linearization (Fig. 6). First, we can notice the effectiveness of the WNCT scheme once the SOA enters the saturated regime (P out -1.5 dBm); we get an improvement of around 2.5 dB at an EVM of 20%, which confirms the simulation results reported in [START_REF] Azou | Performance Improvement of a SOA-based Coherent Optical-OFDM Transmission System via Nonlinear Companding Transforms[END_REF]. Second, we observe a significant decrease in EVM by using a linearization scheme, FLUT being slightly better than LUT; however, at high power, convergence difficulties can degrade the performance. Then, it appears that combining WNCT and linearization leads to interesting results, a lower EVM being obtained over a wide range of P out ; in particular, we can notice that the FLUT gives the best results in this case. The gain in performance is also evident in the constellations (Fig. 6.b&c). We also illustrate in Fig. 7 the AM-AM and AM-PM characteristics for conventional system (at P out = 0.58 dBm) and with WNCT/FLUT implementation (at P out = 0.91 dBm); a huge improvement can be observed, with less scattered points (reduced memory effects) and better linearity. V. CONCLUSIONS Various predistortion schemes have been investigated for nonlinear effects mitigation in CO-OFDM systems with a special focus on the transmitter part, based on a SOA used as a booster. PAPR reduction is performed via WNCT and a FLUT scheme is considered for linearizing the transmitter. The influence of some key parameters of the FLUT block has been studied by considering a physical model of the SOA. Effectiveness of the WNCT is experimentally proven for the first time and we demonstrate that a joint use of WNCT/FLUT approaches leads to a significant performance improvement in terms of power efficiency and EVM, which can result in a reach extension of the optical link.

Figure 1 .

 1 Figure 1. CO-OFDM system with Matlab/ADS co-simulation, including baseband digital predistortion and nonlinear companding blocks (S/P: serial to parallel; IFFT: inverse Fast Fourier Transform; CP: Cyclic Prefix; D/A: Digital to Analog conversion; AWG: Additive White Gaussian).

Figure 2 .

 2 Figure 2. FLUT predistorter (z(n) stands for the baseband OFDM signal)

Figure 3 .

 3 Figure 3. EVM against various parameters of the FLUT (a) look-up table size (b) codebook size (c) FIR length.

Figure 4 .

 4 Figure 4. EVM vs. OFDM signal bandwidth for LUT & FLUT linearizations (left) and OSNR vs. bandwidth (right).

Figure 5 .

 5 Figure 5. Experimental setup. PC: Polarization controller; Att: Attenuator; IQ Mod: IQ modulator; AWG: Arbitrary waveform generator; ISO: Isolator;

Figure 6 .

 6 Figure 6. Experimental peak controlled linearization; (a) EVM vs. output power (Pout) (b) Constellation for conventional system at Pout = 0.58 dBm (c) Constellation for WNCT/FLUT-based system at Pout = 0.91 dBm.

Figure 7 .

 7 Figure 7. Experimental AM-AM (left) and AM-PM (right) characteristics without and with peak controlled linearization (via WNCT/FLUT).


Héctor Ríos

Laurentiu Hetel

Denis Efimov

Vector Lyapunov Function based Stability for a Class of Impulsive Systems

Keywords:

This paper contributes to the exponential stability analysis for impulsive dynamical systems based on a vector Lyapunov function and its divergence operator. The method relies on a 2D time domain representation. The results are applied to analyze the exponential stability of linear impulsive systems based on LMIs. Some examples illustrate the feasibility of the proposed approach.

I. INTRODUCTION

I MPULSIVE systems [1] represent an important class of hybrid systems [START_REF] Goebel | Hybrid dynamical systems[END_REF] containing discontinuities or jumps in the trajectories of the system governed by discrete dynamics. There exist a large variety of phenomena that are characterized by abrupt changes in the system state at certain instants, e.g. power electronics, sample-data systems, models in economics, bursting rhythm models in medicine, etc. According to the manner impulses are triggered several types of impulsive systems can be distinguishable: systems with time dependent impulses, systems with state dependent impulses (reset systems), and the combination of both of them (see, e.g. [START_REF] Bainov | Systems with Impulse Effects: Stability, Theory and Applications[END_REF], [START_REF] Goebel | Hybrid dynamical systems[END_REF] and [START_REF] Zavalishchin | Dynamic Impulse Systems. Theory and Applications[END_REF]). In the current work impulsive systems with time triggered impulses are studied.

In the context of stability analysis of impulsive systems an important effort has been made for linear dynamics. In [START_REF] Briat | A looped-functional approach for robust stability of linear impulsive systems[END_REF], a functional-based approach is developed for stability analysis. This method introduces looped functional leading to LMI conditions which allows to establish dwell-time results. Exponential stability, based on Lyapunov functions with discontinuity at the impulses, is proposed by [START_REF] Naghshtabrizi | Exponential stability of impulsive systems with application to uncertain sampled-data systems[END_REF] for nonlinear time-varying impulsive systems but the stability conditions are only applied to the linear case. In the sampledata system framework, there exist many works where the impulsive systems are used to describe the behavior of aperiodic sample-data systems (see, e.g. [START_REF] Hu | Constrained robust sampled-data control of nonlinear uncertain systems[END_REF] and [START_REF] Sivashankar | Characterization of the L 2induced norm for linear systems with jumps with applications to sampled-data systems[END_REF]). These approaches rely on the equivalent correspondence between continuous and discrete time domains showed in [START_REF] Seuret | A novel stability analysis of linear systems under asynchronous samplings[END_REF]. It is worth mentioning that there exist some works related to the stability for nonlinear 2D systems based on vector Lyapunov † Non-A team @ Inria, Parc Scientifique de la Haute Borne, 40 avenue Halley, 59650 Villeneuve d'Ascq, France; and CRIStAL (UMR-CNRS 9189), Ecole Centrale de Lille, BP 48, Cité Scientifique, 59651 Villeneuved'Ascq, France. Emails: hector.rios_barajas@inria.fr; denis.efimov@inria.fr ‡ CRIStAL (UMR-CNRS 9189), Ecole Centrale de Lille, BP 48, Cité Scientifique, 59651 Villeneuve-d'Ascq, France. Email: laurentiu.hetel@ec-lille.fr function approach, i.e. [START_REF] Emelianov | Stability and stabilization of differential nonlinear repetitive processes with applications[END_REF] and [START_REF] Emelianova | Vector lyapunov function based stability of a class of applications relevant 2d nonlinear systems[END_REF]. However, those results are not directly applicable for impulsive systems.

In this paper, based on a 2D time domain equivalence (see, e.g. [START_REF] Rogers | Control Systems Theory and Applications for Linear Repetitive Processes[END_REF] and [START_REF] Yeganefar | Lyapunov theory for 2-d nonlinear roesser models: Application to asymptotic and exponential stability[END_REF]), the exponential stability notion for a class of impulsive systems is studied using a vector Lyapunov function approach [START_REF] Lakshmikantham | Vector Lyapunov Functions and Stability Analysis of Nonlinear Systems[END_REF]. The proposed method contributes to the development of a stability analysis method for impulsive dynamical systems based on a divergence operator of the vector Lyapunov function in a 2D time domain. The main result is applied to analyze the exponential stability of linear impulsive systems based on LMIs. Some numerical examples illustrate the feasibility of the proposed approach.

The outline of this work is as follows. The problem statement and a preliminary definition of exponential stability are given in Section II. The main result as well as the conditions of exponential stability are described in Section III. Application of the developed theory to the problem of exponential stability for linear impulsive systems is considered in Section IV. Finally, some concluding remarks are discussed in Section V.

II. PROBLEM STATEMENT

Consider the class of impulsive dynamical systems for which the impulse times are external to the system and time dependent, i.e.

ẋ(t) = f 1 (x(t)) , ∀t ∈ R + \ I, (1)
x(t) = f 2 x(t -) , ∀t ∈ I, (2) x(0) = x 0 ,
where x, x 0 ∈ R n are the state of the system and the initial condition, respectively. The set of impulse times I := {t i } i∈N is a countable subset of R + with t i+1t i > 0, for all i ∈ N and t 0 = 0 and lim i→∞ t i = +∞ in order to avoid any zeno phenomena. The state trajectory is assumed to be right continuous and to have left limits at all times. The notation x(t - i+1 ) denotes the left limit of x(t) as t goes to t i+1 from the left, i.e. x(t -) = lim t↑ti+1 x(t). The distance between the impulses, i.e. the dwell-time (see, e.g. [START_REF] Liberzon | Switching in Systems and Control. Systems and Control: Foundations and Applications[END_REF] for more details), is defined as T i := t i+1t i ≤ T max , for all i ∈ N. The Lipschitz nonlinear functions f 1 and f 2 are such that f 1 (0) = 0 and f 2 (0) = 0. It is assumed that for all x 0 ∈ R n , f 1 is such that system (1) has a unique solution over any time interval [t i , t i+1 ) for all i ∈ N.

The proposed Vector Lyapunov function based approach relies on the embedding of system (1)-( 2) into a 2D time domain. Indeed, the entire state trajectory x(t) can be viewed as a sequence of the diagonal dynamics 1 of the following 2D system:

dx dt (t, k) = f 1 (x(t, k)) , ∀t ∈ [t i , t i+1 ) , (3)
x(t, k + 1) = f 2 (x(t, k)) , ∀t = t i+1 , ( 4
) ∀i = k ∈ N, t 0 = 0, x(0, 0) = x 0 0 , where x(t, k) : [t i , t i+1 )×N → R n is the current state vector, x(t i+1 , k + 1) : t i+1 × N → R n represents the reset vector state, x(t i+1 , k) : t i+1 × N → R n denotes the value of x just before the switching k+1, for all i = k ∈ N, i.e. x(t i+1 , k) = lim t↑ti+1 x(t, k), for t < t i+1 .
It is worth saying that the solutions of system (3) are unique for the diagonal dynamics, i.e. for all i = k ∈ N, t ∈ R + ; and they are represented by a sequence of functions {x(

t i + τ, k), τ ∈ (0, t i + 1]} i , for all i = k ∈ N.
Assume that |q| denotes the Euclidean norm of a vector q. The following stability definitions are introduced: Definition 1. A 2D system described by ( 3)-( 4), is said to be diagonal exponentially stable (DES) if there exist positive constants κ 1 , κ 2 , c and ε such that

|x(t i+1 , k + 1)| 2 ≤ cκ k+1 1 x 0 0 2 , 0 < κ 1 < 1 -ε, (5) |x(t, k)| 2 ≤ κ 2 |x(t i , k)| 2 , ∀t ∈ [t i , t i+1 ) , (6)
for all i, k ∈ N, and a small positive ε.

The goal in this work is to find the conditions for the stability of the impulsive systems described by ( 1)-( 2) by means of the 2D time representation (3)-(4).

III. MAIN RESULT: STABILITY OF IMPULSIVE SYSTEMS

In order to give the stability conditions a vector Lyapunov approach is used, i.e.

V (t, x(t, k), x(ti+1, k + 1)) = V1 (t, x(t, k)) V2 (x(ti+1, k + 1)) , (7)
where V 1 (t, x) > 0, V 2 (x) > 0, for all t ≥ 0 and x = 0, and V 1 (t, 0) = 0, V 2 (0) = 0, ∀t ≥ 0. Now, the divergence operator of the candidate function V along the trajectories of system (3)-( 4) is introduced for all t ∈ [t i , t i+1 ): [START_REF] Liberzon | Switching in Systems and Control. Systems and Control: Foundations and Applications[END_REF] i.e. V 1 is differentiable with respect to continuous time t while the difference in V 2 is calculated in discrete time. Based on the previous explanations, the following theorem is established: Theorem 1. Assume that for some T i there exist positive constants ε, c 1 , c 2 , c 3 , c 4 and c 5 such that the vector 1 The diagonal dynamics make reference only to those dynamics given by ( 3)-( 4) corresponding to i = k, for all i, k ∈ N and for all t ∈ R + .

divV (t, x(t, k), x(t i+1 , k + 1)) = dV 1 (t, x(t, k)) dt + V 2 (x(t i+1 , k + 1)) -V 2 (x(t i+1 , k)) ,
Lyapunov function V (t, x(t, k), x(t i+1 , k + 1)) and its divergence along the trajectories of the system (3)-( 4) satisfy, for all t ∈ [t i , t i+1 ), the following inequalities:

c 1 |x(t, k)| 2 ≤ V 1 (t, x(t, k)) ≤ c 2 |x(t, k)| 2 , ( 9
)
c 3 |x(t i+1 , k)| 2 ≤ V 2 (x(t i+1 , k)) ≤ c 4 |x(t i+1 , k)| 2 , ( 10
) divV ≤ -c 5 |x(t, k)| 2 + |x(t i+1 , k)| 2 , ( 11
)
c 2 (c 4 -c 5 ) ≤ c 1 c 5 ∨ T i < c 2 c 5 α, ( 12
)
c 2 c 5 γ < T i , ( 13
)
where γ = -ln c3(1-ε) c2 and α = -ln c2(c4-c5)-c1c5 c2(c4-c5) for all c 2 (c 4 -c 5 ) > c 1 c 5 .
Then, the 2D system ( 3)-( 4) is DES.

Proof. It is omitted due to lack of space.

The result given by Theorem 1 can be perfectly depicted by the following algorithm: Based on the statements given by Theorem 1, the Algorithm 1 provides some notions of minimum and maximum dwell-time depending on the structure of the system dynamics. In particular, the first and third case for diagonal exponential stability (pseudo-code lines: 5 and 13, Algorithm 1) give conditions for minimum dwell-time while the second case (pseudo-code lines: 7, Algorithm 1) provides conditions for maximum dwell-time.

In the following, the statements given by Theorem 1 are applied to prove the exponentially stability of linear impulsive systems.

IV. EXPONENTIAL STABILITY OF LINEAR IMPULSIVE SYSTEMS Consider the following class of linear impulsive systems

dx dt (t, k) = Ax(t, k), ∀t ∈ [t i , t i+1 ) , ( 14
)
x(t, k + 1) = Ex(t, k), ∀t = t i+1 , ( 15
) ∀i = k ∈ N, t 0 = 0, x(0, 0) = x 0 0 ,
where A and E are constant matrices with corresponding dimensions. It is assumed that A is Hurwitz and E is Schur. Then, the direct application of Theorem 1, provide us the following corollary that establishes the conditions to ensure the diagonal exponentially stability, described by Definition 1, for the linear impulsive system ( 14)-(15).

Corollary 1. Consider the vector Lyapunov function V (x(t, k), x(t i+1 , k + 1)) in [START_REF] Lakshmikantham | Vector Lyapunov Functions and Stability Analysis of Nonlinear Systems[END_REF], with V 1 (x(t, k)) = x T (t, k)P 1 x(t, k) and V 2 (x(t i+1 , k + 1)) = x T (t i+1 , k + 1)P 2 x(t i+1 , k + 1). If there exist matrices P 1 = P T 1 > 0 and P 2 = P T 2 > 0, such that

P 1 A + A T P 1 0 0 E T P 2 E -P 2 = -Q, ( 16
)
for Q = Q T > 0, and the constraints (12) and (13) hold; for c 1 = λ min (P 1 ), c 2 = λ max (P 1 ), c 3 = λ min (P 2 ), c 4 = λ max (P 2 ) and c 5 = λ min (Q); then, the linear impulsive system ( 14)-( 15) is DES. 2 Proof. It is clear that the vector Lyapunov function V (x(t, k), x(t i+1 , k + 1)) satisfies inequalities ( 9) and ( 10), i.e.

λmin(P1) |x(t, k)| 2 ≤ V1 (x(t, k)) ≤ λmax(P1) |x(t, k)| 2 , λmin(P2) |x(ti+1, k)| 2 ≤ V2 (x(ti+1, k)) ≤ λmax(P2) |x(ti+1, k)| 2 .
From the divergence definition and the previous inequalities, it is obtained that for all t ∈ [t i , t i+1 ) divV = x T (t, k) P 1 A + A T P 1 x(t, k)

+ x T (t i+1 , k + 1)P 2 x(t i+1 , k + 1) -x T (t i+1 , k)P 2 x(t i+1 , k),
and using the reset equation (15), it follows that

divV = x T (t, k) P 1 A + A T P 1 x(t, k) + x T (t i+1 , k) E T P 2 E -P 2 x(t i+1 , k).
From the matrix equality (16), it is obtained

divV ≤ -λ min (Q) |x(t, k)| 2 + |x(t i+1 , k)| 2 .
Then, the divergence of V satisfies the inequality [START_REF] Seuret | A novel stability analysis of linear systems under asynchronous samplings[END_REF], and according to the statements given by Theorem 1, system ( 14)-( 15) is DES, if the constraints ( 12) and ( 13) hold.

Example 1. Let us consider system ( 14)-( 15) with

A = -1 0 1 -2 , E = 0.5 0 0 0.5 .
Note that both matrices are stable, i.e. A is Hurwitz and E is Schur. Then, Corollary 1 is applied and the following feasible results are obtained:

P 1 =
5.2500 0.7500 0.7500 2.2500

, P 2 = 12I 2 , Q = 9I 4 , 2
The values of λ min (Y ) and λmax(Y ) correspond to the minimum and maximum eigenvalues of a matrix Y , respectively.

with c 1 = 2.0729, c 2 = 5.4271, c 3 = 12, c 4 = 12 and c 5 = 9. Therefore, according to the Algorithm 1, as it was expected, the impulsive system is DES for all T i > 0 (It is easy to check that c 4 > c 5 and c 2 (c 4c 5 ) ≤ c 1 c 5 hold and then the first case for diagonal exponential stability is obtained). The trajectories of the system, for T i = 0.5, are depicted in Fig. 1. In order to deal with linear impulsive systems where matrix A, or matrix E, is not Hurwitz, or not Schur, respectively; it is necessary to improve the structure of the Lyapunov function V (t, x(t, k), x(t i+1 , k + 1)), to be precisely V 1 (t, x(t, k)). In this sense, the following corollaries provide different ways to design V 1 (t, x(t, k)) and the conditions to get diagonal exponential stability.

Corollary

2. Consider the vector Lyapunov function V (t, x(t, k), x(t i+1 , k + 1)) in [START_REF] Lakshmikantham | Vector Lyapunov Functions and Stability Analysis of Nonlinear Systems[END_REF], with V 2 (x(t i+1 , k + 1)) = x T (t i+1 , k + 1)P 2 x(t i+1 , k + 1), and

V 1 (t, x(t, k)) = x T (t, k)P (t)x(t, k), ( 17
)
P (t) = P 1 + Δ(t) R -e A T Δ(t) E T REe AΔ(t) , ( 18
)
where

Δ(t) = T i -τ (t), with τ (t) = 1, and τ (t i ) = 0, ∀t ∈ [t i , t i+1 ). Then, if there exist matrices R = R T > 0, P 1 = P T 1 > 0 and P 2 = P T 2 > 0, such that Υ 0 + T i Υ 1 = -Q 1 , ( 19
) Υ 0 = -Q 2 , ( 20
)
for

Q 1 = Q 1 T > 0 and Q 2 = Q 2 T > 0, with Υ 0 = P 1 A + A T P 1 -R R E T P 2 E -P 2 -E T RE , Υ 1 = RA + A T R -A T RE 0 ,
and the constraints (12) and (13) hold; for c 1 = λ min (P 1 ), c 2 = λ max (P 1 ) + T i λ max (R) 1 + ϕ 2 , with ϕ = sup s∈[0,Ti] λ max (Ee AΔ(s) ), c 3 = λ min (P 2 ), c 4 = λ max (P 2 ) and c 5 = min(λ min (Q 1 ), λ min (Q 2 )); then, the linear impulsive system ( 14)-( 15) is DES.

Proof. First, taking into account that x(t i+1 , k + 1) = Ee AΔ(t) x(t, k), V 1 (t, x(t, k)) can be rewritten as follows

V 1 (t, x(t, k)) = x T (t, k)P 1 x(t, k) + Δ(t)ξ T (t, k)Rξ(t, x), ( 21
)
where

ξ(t, x) = x(t, k) -x(t i+1 , k + 1). Now, from the divergence definition it is obtained that for all t ∈ [t i , t i+1 ) divV = x T (t, k) P 1 A + A T P 1 x(t, k) + Δ(t)ξ T (t, k)R dξ(t, x) dt + Δ(t) dξ T (t, x) dt -ξ T (t, k) Rξ(t, x) + x T (t i+1 , k) E T P 2 E -P 2 x(t i+1 , k), divV = x T (t, k) P 1 A + A T P 1 x(t, k) + Δ(t)ξ T (t, k)RAx(t, k) + Δ(t)x T (t, k)A T -ξ T (t, k) Rξ(t, x) + x T (t i+1 , k) E T P 2 E -P 2 x(t i+1 , k).
Then, it follows that

divV ≤ x T (t, k) x T (t i+1 , k) × Υ 0 + (T i -τ (t)) Υ 1 × x(t, k) x(t i+1 , k) .
Since this LMI is affine in τ , its negative definiteness over the interval [t i , t i+1 ) is given by the negative definiteness at the vertices of the set, i.e. over the finite set τ ∈ {0, T i }; hence, the LMIs (19) and (20) are provided. Then, if these LMIs are feasible, the divergence of V satisfies [START_REF] Seuret | A novel stability analysis of linear systems under asynchronous samplings[END_REF] with

c 5 = min(λ min (Q 1 ), λ min (Q 2 )).
Finally, in order to take advantage of Theorem 1, it is necessary to prove that V 1 (t, x(t, k)) satisfies inequality [START_REF] Naghshtabrizi | Exponential stability of impulsive systems with application to uncertain sampled-data systems[END_REF]. From definition of V 1 (t, x(t, k)), it follows that

V1(t, x(t, k)) ≤ λmax(P1) |x(t, k)| 2 + Tiλmax(R) |ξ(t, k)| 2 , ≤ λmax(P1) |x(t, k)| 2 + Tiλmax(R) |x(t, k)| 2 + |Ex(ti+1, k)| 2 , ≤ λmax(P1) |x(t, k)| 2 + Tiλmax(R) |x(t, k)| 2 + Ee AΔ(t) x(t, k) 2 , ≤ λmax(P1) + Tiλmax(R) 1 + ϕ 2 |x(t, k)| 2 ,
where ϕ = sup s∈[0,Ti] λ max (Ee AΔ(s) ), thus; the upper bound is obtained. For the lower bound, it holds that V 1 (t, x(t, k)) ≥ λ min (P 1 ) |x(t, k)| 2 . Therefore, the lower bound is also obtained, and according to the statements given by Theorem 1, system ( 14)-( 15) is DES, if the constraints ( 12) and ( 13) hold.

Example 2. Let us consider system ( 14)-( 15) with

A = 1 3 -1 2
, E = 0.5 0 0 0.5 .

Note that the continuous dynamics is unstable while the discrete one is stable, i.e. A is not Hurwitz and E is Schur. Corollary 1 is applied together with a bisection-like approach and the following feasible results are obtained: with c 1 = 0.0032, c 2 = 0.7350, c 3 = 2.9415, c 4 = 4.4852 and c 5 = 0.0052. Therefore, according to the Algorithm 1, the impulsive system is DES for all 0.2499 > T i > 0 (It is easy to check that c 4 > c 5 and c 2 (c 4c 5 ) > c 1 c 5 hold and then the second case for diagonal exponential stability is obtained). The trajectories of the system, for T i = 0.24, are depicted in Fig. 2.

A = -1 0 1 -2 , E = 2 1 1 3 .
For this example the continuous dynamics is stable while the discrete one is unstable, i.e. A is Hurwitz and E is anti-Schur. Corollary 1 is applied and the following feasible results are obtained: Therefore, according to the Algorithm 1, the impulsive system is DES for all T i > 1.7092 (it is easy to check that c 5 > c 4 holds and then the third case for diagonal exponential stability is obtained). The trajectories of the system, for T i = 2, are depicted in Fig. 3.

V. CONCLUSIONS

In this paper the exponential stability notion for a class of impulsive systems is studied using a vector Lyapunov function approach. The proposed method contributes to the development of a stability analysis for impulsive dynamical systems based on a vector Lyapunov function and its divergence operator in a 2D time domain. The main result is applied to analyze the exponential stability of linear impulsive systems based on LMIs, and some numerical examples illustrate the feasibility of the proposed approach. The analysis of exponential stability of nonlinear systems is in the scope of the future research.

Algorithm 1 3 : if c 4 > c 5 then 4 :c 2 (c 4 -c 5 ) ≤ c 1 c 5 and T i > c 2

 13442452 Exponential Stability 1: Define the Lyapunov functions V 1 and V 2 2: Calculate the constants c 1 , c 2 , c 3 , c 4 and c 5 if

Figure 1 .

 1 Figure 1. Example 1 (Stable-Stable): Trajectories of the linear impulsive system for T i = 0.5.

Figure 2 .Remark 2 .Example 3 .

 223 Figure 2. Example 2 (Unstable-Stable): Trajectories of the linear impulsive system for T i = 0.24.

with c 1

 1 = 0.8225, c 2 = 2.1564, c 3 = 0.2, c 4 = 0.2 and c 5 = 3.

Figure 3 .Remark 3 .

 33 Figure 3. Example 3 (Stable-Unstable): Trajectories of the linear impulsive system for T i = 0.2.

ACKNOWLEDGMENT

This work was supported in part by Conseil Region Nord-Pas de Calais (ESTIREZ), the Government of Russian Federation (Grant 074-U01), the Ministry of Education and Science of Russian Federation (Project 14.Z50.31.0031), and ANR ROCC-SYS (ANR-14-CE27-0008).


Matteo Guerra
email: matteo.guerra@inria.fr

Carlos Vázquez

Gang Zheng Denis Efimov

Gang Zheng

Leonid Freidovich

Wilfrid Perruquetti

ϵ-Invariant Output Stabilization: Homogeneous Approach and Dead Zone Compensation

published or not. The documents may come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

The problem of uniform stabilization of dynamical systems in the presence of uncertain bounded inputs has a rather long history [START_REF] Tian | From poncelet's invariance principle to active disturbance rejection[END_REF]. By uniformity in this context we understand invariance (exact or approximate) of the closedloop system with respect to disturbing inputs (disturbance rejection or cancellation are another names of that problem). Initiated by a French engineer Jean-Victor Poncelet [START_REF] Preminger | Some considerations relating to control systems employing the invariance principle[END_REF], these ideas received a large attention in Soviet Union following the theory developed by Georgy Vladimirovich Shipanov [START_REF] Schipanov | Theory and methods of designing automatic regulators[END_REF], which is called the theory of ε-invariance (it was supposed to provide invariance up to ε > 0 deviations caused by disturbances of a given class). Next, many different solutions for ε-invariant stabilization have been proposed: time delay control [START_REF] Youcef-Toumi | Input/output linearization using time delay control[END_REF], active disturbance rejection [START_REF] Han | A class of extended state observers for uncertain systems[END_REF], universal integral controls [START_REF] Khalil | Universal integral controllers for minimum-phase nonlinear systems[END_REF], [START_REF] Freidovich | Performance recovery of feedbacklinearization-based designs[END_REF], various sliding-mode control algorithms [START_REF] Levant | Universal siso sliding-mode controllers with finite-time convergence[END_REF], [START_REF] Ferreira | Robust control with exact uncertainties compensation: With or without chattering? Control Systems Technology[END_REF] converging in a finite time, model-free control [START_REF] Fliess | Commande sans modèle et commande à modèle restreint[END_REF] (just to mention a few, there are also many other adaptive/fuzzy/neural control solutions). The statement of ε-invariant control design problem can be given following a recent development [START_REF] Fliess | Commande sans modèle et commande à modèle restreint[END_REF] (model-free control). Consider a SISO uncertain nonlinear system, whose model is given in the implicit form (it is not resolved with respect to the highest derivative):

f [y(t), ẏ(t), . . . , y (n) (t), u(t), d(t)] = 0, t ≥ 0, where y(t) ∈ R is the measured output, u(t) ∈ R is the control input, d(t) ∈ R m is the vector of uncertain parameters/signals, n ≥ 1 is the system dimension, which may be unknown, f : R n+m+1 → R is an unknown nonlinear function ensuring existence of the system solutions at least locally. Fixing k ≥ 1, a local model can be extracted:

y (k) (t) = u(t) + F (t),
where F (t) ∈ R is a new unknown input including y, y (1) , . . . , y (n) , u and d. This model may have sense only locally, but under assumption that the dynamics of y (k+1) , . . . , y (n) are stable (i.e. the system is minimum phase with relative degree k [START_REF] Khalil | Universal integral controllers for minimum-phase nonlinear systems[END_REF], [START_REF] Fliess | Commande sans modèle et commande à modèle restreint[END_REF]) the original stabilization problem for uncertain nonlinear system can be reduced to uniform (ε-invariant) stabilization of a chain of k integrators subjected by unknown matched input F (frequently assumed bounded along with its derivatives). There are many solutions to this problem, which are based on the idea that if it is possible to estimate y (k) (t) then F (t) = y (k) (t) -u(t) can be evaluated and compensated by the control. The difference is mainly in the tools used for estimation of y (k) (t) (highgain observers in [START_REF] Khalil | Universal integral controllers for minimum-phase nonlinear systems[END_REF], [START_REF] Freidovich | Performance recovery of feedbacklinearization-based designs[END_REF], sliding-mode differentiators in [START_REF] Levant | Universal siso sliding-mode controllers with finite-time convergence[END_REF], [START_REF] Ferreira | Robust control with exact uncertainties compensation: With or without chattering? Control Systems Technology[END_REF] or algebraic ones in [START_REF] Fliess | Commande sans modèle et commande à modèle restreint[END_REF]). Time delay is frequently introduced to break the algebraic loop [START_REF] Youcef-Toumi | Input/output linearization using time delay control[END_REF], [START_REF] Back | Adding robustness to nominal output feedback controllers for uncertain nonlinear systems: a nonlinear version of disturbance observer[END_REF], [START_REF] Guo | Anti-disturbance control theory for systems with multiple disturbances: A survey[END_REF], [START_REF] Huang | Active disturbance rejection control: Methodology and theoretical analysis[END_REF], which appears when using the estimate y (k) (t) -u(t) in the control u(t) itself. Another difference between [START_REF] Ferreira | Robust control with exact uncertainties compensation: With or without chattering? Control Systems Technology[END_REF], [START_REF] Fliess | Commande sans modèle et commande à modèle restreint[END_REF], [START_REF] Freidovich | Performance recovery of feedbacklinearization-based designs[END_REF], [START_REF] Khalil | Universal integral controllers for minimum-phase nonlinear systems[END_REF], [START_REF] Levant | Universal siso sliding-mode controllers with finite-time convergence[END_REF] consists in the type of feedback used for the system stabilization. Theoretically sliding-mode controls provide a finite-time exact cancellation of matched disturbances [START_REF] Levant | Universal siso sliding-mode controllers with finite-time convergence[END_REF], [START_REF] Ferreira | Robust control with exact uncertainties compensation: With or without chattering? Control Systems Technology[END_REF], which is better than ε-invariance provided by linear feedbacks from [START_REF] Fliess | Commande sans modèle et commande à modèle restreint[END_REF], [START_REF] Freidovich | Performance recovery of feedbacklinearization-based designs[END_REF], [START_REF] Khalil | Universal integral controllers for minimum-phase nonlinear systems[END_REF], [START_REF] Youcef-Toumi | Input/output linearization using time delay control[END_REF]. But in practice the slidingmode controls suffer from chattering that returns them back to ε-invariance setting. A related difference is robustness with respect to different nonlinearities of y, y (1) , . . . , y (n) hidden in F (for example, linear feedback treats only Lipschitz or linear perturbations). In order to improve robustness and to avoid chattering, an intermediate solution could be represented by homogeneous high-gain controls [START_REF] Bhat | Geometric homogeneity with applications to finite-time stability[END_REF] and observers [START_REF] Perruquetti | Finite-time observers: application to secure communication[END_REF]. Due to homogeneity, local asymptotic stability of this systems implies global one, and robustness with respect to disturbances is inherited next [START_REF] Bernuau | Verification of ISS, iISS and IOSS properties applying weighted homogeneity[END_REF]. Adjusting nonlinear gains in control and estimation algorithms from [START_REF] Bhat | Geometric homogeneity with applications to finite-time stability[END_REF], [START_REF] Perruquetti | Finite-time observers: application to secure communication[END_REF] it is possible to get a needed degree of robustness with respect to F . A development of ε-invariant output control based on [START_REF] Bhat | Geometric homogeneity with applications to finite-time stability[END_REF], [START_REF] Perruquetti | Finite-time observers: application to secure communication[END_REF] is presented in this work and thus applied to control a specific link of an actuated industrial crane used in forestry. The aim is to achieve precision and smooth extension and retraction in a way that standard methods cannot guarantee [START_REF] Han | From pid to active disturbance rejection control[END_REF]. In the treated example the telescopic link of the crane must track a reference trajectory using the proposed approach to compensate incertitudes due to a dead zone of the control input modeled as in [START_REF] Tao | Adaptive control of plants with unknown dead-zones[END_REF] along with other perturbations.

II. PRELIMINARIES

Through the paper the following notations is used: • The notation DV (x)f (x) stands for the directional derivative of a continuously differentiable function V with respect to the vector field f evaluated at point x. Following [START_REF] Dashkovskiy | Input to state stability and allied system properties[END_REF], consider a nonlinear system

• R + = {x ∈ R : x ≥ 0},
ẋ(t) = f [x(t), d(t)], (1)
where x(t) ∈ R n is the state, d(t) ∈ R m is the external input, d ∈ L ∞ , and f : R n+m → R n is a locally Lipschitz (or Hölder) continuous function, f (0, 0) = 0. For an initial condition x 0 ∈ R n and input d ∈ L ∞ , define the corresponding solutions by x(t, x 0 , d) for any t ≥ 0 for which the solution exists.

Definition 1. The system (1) is called input-to-state practically stable (ISpS), if for any input d ∈ L ∞ and any x 0 ∈ R n there are some functions β ∈ KL, γ ∈ K and c ≥ 0 such that

x(t, x 0 , d) ≤ β( x 0 , t) + γ(||d|| [0,t) ) + c ∀t ≥ 0.
The function γ is called nonlinear asymptotic gain. The system is called

ISS if c = 0. Definition 2. A smooth function V : R n → R + is called ISpS Lyapunov function for the system (1) if for all x ∈ R n , d ∈ R m and some r ≥ 0, α 1 , α 2 , α 3 ∈ K ∞ and θ ∈ K: α 1 ( x ) ≤ V (x) ≤ α 2 ( x ), DV (x)f (x, d) ≤ r + θ( d ) -α 3 ( x ). Such a function V is called ISS Lyapunov function if r = 0.
Note that an ISS Lyapunov function can also satisfy the following equivalent condition for some χ ∈ K:

x > χ( d ) ⇒ DV (x)f (x, d) ≤ -α 3 ( x ).
Theorem 1. [START_REF] Dashkovskiy | Input to state stability and allied system properties[END_REF] The system (1) is ISS (ISpS) iff it admits an ISS (ISpS) Lyapunov function.

A. Weighted homogeneity

Following [START_REF] Bacciotti | Lyapunov Functions and Stability in Control Theory[END_REF], for fixed strictly positive numbers r i , i = 1, . . . , n called weights and λ > 0, one can define:

• the vector of weights r = (r 1 , . . . , r n ) T , r max = max 1≤j≤n r j and r min = min 1≤j≤n r j ;

• the dilation matrix function Λ r (λ) = diag{λ ri } n i=1 , note that ∀x ∈ R n and ∀λ > 0 we have Λ r (λ)x = (λ r1 x 1 , . . . , λ ri x i , . . . , λ rn x n ) T . Definition 3. A function g : R n → R is r-homogeneous with degree µ ∈ R if ∀x ∈ R n and ∀λ > 0 we have: λ -µ g(Λ r (λ)x) = g(x).
A vector field f : R n → R n is r-homogeneous with degree ν ∈ R, with ν ≥ -r min if ∀x ∈ R n and ∀λ > 0 we have:

λ -ν Λ -1 r (λ)f (Λ r (λ)x) = f (x),
which is equivalent for i-th component of f being a rhomogeneous function of degree r i + ν.

The system (1) with d = 0 is r-homogeneous of degree ν if the vector field f is r-homogeneous of degree ν.

Theorem 2.

[21] For the system (1) with d = 0 and r-homogeneous and continuous function f the following properties are equivalent:

• the system (1) is (locally) asymptotically stable;

• there exists a continuously differentiable rhomogeneous Lyapunov function V :

R n → R + such that α1( x ) ≤ V (x) ≤ α2( x ), DV (x)f (x, 0) ≤ -α( x ), λ -µ V (Λr(λ)x) = V (x), µ > rmax, ∀x ∈ R n and ∀λ > 0, for some α 1 , α 2 ∈ K ∞ and α ∈ K. Define f (x, d) = [f (x, d) T 0 m ] T ∈ R n+m ,
it is an extended auxiliary vector field for the system (1), where 0 m is the zero vector of dimension m.

Theorem 3.

[4] Let the vector field f be homogeneous with the weights r =

[r 1 , . . . , r n ] > 0, r = [r 1 , . . . , rm ] > 0 with a degree ν ≥ -r min , i.e. f (Λ r (λ)x, Λ r (λ)d) = λ ν Λ r (λ)f (x, d) for all x ∈ R n , d ∈ R m and all λ > 0.
Assume that the system (1) is globally asymptotically stable for d = 0, then the system (1) is ISS.

Therefore, for homogeneous system (1) its ISS property follows asymptotic stability for d = 0 (as for linear systems [START_REF] Dashkovskiy | Input to state stability and allied system properties[END_REF]). The nonlinear asymptotic gain function has been also estimated in [START_REF] Bernuau | Verification of ISS, iISS and IOSS properties applying weighted homogeneity[END_REF].

B. Homogeneous stabilizing control

Consider a nonlinear system

ξi = ξ i+1 , i = 1, . . . , n -1, (2)
ξn = - n i=1 a i ξ i αi ,
where ξ = [ξ 1 , . . . , ξ n ] ∈ R n is the state vector, α i and a i are real parameters. For r i = 1 + (i -1)ν and α i = 1+nν 1+(i-1)ν , i = 1, . . . , n, where ν > -1 n-1 , the system (2) is r-homogeneous of degree ν. Theorem 4. [START_REF] Bhat | Geometric homogeneity with applications to finite-time stability[END_REF] Let a 1 , . . . , a n form a Hurwitz polynomial, then there exists 0 < < 1 n-1 such that for any ν ∈ (-1 n-1 + , 0) the system (2) with α i = 1+nν 1+(i-1)ν , i = 1, . . . , n is globally finite-time stable.

C. Homogeneous observer

Consider a nonlinear system ξi = ξ i+1 -λ i ξ 1 βi , i = 1, . . . , n -1, (3) ξn
= -λ n ξ 1 βn ,
where ξ = [ξ 1 , . . . , ξ n ] ∈ R n is the state vector, λ i and β i are real parameters. For r i = 1 + (i -1)µ and β i = 1 + iµ for all i = 1, . . . , n, where µ > -1 n , the system (3) is rhomogeneous of degree µ. Theorem 5. [START_REF] Perruquetti | Finite-time observers: application to secure communication[END_REF] Let λ 1 , . . . , λ n form a Hurwitz polynomial, then there exists 0 < < 1 n such that for any µ ∈ (-1 n + , 0) the system (3) with β i = 1 + iµ, i = 1, . . . , n is globally finite-time stable.

III. PROBLEM STATEMENT

The following state-space representation will be considered in this work:

ẋi (t) = x i+1 (t), i = 1, . . . , k -1, (4) ẋk (t) = u(t) + F (t), y(t) = x 1 (t),
where x(t) = [x 1 (t), . . . , x k (t)] T ∈ R k is the state space vector of the system (4) at time instant t ≥ 0; u(t) ∈ R and F (t) ∈ R are the control and disturbance inputs, respectively; y(t) ∈ R is the measured output. Since (4) is linear, then the measurement noise can be transferred to the input and included in F . The following restrictions are introduced for (4).

Assumption 1. Let F ∈ L ∞ and Ḟ ∈ L ∞ , in addition a constant f > 0 is given such that ||F || ∞ ≤ f, || Ḟ || ∞ ≤ f.
It is required to design a dynamical output feedback u such that for the given ε > 0 and all initial conditions x 0 ∈ R k , lim t→+∞ x(t) ≤ ε for all F satisfying Assumption 1. The conditions of that assumption can also be relaxed supposing that F is a nonlinear function of x and asking for a semi-global εinvariance.

IV. CONTROL DESIGN First, the vector x has to be estimated. Due to the system structure this problem is equivalent to the estimation of the derivatives y (1) (t), . . . , y (k-1) (t) for the output y(t), for this purpose the following linear filter can be designed

żi = z i+1 + l i (y -z 1 ), i = 1, . . . , k -1, (5) żk = l k (y -z 1 ) + u,
where z = [z 1 , . . . , z k ] T ∈ R k and high-gain tuning parameters l i > 0 for i = 1, . . . , k form a Hurwitz polynomial (more precise restrictions on l i will be given later). Denote ŷ(i) as an estimate of y (i) = x i+1 , then we can select ŷ(i) = z i+1 for i = 0, . . . , k -1 and the filter estimation error e = x -z has dynamics:

ėi = e i+1 -l i e 1 , i = 1, . . . , k -1, ėk = -l k e 1 + F.
From the last equation the following estimate

F = êk + l k e 1
of F can be calculated, where êk is an estimate of ėk . In order to calculate êk a second filter/differentiator should be designed that has to converge faster than exponentially (the rate of decay in the linear one ( 5)). For this purpose a homogeneous high-gain differentiator can be used:

ζi = ζ i+1 -l i e 1 + λ i e 1 -ζ 1 βi , i = 1, . . . , k, (6) ζk+1
˙ i = i+1 -λ i 1 βi , i = 1, . . . , k, ˙ k+1 = -λ k+1 1 β k+1 -l k e 2 + Ḟ , which is r-homogeneous of order µ > -1 k+1 for r i = 1 + (i -1)µ and β i = 1 + iµ for all i = 1, . . . , k + 1 (if µ = -1
k+1 then β k+1 = 0 and (6) reduces to a high-order sliding mode observer, while for all µ > -1 k+1 the filter (6) stays continuous), i.e. a perturbed version of (3). Therefore,

F = ζ k+1 + l k e 1 . (7)
Second, the control u can be introduced

u = - k-1 i=0 a i+1 ŷ(i) αi+1 -F = - k i=1 a i z i αi -F , (8)
where the coefficients a i , i = 1, . . . , k form a Hurwitz polynomial and α i > 0, i = 1, . . . , k will be defined in the next section.

To conclude, the proposed model-free invariance control algorithm includes two filters ( 5) and ( 6) (one for differentiation and another for decoupling the control u and the estimate of F appearing into the same equation), the unknown input F estimate [START_REF] Ferreira | Robust control with exact uncertainties compensation: With or without chattering? Control Systems Technology[END_REF] and the stabilizing control [START_REF] Fliess | Commande sans modèle et commande à modèle restreint[END_REF].

Therefore, the following result can been proven 1 .

Theorem 6. Let Assumption 1 be satisfied. Then for any given > 0 there exist

l ∈ R k , λ ∈ R k+1 , a ∈ R k , µ ∈ (-1 k+1 , 0), ν ∈ (-1 k-1 , 0), β i = 1 + iµ for i = 1, . . . , k + 1, α i = 1+kν
1+(i-1)ν for i = 1, . . . , k, such that in the system (4) with the output regulator ( 5)-( 8) for all initial conditions

x ∈ L k ∞ , z ∈ L k ∞ , ζ ∈ L k+1 ∞ ,
and

lim t→+∞ x(t) ≤ ε.
Moreover, the system (4), ( 5)-( 8) is ISS with respect to the input (F, Ḟ ).

Remark 1. Since dynamics of all variables, x(t), (t) and e(t), are homogeneous, their asymptotic gains can be evaluated as it is proposed in [START_REF] Back | Adding robustness to nominal output feedback controllers for uncertain nonlinear systems: a nonlinear version of disturbance observer[END_REF] and using the parameters γ l , γ λ and γ a . Finally, for given f the value of ε(t) can be estimated.

Since ( 6) and ( 8) contain nonlinear gains, then the asymptotic gain of ( 4), ( 5)-( 8) (see [START_REF] Bernuau | Verification of ISS, iISS and IOSS properties applying weighted homogeneity[END_REF] for an algorithm of its estimation) close to the origin is better than in a pure linear system (i.e. replacing ( 6) and ( 8) by linear filter and feedback, respectively).

V. DEAD ZONE COMPENSATION

We are going to present the compensation of an input nonlinearity, known as a dead zone, which is depicted in Fig. 1. This dead zone model is a static representation of diverse physical phenomena with negligible fast dynamics, see [START_REF] Tao | Adaptive control of plants with unknown dead-zones[END_REF]. One well-known example is the model of an industrial electro-hydraulic valve in which the spool occludes the orifice with some overlap. In this case, system (4) should be rewritten as below:

ẋi (t) = x i+1 (t), i = 1, . . . , k -1, (9) ẋk (t) = D(u(t)) + F (t), y(t) = x 1 (t),
where the dead zone input is represented by D(u(t)) and it has the following structure:

D(u(t)) =    m r (u -b r ) if u ≥ b r , 0 if -b l ≤ u ≤ b r , m l (u -b l ) if u ≤ -b l . ( 10
)
where m i = m 0 + ∆m i and b i = b 0 + ∆b i , with i = l, r; the subscript l stands for "left" and r for "right", m 0 and b 0 are the nominal values while ∆m i and ∆b i are uncertain terms. 1 The Proof is omitted due to space limitation. Let u 0 (t) be the control signal from a model-free invariance control design. Then, our approach is based on the design of a nominal dead zone inverse, DI(u 0 (t)), where the remaining uncertain terms ∆m i and ∆b i will be canceled by the invariance control algorithm. For this aim the nominal parameters m 0 and b 0 are assumed to be known and they are used for the construction of a static nominal dead zone inverse: 11) Substituting ( 11) in ( 9) we obtain the following dynamics:

u(t) = DI(u 0 (t)) = m -1 0 (u 0 (t) + m 0 b 0 sign(u 0 )). (
ẋi (t) = x i+1 (t), i = 1, . . . , k -1, (12) ẋk
(t) = u 0 (t) + F 0 (t), y(t) = x 1 (t),
where u 0 (t) is the final control input and F 0 (t) represents a perturbation containing new terms related with the uncertain parameters of a dead zone. Note that the inverse of a dead zone is a relay-type discontinuity that can be canceled if the inverse is exact, see [START_REF] Tao | Adaptive control of plants with unknown dead-zones[END_REF]. Besides, with a nominal dead zone inverse the structure of system ( 4) is recovered and the uncertain terms are to be compensated by a model free invariance controller.

VI. HYDRAULIC ACTUATOR CASE STUDY

The experimental setup is the telescopic link of a laboratory prototype of a typical industrial hydraulic forestry crane. Such industrial equipment is widely used in forestry and is a subject of many researches aimed at automation of these systems, see [START_REF] Papadopoulos | On modeling, identification, and control of a heavy-duty electrohydraulic harvester manipulator[END_REF]. The telescopic link of the crane, see Fig. 2, consists of a double-acting single-side hydraulic cylinder and a solid load, which is attached to a piston of the cylinder. This link can be described as a restricted 1-DOF mechanical system actuated by a hydraulic force:

mẍ = f h -f grav -f f ric , (13)
where m is the mass, f h is the force generated by the hydraulics, f grav is the gravity and f f ric is the friction force. The force generated by the hydraulics is controlled via a current signal to the valve of the cylinder and it is given by

f h = p a A a -p b A b , if variables are in R, f uncertain , otherwise, (14)
where the validity region

R = {x ∈ (0, xp ), | ẋ| < xv , p a , p b ∈ (p t , p s )} , (15)
is defined by certain constant bounds; p s is the pump pressure, and p t is the return (exit) pressure; the piston areas A a and A b are known geometric parameters; p a and p b are the measured pressures in chambers A and B of the cylinder. The friction and gravity terms are considered as unknown perturbations. Note that the position of the link, x, is limited by geometrical constrains, the velocity is limited by the maximum achievable flow from a pump, pressures are limited through a set of service anti-cavitation and pressurerelief valves, which, in particular, ensure p t ≤ p i ≤ p s , i = a, b. These devices play a fault-preventing role and do not influence a normal operation. Moreover, the initial conditions are within the region R and f uncertain prevent from leaving the region. The dynamics of the pressures can be modeled, see [START_REF] Merrit | Hydraulic Control Systems[END_REF]Sec. 3.8], by:

ṗa = β Va(x) (-ẋ A a + q a ) , if variables are in R, ṗb = β V b (x) ( ẋ A b -q b ) , if variables are in R, ṗi = p uncertain , i = a, b otherwise, (16)
and initial conditions p a,b (0) ∈ (p t , p s ); V a (x) = V a0 +x A a and V b (x) = V b0 -x A b are volumes of the chambers A and B at the given piston position x, V a0 and V b0 are known geometric constants, β is a known bulk modulus, q a and q b are flows to the chamber A and from the chamber B. The flow q a is positive when the oil goes into chamber A, and the flow q b is positive when the oil goes out of chamber B. Following [START_REF] Merrit | Hydraulic Control Systems[END_REF], [START_REF] Niksefat | Designing robust force control of hydraulic actuators despite system and enviromental uncertainties[END_REF], the nonlinear equations describing the fluid flow distribution in the valve can be written, in their simplest forms, as:

q a = c a S a (x s ) √ p s -p a and q b = c b S b (x s ) √ p b -p t for x s ≥ 0; q a = -c a S a (x s ) √ p a -p t and q b = -c b S b (x s ) √ p s -p b for x s < 0.
Here c a and c b are constant coefficients which depend on physical values (fluid density, discharge coefficient and other), S a (x s ) and S b (x s ) are (non-negative) areas of orifices for the ports A and B, x s is a displacement of a spool inside a valve, this spool is actuated by an electromagnetic actuator where an input (current) signal u is applied, u ∈ [-u -, u + ]. Assuming that the valve is symmetric, i.e. ∀x s : S a (x s ) = S b (x s ), we introduce the signed area function S(x s ) given by S(x s ) = S a (x s ) sign(x s ) = S b (x s ) sign(x s ). The absolute value of this function is equal to the area of the orifices and the sign indicates directions of the flows. Then, the flow equations can be rewritten as: q a = c a ϕ a S(x s ) and q b = c b ϕ b S(x s ), with: . Note that in industrial hydraulic systems a nonzero pressure difference through the valve is ensured by a set of service valves, i.e. ϕ a ≥ 0 and ϕ b ≥ 0. Taking into account a high-response servo valve and assuming that the spool displacement is proportional to the input signal, the signed area can be modeled as defined by the input signal u through a nonlinear static relation S(x s ) = D(u). The shape of the function D strongly depends on a type of the valve; for industrial heavyduty systems a dead-zone due to a leakage-preventing closedcenter spool and a saturation due to limiting screws are common, see [START_REF] Aranovskiy | Modeling and identification of spool dynamics in an industrial electro-hydraulic valve[END_REF] and references therein. Taking the derivative of ( 14) and substituting [START_REF] Merrit | Hydraulic Control Systems[END_REF], one obtains:

ϕ a = √ p s -
ḟh = -ϕ 0 ẋ + ϕ 1 D(u), (17)
where

ϕ 0 = β( A 2 a Va(x) + A 2 b V b (x) ) and ϕ 1 = β( caAaϕa Va(x) + c b A b ϕ b V b (x) ) with 0 < ϕ i ≤ ϕ i ≤ ϕ i for i = 0, 1.
Besides, the nonlinear function D(u) can be represented by [START_REF] Guo | Anti-disturbance control theory for systems with multiple disturbances: A survey[END_REF]. Since ḟh is bounded, sytem (17) can be rewritten as:

ẋ = ϕD(u) -ϕ -1 0 ḟh , (18)
where ϕ = ϕ1 ϕ0 and ϕ = φ + ∆ φ. Given an appropriate desired trajectory for (15), x ref , with Lipschitz continuous second derivative, ẍref , the objective of this section is to design a control law for u to achieve the tracking of the cylinder position x.

A. Experiments

The experiments are carried out using a real-time platform dSpace 1401 with sample time of 1ms using forward Euler integration method. The position of the telescopic link is measured with a wire-actuated encoder which provides 2381 counts for the range from 0 to 1.55m with a quantization interval of Q = 0.651mm. The desired trajectory is selected as a sinusoidal signal: x d = 0.8 + 0.4 sin(ωt). In the experiments the values used for the dead zone compensation [START_REF] Han | A class of extended state observers for uncertain systems[END_REF] are m 0 = 1 and b 0 = 0.3; being such values the nominal ones, the controller should be able to compensate the gap between nominal and real values, namely compensate the remaining uncertain terms ∆m i and ∆b i along with other perturbations. Firstly, it is necessary to tune the a 1 , λ 1 and l coefficients, then to examine the performances of the controller varying the parameter α 1 which determines the non-linearity of the controller, lastly on has to evaluate the effect of the F0 estimate. The error in the tracking is smaller decreasing α 1 and even smaller adding the F0 compensation, which proves clearly the efficiency of the method. It is worth to remark thought, that a high frequency oscillatory behavior is present in the error (and control input) for α 1 < 1 which is not desirable in the application considered. An α 1 equal or close to one had to be chosen to avoid it. The closest α 1 to 1 the faster the oscillations, which are not propagated from the control input to the telescopic link guaranteeing a Comparison between the Homogeneous controller and the PI approach. (x-axis in samples, y-axis in meter) smooth extension and retraction. To lay stress on the good performances of the method, Fig. 3 shows the behavior in comparison with a standard PI approach often used to regulate such a system. It can be seen that, although the integral action applied, the link could not follow the reference as well as the homogeneous control which does not have what appears to be a static error with respect to the reference; that it could be caused by the nominal values used for the dead zone compensation. As stated above the homogeneous controller is able to handle the error between nominal and real values, that's not achieved with the PI control despite the tuning to have the best behavior possible.

VII. CONCLUSION

This paper presented the synthesis of a control method for system which can be described as in [START_REF] Bernuau | Verification of ISS, iISS and IOSS properties applying weighted homogeneity[END_REF]. A linear filter allows the decoupling of the control variables and the disturbances, which are afterward estimated starting from the estimation error using a homogeneous high gain observer. The estimation of perturbations is then included in the controller to compensate them. The approach is fully proven theoretically then applied on the telescopic link of a hydraulic actuated robotic crane used in forestry. To deal with such a platform an easy but effective dead zone compensation technique has been implemented; such a compensation relies just on nominal parameters. The proposed controller behaves as expected, compensating perturbations and dead zone incertitudes, and outperforming the standard PI control as it is shown with experiments.

= λ k+1 e 1 -ζ 1 β k+1 + l k l 1 ,

 1 where ζ = [ζ 1 , . . . , ζ k ] T ∈ R k+1 and the tuning parameters β i > 0 and λ i > 0 for i = 1, . . . , k + 1 will be derived later. Denote ē = [e T e k+1 ] T , where e k+1 = -l k e 1 + F, ėk+1 = -l k ė1 + Ḟ = -l k (e 2 -l 1 e 1 ) + Ḟ , then we can select êk = ζ k+1 and the estimation error = ē -ζ has dynamics:

Fig. 1 .

 1 Fig. 1. Dead zone input nonlinearity.

Fig. 2 .

 2 Fig. 2. Industrial hydraulic forestry crane.

4 -

 4 Fig. 3.Comparison between the Homogeneous controller and the PI approach. (x-axis in samples, y-axis in meter)

  A continuous function α : R + → R + belongs to the class K if α(0) = 0 and the function is strictly increasing. The function α : R + → R + belongs to the class K ∞ if α ∈ K and it is unbounded. A continuous function β : R + × R + → R + belongs to the class KL if β(•, t) ∈ K ∞ for each fixed t ∈ R + and lim t→+∞ β(s, t) = 0 for each fixed s ∈ R + .

	where R is the set of real
	number.
	• | • | denotes the absolute value in R, • denotes the
	Euclidean norm on R n .
	• For a (Lebesgue) measurable function d : R + → R m
	define the norm ||d|| [t0,t1) = ess sup t∈[t0,t1) d(t) ,
	then ||d|| ∞ = ||d|| [0,+∞) and the set of d(t) with the
	property ||d|| ∞ < +∞ we further denote as L ∞ (the
	set of essentially bounded measurable functions).

• • • α denotes the following operation | • | α sign(•).

ACKNOWLEGMENT

This work was supported in part by Région Nord-Pas de Calais, by the Government of Russian Federation (Grant 074-U01) and the Ministry of Education and Science of Russian Federation (Project 14.Z50.31.0031).