Mokrane Boudaoud
email: mokrane.boudaoud@isir.upmc.fr

Marcelo Gaudenzi De Faria

Yann Gorrec
email: legorrec@femto-st.fr

Yassine Haddab
email: yassine.haddab@femto-st.fr

Philippe Lutz
email: philippe.lutz@femto-st.fr

An output feedback LPV

Keywords: Microelectromechanical systems, Microrobotics, Linear Parameter Varying (LPV), Singular implicit modeling, H∞ control, Linear Matrix Inequality (LMI)

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

An Output feedback LPV control strategy of a nonlinear electrostatic microgripper through a singular implicit modeling Mokrane Boudaoud * , * * Marcelo Gaudenzi De Faria * Yann Le Gorrec * Yassine Haddab * Philippe Lutz *

INTRODUCTION

The manipulation of micrometer sized and nanometer sized objects is very challenging in microrobotics [START_REF] Chaillet | Microrobotics for Micromanipulation[END_REF]). For instance, single living cells manipulation has a great interest for a range of biomedical disciplines [START_REF] Sun | Biological cell injection using an autonomous microrobotic system[END_REF]; Kimuraet al. [1995]). The size of biological cells is ranging from few nanometers to several tens of micrometers. Microgrippers are therefore an interesting solution to perform manipulation tasks at such scales. MEMS (Micro Electro Mechanical systems) based microgrippers using electrostatic comb drive actuators are one of the most popular micromanipulation systems in biological research [START_REF] Muntwyler | Monolithically Integrated Two-Axis Microtensile Tester for the Mechanical Characterization of Microscopic Samples[END_REF]; [START_REF] Yamahata | Silicon Nanotweezers with Adjustable and Controllable Gap for the Manipulation and Characterization of DNA Molecules[END_REF]). Although, comb drive actuators are highly nonlinear when exceeding few micrometer displacements [START_REF] Legtenberg | Comb drive actuators for large displacements[END_REF]), they allow positioning gripping arms over several tens of micrometers, they have no hysteretic behavior and they do not generate heating.

Nonlinearities of comb drive actuators are due to electrostatic potentials [START_REF] Lee | Nonlinear behaviors of a comb drive actuator under electrically induced tensile and compressive stresses[END_REF]) and axial forces acting on compliant mechanisms [START_REF] Legtenberg | Comb drive actuators for large displacements[END_REF]). The side instability is also a non-linearity that occurs when the gap spacing between fixed and movable electrodes of the actuator is no longer constant [START_REF] Alwi | Sliding Mode Control Design of an Electrostatic Microactuator using LPV schemes[END_REF]). It is therefore very difficult to control such actuators with robust performances in a wide operating range. This is nevertheless fundamental for an efficient manipulation of fragile samples such as biological cells regardless their size.

An extensive literature is available for the control of comb drive actuators [START_REF] Bryzek | Control Issues for MEMS[END_REF]; [START_REF] Borovic | Control Issues for Microelectromechanical Systems[END_REF]). Open loop control strategies are very useful for MEMS [START_REF] Acar | Concept, implementation, and control of wide bandwidth MEMS gyroscopes[END_REF]) but they are most of time not robust w.r.t model uncertainties and nonlinear dynamics. Linear closed loop controllers are often sufficient to ensure the stability and robust performances when the actuator is driven for few micrometers only. The linear quadratic Gaussian control (Boudaoud et al. [2013]) is powerful to filter dynamic noises but requires a full state estimate and is not robust w.r.t model uncertainties and environmental disturbances. Linear time invariant H∞ control [START_REF] Haddab | Commande robuste d'une pince microfabrique actionnement lectrostatique[END_REF]) can improve both positioning precision and bandwidth of the actuator but is limited to satisfy robust performances in a wide operating range due to fixed weights. The H∞ loop shaping controller is useful to control multiple degrees of freedom of MEMS by decoupling several input/output transfers functions and satisfying robust performances [START_REF] Pang | fabrication, sensor fusion, and control of a micro X-Y stage media platform for probe-based storage systems[END_REF]). Direct and indirect sliding mode control with an adaptive gain [START_REF] Fei | Robust Adaptive Sliding Mode Controller for Triaxial Gyroscope[END_REF]) or with an observer [START_REF] Fei | Adaptive sliding mode control with sliding mode observer design for a MEMS vibratory gyroscope[END_REF]) proved its efficiency to control comb drive actuators in MEMS gyroscopes. Closed loop stability and robust performances are satisfied with an on-line identification of the upper bounds of uncertainties [START_REF] Fei | Robust Adaptive Sliding Mode Controller for Triaxial Gyroscope[END_REF]) or by considering a single non-linearity such as the damping [START_REF] Ebrahimi | Robust sliding-mode control of a MEMS optical switch[END_REF]) for the controller synthesis. To deal with the side instability, linear and nonlinear control approaches have been proposed. They are based on input/output linearization [START_REF] Maithripala | Capacitive stabilization of an electrostatic actuator : output feedback viewpoint[END_REF]), differential flatness [START_REF] Zhu | Improving the performance of an electrostatically actuated MEMS by nonlinear control : advances and comparison[END_REF]), backstepping [START_REF] Salah | Robust Backstepping Nonlinear Control for Parallel-Plate Micro Electrostatic Actuators[END_REF]) or robust PID [START_REF] Vagia | Design of a robust PID-control switching scheme for an electrostatic microactuator[END_REF]). Recent works demonstrated the interest of Linear Parameter Varying (LPV) controllers to deal with this non-linearity [START_REF] Shirazi | An LPV Design Approach for Voltage Control of an Electrostatic MEMS Actuator[END_REF]; [START_REF] Alwi | Sliding Mode Control Design of an Electrostatic Microactuator using LPV schemes[END_REF]). Such controllers allow driving the actuator beyond the side instability (generally less than 10 µm displacement). Nevertheless, they do not provide degrees of robustness to matched nonlinearities such as damping and stiffness variation. The problem of positioning comb drive actuators from few hundred nanometers to several tens of micrometers is not addressed. Moreover, the design of the controllers assumes that the MEMS is working in ideal conditions and the effect of environmental noises on the positioning resolution is not considered. Such issues have not been analyzed previously and remain an open research area although it is crucial for MEMS based microgrippers.

The aim of this paper is to provide a new approach to control MEMS based microgrippers using comb drive actuators. A multivariable controller is designed through the LPV theory in order to ensure the stability and some closed loop performances (precision, nanometer resolution, no overshoot, fast response time) over a very wide operating range despite significant noises and nonlinearities. Nonlinearities arise from a large variation of the stiffness and the damping when increasing the displacement of the actuator. Nonlinear parameters, namely the stiffness and the damping, are matched. Therefore instead of using them "independently" as scheduling variables, the use of only one scheduling parameter is proposed. To do so, a multivariable polynomial LPV model of the system is introduced. The polynomial parametric dependence of the model causes some issues for the derivation of the controller. As such, an affine LPV descriptor model and a finite number of linear matrix inequalities are used for the controller synthesis. In order to achieve resolution requirement, the main source of noise to be considered for the controller synthesis is defined taking into account results from previous works [START_REF] Boudaoud | Noise characterization in millimeter sized micromanipulation systems[END_REF]; [START_REF] Boudaoud | Study of thermal and acoustic noise interferences in low stiffness AFM cantilevers and characterization of their dynamic properties[END_REF]). A specific weighting function is introduced to ensure a high positioning resolution.

The relevance of the control approach is demonstrated considering a commercial microgripper (FT-G100 from FemtoTools Company) and shows that large positioning with nanometer resolution can be achieved using two sensors only and does not require the use of expensive vibration isolation tables commonly employed in microrobotics which is a great advantage from application point of view. The paper is organized as follow. In section 2, the global architecture of the microgripper is presented. Sections 3 and 4 deal with the nonlinear modeling and the LPV modeling of the actuation mechanism. The proposed control strategy is presented in the fifth section. In section 6, experimental control results are presented.

FEATURES OF THE MICROGRIPPER

The FT-G100 microgripper is designed to handle objects ranging from 1µm to 100µm. To pick up an object, an actuated arm is pushed toward closure thanks to a comb drive actuator. The base of the actuated arm is fixed to a flexure joint Fig. 1 and a suspension mechanism including two pairs of clamped-clamped beams holds the movable part of the actuator. While the gripping arms are closed around an object, the deflection of the sensing arm is detected by a capacitive sensor. The sensor consists of a transverse comb-drive with a differential capacity proportional to the displacement of the movable fingers. This displacement (due to the gripping force) is translated into analog voltage V out throughout a MS3110 readout chip (Irvine Sensors). Two pairs of clamped-clamped beams and a flexure joint are also attached to the sensor and the base of the sensing arm respectively. The FT-G100 microgripper is used in many robotics laboratories. Very few control solutions are nevertheless provided by the manufacturer. In this work, the designed controller can be applied to any MEMS based microgripper using a comb drive actuator [START_REF] Yamahata | Silicon Nanotweezers with Adjustable and Controllable Gap for the Manipulation and Characterization of DNA Molecules[END_REF]; [START_REF] Beyeler | Monolithically fabricated microgripper with integrated force sensor for manipulating microobjects and biological cells aligned in an ultrasonic field[END_REF]; [START_REF] Muntwyler | Monolithically Integrated Two-Axis Microtensile Tester for the Mechanical Characterization of Microscopic Samples[END_REF]; [START_REF] Vijayasai | Haptic controlled threeaxis MEMS gripper system[END_REF]) and leads to new perspectives for high precision micro/nano-manipulation tasks.

NONLINEAR MODELING AND IDENTIFICATION

The aim is the accurate position control of the actuated arm tip y a (L) in the direction of y axis (see Fig. 2). In this section, a nonlinear dynamic model of the actuation mechanism is proposed. The model describes the dynamic transfer between the input voltage V in and the position y a (L). Input noises are also considered. Static and dynamic parameters of the model are identified experimentally and a final nonlinear model is derived for the control.

Nonlinear dynamic modeling

The modeling approach is performed under the following assumptions (see Boudaoud et al. [2013] for assumptions verification): (i) the actuated arm is rigid in the direction of y axis, (ii) the slender beam between the shuttle and the actuated arm is not subject to buckling and (iii) the pseudo rigid beams behaves as a hinge joint. Therefore, let us consider a nonlinear Euler Bernoulli equation describing the undamped bending of the suspension mechanism at the midpoint of the doubly clamped beams y a (x ea ) in the direction of y axis (see Fig. 2.a):

EI s ∂ 4 y a (x ea ) ∂x 4 + ρl s e s ∂ 2 y a (x ea ) ∂t 2 (1) + N ∂ 2 y a (x ea ) ∂x 2 = F elec δ (x -L s )
where E is the Young's modulus, I s is the area moment of inertia of the doubly-clamped beams, ρ is the mass density, F elec is the electrostatic force and δ is the Dirac delta function. Moreover, L s , l s and e s are dimensions as described in Table .1.

The two first terms of equation (1) relates to the linear Euler Bernoulli beam equation. The third term introduces a non-linearity due to the so called inplane tension located at the clamped parts of the suspensions (see Fig. 2.a). The inplane tension N is given by (Kovacicv et al. [2011]):

N = El s e s 4L s 2Ls 0 ∂y a ∂x 2 ∂x (2)
The electrostatic force is governed by the following equation [START_REF] Moussa | Investigating the reliability of electrostatic comb-drive actuators used in micro fluidic and space systems using finite element analysis[END_REF]):

F elec = N a εh z 2g .V 2 in (3)
where V in is the actuation voltage, N a = 1300 is the total number of fingers (movable and fixed) in the actuator, ε = 8.85 pF/m is the permittivity of the dielectric material (air), h z = 50 µm is the thickness of comb fingers, and g = 6 µm is the gap spacing between two fingers.

Considering the assumptions of the modeling, the variable y a (L) is introduced in equation (1) leading to:

EI s D a ∂ 4 y a (L) ∂x 4 + ρl s e s D a ∂ 2 y a (L) ∂t 2 (4) + N D a ∂ 2 y a (L) ∂x 2 = F elec δ (x -L s )
where D a = L / x ea is an amplification parameter (x ea = 1100 µm).

Substituting equation ( 2) into (4) gives the partial differential equation:

EI s D a ∂ 4 y a (L) ∂x 4 + ρl s e s D a ∂ 2 y a (L) ∂t 2 (5) +   El s e s 4L s D a 2Ls 0 ∂y a ∂x 2 ∂x   ∂ 2 y a (L) ∂x 2 = F elec δ (x -L s )
Thus, using the perturbation method (see (Kovacicv et al. [2011]) for more information), equation ( 5) can be simplified into a lumped parameter equation of the first Eigen mode of the system as follows:

m a0 ∂ 2 y a (L) ∂t 2 + k 1a y a (L) + k 3a y 3 a (L) = 1 D a F elec (6)
with:

m a0 = ρl s e s φ 2 1 (2L s )D 2 a 2Ls 0 φ 2 1 (x).∂x k 1a = 2EI s φ 2 1 (2L s )D 2 a 2Ls 0 (φ 1 (x)) 2 .∂x k 3a = E φ 2 1 (2L s )D 2 a ( 2Ls 0 (φ 1 (x)) 2 .∂x) 2
The parameter φ 1 is the shape mode of the first Eigen mode of the suspensions. k 1a is the effective linear stiffness of the actuation mechanism. The nonlinear term k 3a is the cubic stiffness (it is due to the inplane tension N ).

In the partial differential equation (1) and in equation ( 6) , the damping of the system is not taking into account and only the mass of the suspension mechanism (i.e. doubly clamped beams) is considered.

Hence, in order to account the mass of the complete actuation mechanism, the mass M sh is added to m a0 . The parameter M sh is the mass of the shuttle and that of the actuated arm. Moreover, to account the real dynamic behavior of the system, the term d a is introduced to define the damping of the actuation mechanism.

As such, equation (6) becomes:

m a ∂ 2 y a (L) ∂t 2 + d a ∂y a (L) ∂t + k 1a y a (L) + k 3a y 3 a (L) = 1 D a F elec (7)
where m a = m a0 + M sh is the mass of the actuation mechanism.

In the next section, both static and dynamic parameters of the nonlinear model are identified using experimental data from a commercial FT-G100 microgripper.

Identification

A high resolution (0.01 nm) laser interferometer (SP-120 SIOS Metechnik GmbH) is used to perform measurements of the displacement y a (L) at the tip of the actuated arm (see Fig. 3). Input voltages are generated using a dSPACE control board (DS1103). First, k 1a and k 3a are identified through the experimental static characteristic y a (L) / V in for 0 < V in < 200 V . Experimental data are then fitted using equation ( 7) in static mode. Let us recall that in the modeling part, equation ( 6) is related to the Euler Bernoulli equation (1). This modeling approach justify and explain why a cubic stiffness term k 3a appear in the lumped parameter equation ( 6). However, from experimental data, it has been shown that for the FT-G100 actuation mechanism, a third order polynomial (related to the nonlinear stiffness with a cubic term) is not sufficient to describe accurately the static characteristic y a (L)/V in . The mean error of the static characteristic in the whole operating range of the system is found to be equal to 17.11 %. . The reason is that in addition to the inplane tension (which produce the cubic stiffness k 3a ), other sources of nonlinearities such electrostatic ones can produce nonlinear stiffness terms of higher orders [START_REF] Ye | Optimal Shape Design of an Electrostatic Comb Drive in Microelectromechanical Systems[END_REF]). Therefore, in this study, we increase the order of the nonlinear stiffness in order to reduce the mean fitting error to less than 10 %. The nonlinear stiffness is described by a sixth order polynomial leading to the following expression:

k a = 6 i=1 k ia y i-1 a (L) = 1 D a y a (L) F elec (8)
As such, in static mode, the equation governing the displacement y a (L) becomes:

k 1a y a (L) + k 2a y 2 a (L) + k 3a y 3 a (L) + k 4a y 4 a (L) +k 5a y 5 a (L) + k 6a y 6 a (L) = 1 D a F elec (9)
Equation ( 9) is not deduced from equation (6) when t → +∞. Only the parameters y a (L),F elec ,k 1a ,k 3a and D a are the same in the two equations. Other terms are due to the polynomial order extension.

The mean fitting error is then reduced to 3.66 % and results are presented in Fig. 4.a. The identified linear, quadratic, cubic and higher order stiffness terms are given in Table.2

(L)/V in , (b) k a /V in and (c) d a /V in .
To describe accurately the damping of the system, it has been sufficient to model it by a fourth order polynomial (using identification results) of the form:

d a = 4 i=0 d ia y i a (L) (10)
Taking into account the nonlinear terms of the stiffness and the damping, the nonlinear state space model of the actuation mechanism is given as:

                             ẏa(L) ÿa(L) =   0 1 -1 ma 6 i=1 k ia y i-1 a (L) -1 ma 4 i=0 d ia y i a (L)   ya(L) ẏa(L) + 0 Ke maDa .V 2 in ya(L) = 1 0 ya(L) ẏa(L) Ke = Naεhz 2g (11)
To validate the modeling approach based on nonlinear polynomial terms, frequency responses of the model are compared with experimental ones. For each operating point, experimental and simulated step responses are recorded and a Fast Fourier Transform (FFT) is computed. Results are presented in Fig. 5.

The nonlinear model describes accurately the increase of the fundamental resonance frequency with increasing the input voltage due to the nonlinear stiffness. The fundamental resonance frequency increases from 827 Hz to more than 2 kHz for 5 V < V in < 175 V . For V in > 70 V , higher order eigenmodes can be observed experimentally. The model is able to capture this characteristic. For V in > 155 V , the value of the fundamental resonance frequency obtained by simulation does not agree with experiments. This is due to the fact that the actuated arm is close to the sensing arm. Hence, the overshoot of the step response causes a contact between the two arms. This contact is not considered in the model.

Noise influences

In [START_REF] Boudaoud | Noise characterization in millimeter sized micromanipulation systems[END_REF]) and [START_REF] Boudaoud | Study of thermal and acoustic noise interferences in low stiffness AFM cantilevers and characterization of their dynamic properties[END_REF]), a detailed study on environmental noises and their effect on the positioning resolution of micromanipulation systems is performed. Three main sources of noise have been studied: cultural noises (ground vibrations), acoustic noises and thermal noises. Models and experiments allowed the determination of the frequency range and the root mean square (r.m.s) amplitude of each source of noise:

(1) The power spectrum of the cultural noise has a 10 -11 f 4 shape (f is the frequency) (see Fig. 6). In the frequency range [1 Hz 100 Hz], the r.m.s amplitude of the cultural noise is about 223 nm. The r.m.s amplitude reaches 6 nm and 1.2 nm in the frequency ranges [100 Hz 250 Hz] and [300 Hz 450 Hz] respectively.

(2) Acoustic noise can amplify the r.m.s amplitude of the cultural noise by a factor of ten beyond 35 dB Sound Pressure Level (SPL). This threshold of acoustic noise can be easily reached in laboratories environments.

(3) The thermal noise is a pseudo white noise. Its r.m.s amplitude is related to the damping of a mechanical system. The thermal noise is often considered as the main limitation of the positioning resolution in MEMS. In [START_REF] Boudaoud | Study of thermal and acoustic noise interferences in low stiffness AFM cantilevers and characterization of their dynamic properties[END_REF]), it has been found that the positioning resolution of AFM (Atomic Force Microscopy) based microgrippers is ranging from 365 pm to 675 pm. It is shown that resolutions are limited by both thermal and acoustic noises. PSD of the cultural noise on a vibration isolation table 10 -11 /f 4 noise Fig. 6. PSD of the cultural noise measured on and out of a vibration isolation table [START_REF] Boudaoud | Noise characterization in millimeter sized micromanipulation systems[END_REF]).

The cultural noise b(t) considered in this paper relates to the noise measured out of the vibration isolation table . Its PDS has a 10 -11 /f 4 shape.

Both the cultural noise and the thermal noise will be considered as input signals (Fig. 7) to define their effect on the positioning resolution of the actuation mechanism.

Cultural noise

The cultural noise b(t) is introduced at the input of the nonlinear model ( 11) such that the dynamic transfer G y a (L) / b satisfies [START_REF] Younis | MEMS Linear and Nonlinear Statics and Dynamics[END_REF]):

     lim f →0 G y a (L) / b (f ) = 1 G y a (L) / b (f = f ra ) = Q a = d -1 a m a k a ( 12
)
where f ra and Q a are respectively the fundamental resonance frequency and the quality factor of the actuation mechanism.

These conditions allow introducing b(t) in (11) as:

                       ẏa(L) ÿa(L) =   0 1 -1 ma 6 i=1 k ia y i-1 a (L) -1 ma 4 i=0 d ia y i a (L)   ya(L) ẏa(L) + 0 Ke maDa .V 2 in + 0 K bt ma .b (t) ya(L) = 1 0 ya(L) ẏa(L) (13)
with:

K e = Naεhz 2g , K bt = 6 i=1 k ia y i-1 a (L)
Therefore, around an operating point δ a , the Power Spectral Density (DSP) of the displacement ỹa (L) = y a (L)-δ a in response to the cultural noise can be deduced from the following relation:

S cult (f, δ a ) = |K bt G ỹa (f, δ a )| 2 S b (f ) (14)
where

G ỹa (f, δ a ) = 1 m a s 2 + d a (δ a ) s + k a (δ a ) (15)
with:

k a (δ a ) = 6 i=1 ik ia δ i-1 a , d a (δ a ) = 4 i=0 d ia δ i a S cult (f, δ a )
and S b (f ) are respectively the PSD of the displacement ỹa (in response to the cultural noise) and the PSD of the cultural noise.

Thermal noise

The thermal noise is considered as the main source of noise that limits the absolute resolution of MEMS [START_REF] Younis | MEMS Linear and Nonlinear Statics and Dynamics[END_REF]). For this reason, this source of noise is initially considered in the study. In the case of the FT-G100 actuation mechanism, one can describe the PSD of the displacement ỹa (L) = y a (L) -δ a in response to the thermal noise as follow:

S therm (f, δ a ) = |G ỹa (f, δ a )| 2 S ξ (f, δ a ) ( 16
)
where S therm (f, δ a ) and S ξ (f, δ a ) are respectively the PSD of the displacement ỹa (in response to the thermal noise) and the PSD of the thermal noise.

The thermal noise ξ, also called Langevin force, is a white noise. Its PSD is related to the damping of the mechanical system, the temperature T expressed in Kelvin and the Boltzman constant K B = 1.38 × 10 -23 J/Kelvin [START_REF] Callen | [END_REF]).

In the case of the actuation mechanism:

S ξ (f ) = 4.K B .T.d a (δ a ) (17)
Equations ( 15) and ( 16) are then used to define the PSD of the displacement ỹa when the actuation mechanism is subject to the thermal noise:

S therm (f, δ a ) = 4K B T m a k a (δ a )d a (δ a ) × 1 -4π 2 m a k a (δ a ) f 2 2 m a k a (δ a ) d a (δ a ) + 4π 2 m a k a (δ a ) f 2 -1 (18)
To define the effect of the cultural noise and that of the thermal noise on the resolution of the system, S cult (f, δ a The r.m.s. values of ỹa are derived from the PDSs according to:

r.m.s ỹa(L) (f 1 , f 2 ) = f2 f1 P SD (f ) ∆f (19)
With ∆f = f 2 -f 1 = 1 kHz. The r.m.s of ỹa is the absolute positioning resolution.

When the system is subject to the thermal noise only, the absolute positioning resolution of the actuation mechanism is 59.8 pm. Such a displacement cannot be measured by the laser interferometer due to its resolution (0.01 nm). On the other hand, the absolute positioning resolution of the system subject to the cultural noise only is much higher. Results show that r.m.s ỹa(L) (f 1 , f 2 ) is in this case equal to 230.32 nm. As such, the cultural noise is the most relevant source of noise to be considered. Therefore, the final model for the controller synthesis is that given in (13). The controller will be designed taking into account frequential properties of the cultural noise in order to reduce r.m.s ỹa(L) to less than 80 nm in the frequency range [1 Hz 1 kHz].

LPV MODELING

Polynomial LPV model

To design a gain scheduled controller based on LMI and the H∞ theory (see [START_REF] Apkarian | Convex characterization of gain scheduled H∞ controllers[END_REF]; [START_REF] Chilali | H∞ design with pole placement constraints: An LMI approach[END_REF]; [START_REF] Gahinet | Affine parameterdependent Lyapunov functions and real parametric uncertainty[END_REF]; [START_REF] Scherer | Multiobjective output-feedback control via LMI optimization[END_REF]; [START_REF] De Souza | Gain scheduled H 2 controller synthesis for linear parameter varying systems via parameter-dependent lyapunov functions[END_REF])) a LPV model is derived from the nonlinear model ( 13). This nonlinear model can be reformulated into an affine (i.e. linear) LPV model if the stiffness k a and the damping d a (see equations ( 8) and ( 10)) are selected as independent varying parameters for which only the lower and the upper bounds of variation are considered. In this case, well known gain scheduled control strategies based on the LPV/H∞ methodology can be used (see an example in [START_REF] Poussot-Vassal | A new semi-active suspension control strategy through LPV technique[END_REF])).

In the case of comb drive actuators, the nonlinear parameters k a and d a are matched and they have a polynomial structure. Such properties are not taken into account with a classical affine LPV/H∞ design which is very conservative. In the present work, the use of the operating point δ a as the only varying parameter is proposed. Using a Jacobian linearization, the nonlinear plant ( 13) is formulated into a polynomial LPV model of the form:

G a (s, δ a ) :    Ẋp (t) = A p (δ a ) X p (t) + B p U (t) b(t) y a (L, t) = C p X p (t) (20) A p (δ a ) = 0 1 -k1a ma -d0a ma + δ a 0 1 -2k2a ma -d1a ma + • • • • • • + δ 2 a 0 1 -3k3a ma -d2a ma + δ 3 a 0 1 -4k4a ma -d3a ma + • • • • • • + δ 4 a 0 1 -5k5a ma -d4a ma + δ 5 a 0 1 -6k6a ma 0 B p = 0 0 Ke maDa K bt ma , C p = 1 0 T , X p = y a (L) ẏa (L)
where A p ∈ R na×na , B p ∈ R na×m and C p ∈ R p×na , with n a = 2, m = 2 and p = 1.

To simplify notations, from now on, the variable ỹa (L) will be designed by y a (L). The non-linearity arising from the square voltage is overcome by considering U = V 2 in as the input of the model for the controller design. This consideration is performed since only positive values of the input voltage are applied to the actuator.

By referring to Fig. 5, the amplitude of the 2nd eigenmode is at least 18 times smaller than that of the first eigenmode. It frequency is two times higher than that of the first eigenmode. Therefore, we assume that unmodeled dynamics due to the linearization have a very little influence on the closed loop dynamic of the nonlinear system.

In LMI based H∞ theory, controllers are derived form decision variables (i.e. solutions) of a set of LMI conditions depending on the varying parameters. Since varying parameters have an infinite number of values between their lower and upper bounds, deriving a LPV controller would require solving an infinite number of LMIs. To deal with this constraint, a commonly used solution consists in solving LMIs only on the vertices of a polytopic model. This solution can be easily applied when the LPV model of the system is affine w.r.t the varying parameters (see examples in (Biannicet al. [1996]; [START_REF] Poussot-Vassal | A new semi-active suspension control strategy through LPV technique[END_REF])).

For model ( 20), the varying parameter enters in a polynomial way into the system definition which causes some issues for the controller derivation. To deal with this problem, a Linear Fractional Transformation of the model ( 20) is proposed. The polynomial LPV model is translated into an affine LPV descriptor model using a singular implicit modeling by a LFT. Afterwards, a polytopic controller is derived considering control specifications and specific LMI conditions (Fig. 9).

Polynomial quasi LPV model

Affine LPV descriptor model

Gain scheduled polytopic controller LFT LMIs

Noise consideration Fig. 9. Derivation of the gain scheduled polytopic controller from the polynomial LPV model.

Derivation of an affine LPV descriptor model

In order to simplify the parametric dependence of the state matrix A p (δ a ) and to define an affine LPV descriptor model, a LFT (see Fig. 10) is performed from the model G a (s, δ a ) such that:

  Ẋp y a z p   = N 11 N 12 N 13 N 21 N 22 N 23 N 31 N 32 N 33 X p U b v p (21)
where z p ∈ R nz and v p ∈ R nv are the input and output signals of the uncertainty block (i.e. varying parameter δ a ) as shown in Fig. 10, with v p = δ a z p .

The sub-matrices {N i,j } i,j=1,...,3 are defined from the interconnection of the LFT. The sizes of the vectors v p and z p are related to the polynomial orders of the stiffness and the damping. As such, n z and n v are defined by n z = n v = 7 (i.e. order of the highest polynomial + 1).

y a v p z p p X & p X ∫ ( ) , a a G s δ a δ b U U b   =    
G des (s, δ a ) :    E Ẋdes (t) = A des (δ a ) X des (t) + B des U (t) b(t) y a (L, t) = C des X des (t) (22) A des = N 11 N 13 δ a N 31 δ a N 33 + I nv , B des = N 12 N 12 δ a N 32 δ a N 32 C des = [ C p 0 1×nv ] , E = diag (I na , 0 nv )
E is a singular matrix (i.e. non invertible matrix). Moreover:

N 11 = 0 1 -k1a ma -d0a ma + δ a 0 1 -2k2a ma -d1a ma , N 12 =   0 K e D a m a   , N 12 =   0 k a (δ a ) m a   N 13 = δ a 0 0 0 0 0 0 0 -3k3a ma -4k4a ma -5k5a ma -6k6a ma -d2a ma -d3a ma -d4a ma , N 21 = [ 1 0 ] , N 22 = 0, N 23 = 0 1×nv , N 31 =         1 0 0 0 0 0 0 0 0 1 0 0 0 0         + δ a         0 0 1 0 0 0 0 0 0 0 0 1 0 0         , N 32 = 0 nv×1 , N 33 = δ a        
0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0

       
The vector N 12 describes the transfer between the cultural noise b(t) and the state vector X p . The PSD of the cultural noise is predominant at low frequencies (see Fig. 6). As such, the effect of this source of noise on the positioning resolution is the most significant in the linear range of the actuation mechanism. Therefore, for the controller synthesis, N 12 is defined considering a frozen value of δ a which is equal to 5 µm.

The descriptor model is used only for the synthesis of the scheduled controller via the resolution of a finite number of LMIs conditions. The controller is thereafter applied on the real microgripper.

LPV BASED ROBUST CONTROL DESIGN

The controller is designed based on the multivariable descriptor model ( 22). In order to achieve a fast and accurate positioning of the actuated arm in a large operating range (from few hundred nanometers to several tens of micrometers) with a nanometer resolution, control performances are defined as follow:

(1) The closed loop response time of the system must be lower than 30 ms (5 times smaller than the response time of the open loop system).

(2) No overshoot is admitted. An overshoot of few micrometers can destroy an object to be manipulated.

(3) The maximum static error lower than 0.1 %.

(4) The r.m.s of y a (L) must be lower than 80 nm in the frequency range [1 Hz 1 kHz].

In order to satisfy control specifications, two weighting transfer functions W 1 and W 2 are introduced. The weight W 1 allows tracking performances (1), ( 2) and (3) by applying constraints to the sensitivity function of the closed loop system (see Fig. 11). The function W 2 is used for the cultural noise rejection.

ε(t) + _ ( , ) a a G s δ W 2 b (t) W 1 e 1 (t) d (t) ( , ) a K s δ H U ∞ ∫ δ a p X p X & p v p z
Fig. 11. Control scheme with the augmented LPV descriptor model considering the input noise b(t).

The controller is designed based on the multivariable descriptor model ( 22). The cultural noise can be measured in real time by fixing the laser spot of the interferometer on a solid block (see [START_REF] Boudaoud | Noise characterization in millimeter sized micromanipulation systems[END_REF]) for more details). The noise measurement is used by the controller to reduce the effect of b(t) on the position y a (L) of the actuated arm. In this work, the use of a shaped form of the signal b(t) as a second input of the controller for the synthesis is proposed (see Fig. 11). This allows obtaining nanometer resolution while keeping a high bandwidth in closed loop.

Control design

The LPV control problem consists in finding a dynamic multivariable scheduled controller K (s, δ a ) that stabilize the closed loop system and minimize the H∞ norm γ of a performance channel ω → z for any value of the varying parameter δ a between its lower bound δ a -and its upper bound δa . Here, for the LPV design presented in Fig. 11, the control problem is defined considering ω = [y c b(t)] T as the external input, z = e 1 as the controlled output, U H∞ as the control input and [y a (L, t) d(t)] T as the measured variable.

From now on, let us consider the following transfer functions:

Gy a / U H∞ = ya(L) U H∞ , Gy a / d = ya(L) d , K U H∞ / ε = U H∞ ε and K U H∞ / d = U H∞ d .
The relationship between the controlled output z and the external input ω is defined by:

e 1 = Sε / y c W 1 y c -Sε / y c Gy a / d + Gy a / U H∞ K U H∞ / ε W 1 W 2 d ( 23
)
where

Sε / y c = 1 + Gy a / U H∞ K U H∞ / ε -1
is the sensitivity function.

Taking into account the relations ( 23), the control problem described above can be solved if the following conditions are satisfied for any admissible value of the varying parameter:

Sε / y c ∞ γ W 1 ∞ Sε / y c Gy a / d + Gy a / U H∞ K U H∞ / ε ∞ γ W 1 W 2 ∞ (24)
To design W 1 , a reference closed loop transfer function T d is defined based on control specifications:

T d = K d 1 ω 2 d s 2 + 2ξ d ω d s + 1 ( 25
)
where K d = 0.9995, ω d = 500 rad/s and ξ d = 1.5 N s/m.

The weighting function W 1 is designed considering the desired closed loop dynamic behavior of the system. Hence:

W 1 = 1 1 -T d (26)
For the noise rejection, the weight W 2 is selected such that the H∞ norm of the transfer ε / d is reduced by 10 dB in the frequency range [1 Hz 10 Hz] where the cultural noise is the most significant (see Fig. 6) and by 6 dB in the frequency range [10 Hz 350 Hz]. This ensures a positioning resolution of less than 80 nm despite of the cultural noise. Thus, for nominal performances, W 2 is defined as:

W 2 = 6.39 × 10 -8 s 2 + 8.75 × 10 -7 s + 3 1 × 10 -7 s 2 + 1 × 10 -7 s + 1.5 (27)
Now, considering the LPV descriptor model ( 22) and the weight functions, a generalized augmented plant description is introduced:

P des (s, δ a ) :        Eg Ẋg(t) = Ag (δa) Xg(t) + B 1g ω(t) + B 2g U (t) z(t) = C 1g Xg(t) + D 11g ω(t) + D 12g U (t) ya(L, t) d(t) = C 2g Xg(t) + D 21g ω(t) (28)
where :

A g ∈ R r×r , B 1g ∈ R r×nw , B 2g ∈ R r×1 , C 1g ∈ R nz×r , C 2g ∈ R 2×r , D 11g ∈ R nz×nw , D 12g ∈ R nz×r , D 21g ∈ R 2×nw , with r = n a + n v + n w
The real n w is equal to the order of W 1 + the order of W 2 . Moreover,

E g = diag [ E, I nw ].
In order to derive the controller K (s, δ a ) satisfying required closed loop performances, it is necessary to use specific LMI conditions that takes into account that the matrix E g is singular. Among LMIs that meet this problem, one can use those proposed in [START_REF] Masubuchi | Gain scheduled controller design based on descriptor representation of LPV systems: Application to flight vehicle control[END_REF]) given as:

Y E T g E g E T g E T g χ = Y E T g E g E T g E T g χ T 0   M A + M T A M B M T C M T B -γI M T D M C M D -γI   < 0 ( 29
)
where

M A = A g (δ a ) Y T + B 2g F T A g (δ a ) H T χ T A g (δ a ) + G T C 2g , M B = B 1g χ T B 1g + G T D 21g , M c = C 1g Y T + D 12g F T C 1g et M D = D 11g
The matrices χ and Y are Lyapunov matrices which are considered constant herein. Decision variables of LMIs (29) are the matrices χ, Y , F , G and H with appropriate dimensions.

As emphasized in Section.4.1, the LMI problem (29) results in an infinite set of LMIs to solve since δ a ∈ δ a δa . Nevertheless, thanks to the affine dependence of A des (δ a ) w.r.t the varying parameter, the LMI can be solved only on the vertices (see [START_REF] Zin | An LPV/Hinf Active Suspension Control for Global Chassis Technology: Design and Performance Analysis[END_REF]) for an example) of the set δ a ∈ δ a δa . Here, two sets of LMIs must be solved because one varying parameter is considered. If the LMI problem is feasible, the solution leads to a LPV polytopic controller with the following state space representation:

K (s, δ a ) = α 1 (δ a ) A k1 B k1 C k1 0 +α 2 (δ a ) A k2 B k2 C k2 0 (30)
The couple {α 1 (δ a ) , α 2 (δ a )} is defined by:

α 1 (δ a ) = δa-δa δa-δa - et α 2 (δ a ) = 1 -α 1 (δ a )
For the computation of the matrices {A ki , B ki , C ki } i=1,2 , the generalized augmented descriptor plant (28) is translated back into a polytopic form. At each vertex of the polytope, the affine descriptor model must be translated into a polynomial form. To address this requirement, the following transformation is performed: for a vertex i, the affine LPV descriptor model that can be represented by the form

P desi : [A gi ] [B 1gi ] [B 2gi ] [C 1gi ] [D 11gi ] [D 12gi ] [C 2gi ] [D 21gi ] 0 =    A 11gi A 12gi A 21gi A 22gi B 11gi B 12gi B 21gi B 22gi [ C 11gi C 12gi ] [D 11gi ] [D 12gi ] [ C 21gi C 22gi ] [D 21gi ] 0    (31)
is equivalent to the polynomial LPV model

P poli : A si B 1si B 2si C 1si D 11si D 12si C 2si D 21si 0 = A 11gi B 11gi B 21gi C 11gi D 11gi D 12gi C 21gi D 21gi -   A 12gi C 12gi C 22gi   A -1 22g i [ A 21gi B 12gi B 22gi ] (32)
This transformation is feasible if the following conditions are satisfied:

(1) B 12gi = 0 and C 12gi = 0, (2) B 22gi = 0 and C 22gi = 0;

(3) A 22g i is an invertible matrix.

The conditions (2) and ( 3) are satisfied in this work. However, the first condition is not satisfied. Therefore, in the generalized state vector X g of ( 28), the state vector X W 1 of the weight function W 1 and the state vector X W 2 of the weight function W 2 are placed from the line n a + 1. In other words, initially X g was defined by

X g = X des X W 1 X W 2 =     X p v p X W 1 X W 2    
and now, it becomes

X g =     X p X W 1 X W 2 v p    
Finally, the matrices {A ki , B ki , C ki } i=1,2 are computed from the decision variables of the LMIs and from the matrices {A si , C 2si , B 2si } i=1,2 of the polynomial model. The controller is implemented into a transfer function form.

Control implementation

In this section, we first compare the performance of two LTI controllers with the performance of the proposed gain scheduled controller in ideal conditions (low noise condition). This work shows the importance of using a scheduled controller in the case of the microgripper. Thereafter, the input noise b(t) is taken into account for the design of the multivariable scheduled controller aiming at satisfying precision requirements.

Low noise condition

Control strategies are implemented considering that the microgripper is working on a vibration isolation table. Single input/output type controllers are sufficient in this case (the noise input is not considered). In addition to the gain scheduled controller, a LTI/H∞ controller and a PID controller are designed for comparison.

For the scheduled controller, the control problem is solved using the LMIs (29). The gain K bt is set at zero and only the weight W 1 is taken into account in the augmented LPV descriptor model (see Fig. 12). Here, the control problem is defined considering ω = y c as the external input, z = e 1 as the controlled output, U H∞ as the control input and y a (L, t) as the measured variable (see Fig. 12). The LMIs are then solved on each vertex of the polytopic model using the Yalmip interface [START_REF] Lofberg | YALMIP : A Toolbox for Modeling and Optimization in MATLAB[END_REF]) with the solver SeDuMi [START_REF] Sturm | Using SeDuMi 1.02, a matlab toolbox for optimization over symmetric cones[END_REF]). For the resolution, the bounds of δ a have been defined such that δ a -= 1 µm and δa = 80 µm. When considering this set, the optimal obtained performance indice γ is equal to 2.73 ensuring a modulus margin greater than 1 γ = 0.36. Therefore, two 4 th order controllers

G ki {A ki , B ki , C ki , 0} i=1,2 are obtained.
The frequency responses of the sensitivity function Sε / y c = ε(t) yc(t) is computed and is compared with the frequency response of 1 W1 . The computation is performed considering 5 µm < δ a < 80 µm. Results are pre- The LTI/H∞ controller is designed using the polynomial LPV model ( 20) for a frozen value δ a = 5 µm and for K bt = 0. The controller is based on the scheme presented in Fig. 12 (without the scheduling part). In other words, the LTI/H∞ control problem is defined considering also ω = y c as the external input, z = e 1 as the controlled output, U H∞ as the control input and y a (L, t) as the measured variable. The same weighting function W 1 is used for the controller synthesis. As a result, a fourth order controller is obtained. The PID controller is designed using the nonlinear model ( 11) and its parameters are identified using Ziegler-Nichols method. For an input reference y c = 5 µm, the ultimate gain K u for which the output of the closed loop system oscillates with a constant amplitude is K u = 8.91. The period of such oscillations is T u = 0.9 ms. Therefore, the proportional, the derivative and the integral gains of the controller are respectively K p = 5.34, K d = 1.12×10 -4 and K i = 2.22 × 10 3 . Real time implementation results of the three controllers are presented in section 6.2.

y c y a (t) ε(t) + _ W 1 e 1 (t) ( , ) a K s δ H U ∞ ( , )

Significant noise condition

In this condition, the microgripper is considered to be working out of a vibration isolation table. Therefore, we consider the multivariable control problem taking into account noise influence. The LMIs (29) are solved on each vertex of the polytopic model derived from (28). For the resolution, the bounds of δ a have been defined such that δ a -= 1 µm and δa = 85 µm.

When considering this set, the optimal obtained performance indice γ is equal to 5.65 ensuring a modulus margin greater than 1 γ = 0.17. Therefore, two 6 th order controllers G ki {A ki , B ki , C ki , 0} i=1,2 are obtained. Each controller has two inputs [ε(t) b(t)] and one output U H∞ . In their transfer function form, the controllers are given in equations ( 33), ( 34), ( 35) and (36) (see also Fig. 18).

G k11 = 7.42 × 10 5 s 5 + 1.25 × 10 9 s 4 + 6.89 × 10 13 s 3 + 1.09 × 10 17 s 2 + 1.6 × 10 21 s + 2.38 × 10 24 s 6 + 2519s 5 + 4.65 × 10 9 s 4 + 7.64 × 10 12 s 3 + 2.13 × 10 17 s 2 + 3.19 × 10 20 s + 2.66 × 10 19

(33) G k12 = 7.95 × 10 5 s 5 + 4.33 × 10 9 s 4 + 2.44 × 10 14 s 3 + 5.27 × 10 17 s 2 + 9.6 × 10 21 s + 1.4 × 10 25 s 6 + 4.26 × 10 5 s 5 + 5.89 × 10 9 s 4 + 2.82 × 10 13 s 3 + 2.68 × 10 17 s 2 + 3.57 × 10 20 s 1 + 2. The frequency responses of the function

Sε / b = ε(t) b(t)
is computed and is compared with the frequency responses of 1 W1W2 . The computation is performed considering 100 nm < δ a < 90 µm. Results are presented in Fig. 14. Real time implementation results of the controller are presented in section 6.3.

EXPERIMENTAL SETUP AND REAL TIME IMPLEMENTATION

Experimental setup

For the real time implementation of the controllers, the experimental setup of Fig. 15 is used. It is composed of:

-The FT-G100 microgripper.

-A laser interferometer sensor (SP-120 SIOS Mebtechnik GmbH) with two heads. The fist head (head 1) is used for the measurement of the actuated arm tip position. The second head (head 2) allows the real time measurement of the cultural noise.

-A 3 degrees of freedom (dof) positioning system which allows the alignment of the laser spot (head 1) with the gripping arm tip.

-A vibration isolation table.

-A controller board (dSPACE DS1103) with a Real Time Interface (RTI).

-An amplifier used to amplify the voltage signal at the output of the controller board which is limited to 10 V . Each controller is designed through the Matlab/Simulink software and is implemented into the controller board. The RTI allows applying different kind of reference trajectories y c . The laser spot of the interferometer head 1 is fixed on the actuated arm tip. The laser spot of the interferometer head 2 is fixed on a solid block. Measured signals are taken into account in real time by the controller board through a digital/digital interface. The root square of the controller output (i.e. √ U H∞ ) is applied to the comb drive actuator.

Controllers implementation in low noise condition

The controllers designed in section 5.3 are implemented in real time on the FT-G100 microgripper with 20 kHz sampling frequency. In such condition, the 3 dof positioning system is fixed on the vibration isolation table. Only the head 1 of the interferometer is used for the measurement. The measured signal is used as a feedback signal.

The sub figures (a), (b) and (c) of Fig. 16 show the normalized (to unity) controlled positions y a (L) for step references y c in the set [5 µm 65 µm]. Results are obtained with the scheduled controller Fig. 16.a, the LTI/H∞ controller Fig. 16.b and the PID controller Fig. 16.c. With the scheduled controller, the worst (i.e. highest) response time is 13.64 ms, there is no overshoot and the highest static error is equal to 0.06 %. With the LTI/H∞ controller, the worst response time is much bigger (38.92 ms), there is no overshoot and the highest static error is equal to 0.02%. With the PID controller, the worst response time is 81.12 ms.

It is clear that the scheduled controller is essential to satisfy robust performances over a large operating range of the system. In Fig. 17, position of the gripping arm y a (L) in response to a reference trajectory are presented (experimental results). The proposed gain scheduled control strategy gives in this case satisfactory results.

In the next section, experimental controlled positioning results of the actuated arm are shown with the LPV controller taken into account the noise influence. Experiments are performed out of the vibration isolation table.

Significant noise condition

The scheduled controller designed in Section.5.4 is implemented on the microgripper which is located out of the vibration isolation table. The experimental setup is as shown in Fig. 15. The position y a (L) is used as a feedback signal and also as scheduling variable for the controller (see Fig. 18). The signal b(t) is the experimental cultural noise measured in real time by the interferometer head 2. The multivariable controller faced real time implementation constraints. Since two digital/digital interfaces are used and since the controller is of high order (transfer functions of 6th order), it was not possible to implement the controller with a sampling frequency higher than 5 kHz. As such, the transfer functions (33), ( 34), ( 35) and ( 36) have been reduced into a third order. As such, the maximum sampling frequency that can be used is equal to 17 KHz. Experiments show that the response time of the closed loop system is higher than the one measured with the single input single output LPV controller. However, for high amplitude displacements, the response time of the multivariable LPV controller remains smaller than of the LTI controllers.

To assess the effectiveness of the scheduled controller for small displacements, input references y c have been applied with amplitudes equal to 100 nm, 250 nm, 500 nm, 750 nm, 1 µm, 5 µm, 10 µm, 15 µm and 20 µm. Results of the controlled displacements y a (L) are presented in Fig. 19 and Fig. 20. For each controller, experiments have been performed considering the case when the input b(t) of the controller is taken into account (referred to "with noise rejection"), and the case when this input is not considered (i.e. the input b(t) of the controller is set to zero) which is referred to "without noise rejection". Results show that it is crucial to consider the noise signal b(t) when the system is subjected to important noise perturbation which is the case in robotics laboratory. This is fundamental for nanopositioning operations. When considering the noise input in the controller, some vibrations are undamped in the controlled position of the gripping arm. This is due to the measurement delay of the digital/digital interface. In Fig. 21, nanopositioning simulation results are demonstrated using the multivariable LPV controller. The input noise is that measured experimentally in [START_REF] Boudaoud | Noise characterization in millimeter sized micromanipulation systems[END_REF]). In this case, since no delay is introduced, vi-

CONCLUSIONS

In this paper, a new control approach of MEMS based electrostatic microgrippers is proposed. The controller allows positioning the gripping arm from few hundred nanometers to several tens of micrometers with a nanometer resolution and a fast response time despite of nonlinearities and significant noises. Such performances are greatly needed in biomedical microrobotics applications and have never been achieved previously with MEMS based microgrippers. To this end, a nonlinear dynamic model of the actuation mechanism has been proposed. First, it has been shown experimentally that the damping of the system is nonlinear and that the use of a third order polynomial to describe the nonlinear stiffness of a comb drive actuator is limited. Therefore, both the stiffness and the damping have been described by polynomials of higher orders. Afterward, two sources of noise that reduce the resolution of MEMS have been studied and the most relevant source of noise to be considered for the control design has been identified and introduced into the nonlinear model. Therefore, the nonlinear model has been reformulated into a polynomial LPV model. To deal with the particular structure of the LPV model, the later has been reformulated into an affine LPV descriptor model and a multivariable gain scheduled control strategy has been proposed. The effectiveness of the controller for noise rejection has been demonstrated experimentally and simulations showed that control results can be improved for nanopositioning if the measurement delay is canceled. The effectiveness of the controller for high amplitude displacements, in terms of precision and closed loop response time, has been demonstrated with experimental arguments. The control strategy allows new perspectives for high precision and high speed micromanipulation tasks in noisy laboratories environments.

Fig. 1 .

 1 Fig. 1. Structure of the FT-G100 microgripper (Femto-Tools GmbH).

Fig. 2 .

 2 Fig. 2. Simplified scheme of the actuation mechanism (a) and equivalent scheme of suspensions (b).

Fig. 3 .

 3 Fig. 3. Measurement of the actuated arm position using a laser interferometer sensor.

Fig. 4 .

 4 Fig. 4. Nonlinear characteristics of the FT-G100 actuation mechanism (a) y a (L)/V in , (b) k a /V in and (c) d a /V in .

Fig. 5 .

 5 Fig. 5. Frequency responses of the nonlinear actuation mechanism: experimental data (continuous line) and simulation data (dashed line).

  /Hz) S b (f) : PSD of the cultural noise out of a vibration isolation table

Fig. 7 .

 7 Fig. 7. The actuation mechanism subject to an input voltage V in , an input 10 -11 /f 4 noise (cultural noise b) and an input white noise (thermal noise ξ).

4 .Fig. 8 .

 48 Fig. 8. PSD of the cultural noise (S b (f )) and PSD of ỹa (L) in response to the cultural noise (S cult (f, δ a )) and in response to the thermal noise (S therm (f, δ a )) for δ a ∈ [5 µm 100 µm].

Fig. 10 .

 10 Fig. 10. Linear Fractional Transformation of the polynomial LPV model. To transform the polynomial LPV model to an affine descriptor model, a new state vector X des = X T p v T p

Fig. 12 .

 12 Fig. 12. Control scheme with the augmented LPV descriptor model considering the input noise b(t) =0.HereK bt = 0 in G a (s, δ a ).sented in Fig.13. The obtained bandwidths are of about 500 rad/s.

Fig. 13 .

 13 Fig. 13. Frequency responses of the sensitivity function (dynamic transfer ε/y c ) and that of the function 1 W1 .

Fig. 14 .

 14 Fig. 14. Frequency responses of the dynamic transfer ε/b and that of the function 1 W1W2 .

Fig. 15 .

 15 Fig. 15. Experimental setup when the microgripper is working out of the vibration isolation table.

Fig. 16 .Fig. 17 .

 1617 Fig. 16. Experimental normalized (to unity) step responses of the controlled gripping arm at different operating points (from 5 µm to 65 µm) with the scheduled controller (a), the LTI/H∞ controller (b) and with the PID controller (c).

Fig. 18 .

 18 Fig. 18. Block diagram of the closed loop control system.

Fig. 19 .

 19 Fig. 19. Controlled positioning (experimental results) of the gripping arm for different input reference (100 nm, 250 nm, 500 nm, 750 nm, 1 µm, 5 µm, 10 µm, 15 µm and 25 µm) with the multivariable scheduled controller.

Fig. 20 .Fig. 21 .

 2021 Fig. 20. Controlled positioning (experimental results) of the gripping arm for 5 µm input reference with the multivariable scheduled controller. brations are damped. Experimental results (with the delay of the measurements) show nevertheless very good results comparing with results obtained without considering the noise input.

Table 1 .

 1 Dimensions of the actuation mechanism

Actuated arm Doubly-clamped beams Length L = 5150 µm Ls = 915 µm (half length) Width l = 150 µm ls = 50 µm Thickness e = 50 µm es = 8 µm

ACKNOWLEDGEMENTS

This work has partially been supported by the French national project NanoRobust (ANR 2011 NANO 006) and by Hamecmopsys.


Philippe Darreau
email: .philippe.darreau@unilim.fr

Francois Pigalle
email: francois.pigalle@unilim.fr

Shelby Woods

Pierre-Olivier Gourinchas

Converge or integrate? A note on Gourinchas and Jeanne : The elusive gains from international financial integration

explain that the gains from capital market integration are small because the natural convergence of economies would have "done the work" of integration if it had not occurred. We provide a simple illustration of this standard theoretical argument using the simplest Solow model in a small open economy.

Introduction

"Standard theoretical arguments tell us that countries with relatively little capital benefit from financial integration as foreign capital flows in and speeds up the process of convergence". In a very influential paper [START_REF] Gourinchas | The elusive gains from international financial integration[END_REF] explain that these gains from capital market integration are low because the natural convergence of economies would have "done the work" of integration if it had not occurred. We provide a simple illustration of their argument using the graph of [START_REF] Schröder | Neo-classical Growth in an Open Economy: A Note[END_REF], the simplest Solow model in a small open economy. Section 2 shows that the gains of integration are small on average. Section 3 shows that the distortion of autarky does not disappear much more quickly, on average, through pure convergence than through integration.

Gains from integration are not very large on average

Consider the Solow model of manuals. The production function is neoclassical, and inputs are capital and labor. The labor growth rate is n, the technical progress growth rate is x, capital depreciation is δ and savings rate is s. There is a parametric similarity and thus all countries of the world have the same steady state. GDP per unit of effective labor is q = f (k) where k is the capital per unit of effective labor. National wealth in units of effective labor is a = k + e where e is borrowing (e < 0) or loaning (e > 0) abroad. GNP per effective unit of effective labor is y = q + Re = w + Rk + Re = w + Ra where R = f ′ (k) is the gross interest rate.

Like Gourinchas and Jeanne we takes a model of small open economy. This framework was illustrated by Shröeder (1972). Consider Figure 1 a small open economy with inital wealth k 0 and rest of the world at steady state k * . The small economy can either remain in autarky and converge "naturally" to its steady state, or financially integrate with the rest of the world and transit to its steady state as a small open economy.

If it does not integrate, the country will converge toward the steady state k * . Its capital will increase gradually to k 0 to k * . Its consumption will increase to BD to F E (see Figure 1). Gradually the gross interest rate R 0 , will decrease to reach R * .

If it is integrated, the country benefits from R * immediately. The straight line y = w * + R * a, tangent to q = f (k) in F has the slope R * and ordinate w * . The straight line sy tangent to sq in E has the slope sR * and ordinate sw * . As R 0 > R * the small economy borrows e 0 < 0, and its capital immediately becomes k * . Its GDP is immediately q * , and its GN P is y 0 . The country would pay -R * e 0 in interest. The country benefits from higher wages but suffers from a lower interest rate, but the immediate effect on income is positive since y 0 > q 0 . This result is due to the concavity of the production function. Consumption increases immediately from BD to AC. Using a Cobb-Douglas function q = k α like Gourinchas and Jeanne, we can calculate, at the moment of opening, the percentage increase in consumption :

c int 0 -c aut 0 c aut 0 = (1 -s)(y 0 -q 0 ) = (1 -s)q 0 = q * + R * e 0 -q 0 q 0 = k * α + αk * α-1 (k 0 -k * ) k α 0 -1 c int 0 -c aut 0 c aut 0 = (1 -α) k * k 0 α + α k * k 0 α-1 -1 (1)
Figure 1: Gain from integration Assuming, as Gourinchas and Jeanne, that α = 0.3 and that steady state capital k * is 2.461 times the initial capital k 0 , we get the instant gain: (1-0.3)(2.46) 0.3 +0.3(2.46) -0.7 -1 = 7.68%. For a given deviation from the steady state, this gain is an increasing then decreasing function of α.

Over the whole period consumption increases by the difference between the integral of the triangle ABF (increased income) and the integral of the triangle CDE (increased savings). The gain is on average :

average gain = 1 k * -k 0 k * k 0 (1 -α) k * k 0 α + α k * k 0 α-1 -1 dk average gain = k * k 0 α-1 2 k * k 0 1-α -α k * k 0 -1 -(2 -α) 1 -k * k 0 -1 (2 -α) -1 (2)
Assuming, as G&J that α = 0.3 and that the steady state capital is 2.46 times the initial capital, we get an average gain of 2.01% 2 .

Figure 1 displays the G&J argument: even if the benefit of integration is initially significant, it is small on average. The area of the triangle ABF is small. The surface is small especially if k is close to the steady state (k * /k 0 low) and depends on the concavity of the production function (α). This graph also shows that by integrating, the economy increases its savings by CDE with respect to autarky. This is true in the Solow model where the savings rate is constant, but not in the Ramsey model used by G&J. As savings increase, convergence to the steady state will be faster.

The speed of convergence is not much higher

The dynamic is different in autarky than in an open economy. In autarky "convergence" is done according to the dynamic equation: Dk = s.q -(n + x + δ)k. In a small open economy3 , dynamics will be done according to the equation : Da = s.y -(n + x + δ)a with y = w(k * ) + R(k * )a. Figure 2 shows the capital growth rate under autarky and under integration. We represent s.q/k and s.y/a and the constant (n + x + δ). The difference represents the growth rate. The curve s.q/k is decreasing and convex since the average product of capital is decreasing. The curve s.y/a = s(w * /a) + sR * is decreasing and convex since the denominator of constant wage w * , increases. The slope of s.q/k is (-s.w/k 2 ), and the slope of s.y/a is (-s.w * /a 2 ). The two slopes are the same for k = a * ; the first is less than the second for k < a * and higher for k > a * .

Figure 2: Growth rates Figure 2 shows that the transitional growth rate is higher when there is integration (AD) than under autarky (BD). However, the linear approximation around the steady state is identical for both processes (CD). The speed of convergence is not much higher; it is in both cases equal to β = (1 -α)(x + n + δ).

Conclusion

Gourinchas and Jeanne calibrate the Ramsey model to show small gains from the integration of capital for a poor country compared to the natural convergence. The Solow model does not measure the gain in utility and does not reflect the change in the savings rate. However it provides a simple illustration of the argument.

G&J consider the ratio k/y. Their benchmark is k * /y * =

2.63 and r = Rδ = 5.42%. Note that k * /k 0 = ((k * /y * )/(k 0 /q 0 )) 1/(1-α) and average capital-output ratio (from PWT 6.1) equal to 1.40 so k * /k 0 ≈ 2.46.2 G&J find that the gain is less than 2% between k/y = 1.29 and k/y = 4.38. In Solow's model gains are less than 2% between k/y = 1.41 and k/y = 4.21. See Appendix.

In a small open economy dynamics is very simple. For an explanation of the dynamics in a general equilibrium model, see[START_REF] Darreau | Long-run effects of capital market integration using Solow's model[END_REF].

Appendix

Same figure than G&J p721 for instant and average gains :


Didace Habineza
email: didace.habineza@femto-st.fr

Micky Rakotondrabe

Yann Le Gorrec
email: legorrec@femto-st.fr

Multivariable Generalized Bouc-Wen modeling, identification and feedforward control and its application to multi-DoF piezoelectric actuators

Keywords: Feedforward Control, Multi-Degrees of Freedom Systems, Asymmetric Hysteresis, Couplings, Generalized Bouc-Wen Model

In the literature, the generalized Bouc-Wen model can track precisely asymmetric hysteresis nonlinearity. In this paper, we propose to extend this generalized model to multivariable hysteresis model that can track the nonlinearities in multi-degrees of freedom (multi-DoF) hysteretic actuated systems. In particular, these systems are typified by strong hysteresis couplings. Then, a method of identification of the multivariable hysteresis model is proposed. Finally, based on the inverse multiplicative structure, we propose a multivariable feedforward compensator of the nonlinearity. The proposed approach has been applied to a multi-DoF piezoelectric tube (piezotube) used in scanning probe microscopy and the experimental verification demonstrated its validity in terms of model precision and compensation efficiency.

INTRODUCTION

Hysteresis is very known to be one of the piezoelectric actuators behaviours that compromise the positioning accuracy. Hysteresis control has been largely studied either by using feedback or feedforward control techniques.

In some applications such as precise positioning and micromanipulation, the use of feedback control techniques to remove hysteresis are limited by the lack of convenient sensors having the required performances and convenient dimensions. In fact, embeddable sensors such as strain gauges do not possess the required performances while performant sensors such as optical sensors are very bulky [START_REF] Rakotondrabe | Bouc-Wen modeling and inverse multiplicative structure to compensate hysteresis nonlinearity in piezoelectric actuators[END_REF] [START_REF] Devasia | A survey of Control Issues in Nanopositioning[END_REF]. For this reason, feedforward (open-loop) techniques have been adopted.

Two approaches have used feedforward control hysteresis in piezoelectric systems. The first approach consists in modeling the hysteresis and then apply a model inversion, in order to establish a linear relation between the input and output displacement [11][8]. This approach is called voltage control. The second approach is the charge control where, instead of inverting the hysteresis model, the input to the system is no longer the voltage, but the electric charge provided by a designed electronic circuit [START_REF] Clayton | Hysteresis and vibration compensation in piezoelectric actuators by integrating charge control and inverse feedforward[END_REF]. This latter approach benefits from the linearity between the input charge and the output displacement of the piezoelectric actuator.

⋆ This work is supported by the national ANR-JCJC C-MUMSproject (National young investigator project ANR-12-JS03007.01: Control of Multivariable Piezoelectric Microsystems with Minimization of Sensors).

Hysteresis modeling for smart materials actuators has been widely studied [START_REF] Banks | Hysteresis modeling in smart material systems[END_REF] [START_REF] Sain | Models for Hysteresis and Application to Structural Control[END_REF]. For piezoelectric actuators, the Preisach, the Prandtl-Ishlinskii and the Bouc-Wen models are the most used if the further aim is to calculate a feedforward controller. The two first are based on the superposition of many elementary hysteresis operators (hysterons) [START_REF] Kuhnen | Inverse feedforward controller for complex hysteretic nonlinear-ities in smart-material systems[END_REF] [START_REF] Iyer | Approximate inversion of the Preisach hysteresis operator with application to control of smart actuators[END_REF] [START_REF] Janaideh | Compensation of rate-dependent hysteresis nonlinearities in a piezo micro-positioning stage[END_REF] [START_REF] Ang | Feedforward controller with inverse rate-dependent model for piezoelectric actuators in trajectory-tracking applications[END_REF], while the third lies on a set of two equations made of a number of parameters controlling the shape and the scale of hysteresis [START_REF] Rakotondrabe | Bouc-Wen modeling and inverse multiplicative structure to compensate hysteresis nonlinearity in piezoelectric actuators[END_REF] [START_REF] Ismail | The hysteresis Bouc-Wen model, a survey[END_REF].

In the classical models of the above techniques (Preisach, Prandtl-Ishlinskii, Bouc-Wen), the hysteresis modelled are symmetrical. However, many hysteresis behaviour encountered in piezoelectric actuators like piezo-stack actuators [START_REF] Merry | Identification, control and hysteresis compensation of a 3 DOF metrological AFM[END_REF], piezoelectric multimorph cantilevers [START_REF] Rakotondrabe | Robust Feedforward-Feedback Control of a Nonlinear and Oscillating 2-DOF Piezocantilever[END_REF], piezotube scanners, etc, are non-symmetrical. To model this class of hysteresis, additional models have been introduced (Dahl [START_REF] Xu | Dahl model-based hysteresis compensation and precise positioning control of an XY parallel micromanipulator with piezoelectric actuation[END_REF]), and some of the existing models have been extended or modified: Extended Coleman-Hodgdon [START_REF] Merry | Modeling and compensation of asymmetric hysteresis in a piezo actuated metrological AFM[END_REF], Modified Preisach [START_REF] Jang | Modeling and control of a piezoelectric actuator driven system with asymmetric hysteresis[END_REF], Modified Prandtl-Ishlinskii [START_REF] Jiang | A modified prandtl-ishlinskii model for modeling asymmetric hysteresis of piezoelectric actuators, Ultrasonics, Ferroelectrics and Frequency Control[END_REF] and Generalized Bouc-Wen [START_REF] Song | Generalized Bouc-Wen model for highly asymmetric hysteresis[END_REF].

The asymmetric hysteresis models cited above are not easily extensible to multi-DoF systems. In addition, in the case of multiple degrees of freedom (multi-DOF), the number of parameters to be identified can be very high leading to a very complex calculation and implementation of the compensator. In fact, the Preisach and the Prandtl-Ishlinskii techniques present these disadvantages relative to the Generalized Bouc-Wen technique. Indeed, the two formers being based on the superposition of many elementary hysteresis, extending them to multivariable implies a substantial increase of the parameters. In the Generalized Bouc-Wen technique, the number of parameters will still be conveniently limited if extended in multivariable offering a simplicity of model identification and of compensator calculation and implementation.

In the literature, the study of multivariable hysteresis is not yet well addressed. [START_REF] Yong | Design, identification, and control of a flexure-based XY stage for fast nanoscale positioning[END_REF][16] [START_REF] Xu | Dahl model-based hysteresis compensation and precise positioning control of an XY parallel micromanipulator with piezoelectric actuation[END_REF] focused on the modeling of multivariable hysteresis which are symmetrical. In [START_REF] Merry | Modeling and compensation of asymmetric hysteresis in a piezo actuated metrological AFM[END_REF], the compensation of hysteresis in a multi-DoF piezoelectric actuator has been considered. However, only the direct transfers have been controlled and the couplings were not considered. The works come back to the feedforward control of several monovariable hysteresis. In [START_REF] Merry | Identification, control and hysteresis compensation of a 3 DOF metrological AFM[END_REF][17], couplings were taken into account but their suppression were performed by combining the feedforward and the feedback schemes. In this paper we propose the modeling, the parameters identification and the feedforward control of asymmetrical hysteresis in multi-DoF piezoelectric actuators. The proposed hysteresis model and compensator are multivariable because they account both the direct hysteresis and the coupling hysteresis. In order to benefit from the model simplicity of the Generalized Bouc-Wen model [START_REF] Song | Generalized Bouc-Wen model for highly asymmetric hysteresis[END_REF], we propose in this paper to extend such monovariable hysteresis model into a multivariable hysteresis model. Then, we propose a procedure of identification of its parameters. Afterwards, a new compensator of the multivariable hysteresis is proposed. The compensation technique is based on the inverse multiplicative structure and on the identified model. The advantage of this proposed compensator is: as soon as the model is identified, the compensator is yielded without additional calculation thanks to the structure. Finally, experiments on a piezotube scanner validate the efficiency of the proposed modeling, identification and compensation approach.

The remainder of this paper is as follows. Section-2 reminds the existing monovariable generalized Bouc-Wen asymmetrical hysteresis model. In section-3, its extension into multivariable hysteresis model is presented. The proposed procedure of parameters identification is also detailed in the same section. Section-4 is devoted to the compensator calculation for the multivariable asymmetrical hysteresis model. Finally, the experimental applications and results with a piezotube actuator are detailed in section-5.

MODELING AND IDENTIFICATION OF MONOVARIABLE GENERALIZED BOUC-WEN

Remind of the monovariable Bouc-Wen and Wang-Wen modeling

Hysteresis in structural materials can be modelled by using the set of equations [START_REF] Song | Generalized Bouc-Wen model for highly asymmetric hysteresis[END_REF] which is relation between force and material deformation [START_REF] Song | Generalized Bouc-Wen model for highly asymmetric hysteresis[END_REF].

f (x, ẋ, z) = αk 0 x + (1 -α)k 0 z ż = ẋ[A -|z| m ψ(x, ẋ, z)] (1)
f is the resisting force, x and ẋ represent respectively the deformation and the velocity of deformation, α is the structure stiffness ratio, k 0 is the initial stiffness, and z is a variable denoting the structure inelastic behaviour.

ψ is the hysteresis shape control function which is defined according to the hysteresis model. For Bouc-Wen and Wang-Wen models, ψ functions are represented respectively by ( 2) and (3):

ψ BW = γ + βsgn( ẋz) (2) ψ W W = γ + βsgn( ẋz) + φ[sgn( ẋ) + sgn(z)]
(3) The shape control functions for both Bouc-Wen and Wang-Wen models depend only on two variables ( ẋ and z). This allows them to provide four different phases for the shape control, but with one additional parameter φ for Wang-Wen model.

Remind of the monovariable Generalized Bouc-Wen

In [START_REF] Song | Generalized Bouc-Wen model for highly asymmetric hysteresis[END_REF], authors extend the number of variables on which depend the hysteresis shape and they include also x in addition to ẋ and z. This permits them to propose a novel shape control function depending on three variables and to control hysteresis with six shape phases instead of four as in Bouc-Wen and Wang-Wen modelling. Such control function is presented in (4) where β i are parameters to be identified.

ψ GBW = β 1 sgn( ẋz) + β 2 sgn(x ẋ) +β 3 sgn(xz) + β 4 sgn( ẋ) +β 5 sgn(z) + β 6 sgn(x) (4)
For the case of piezoelectric based actuators, the first equation of (1) relates the driving signal which is the voltage U with the bending displacement y of the actuator.

In addition, the parameter m can be assumed to be equal to 1 due to their elastic structure [START_REF] Low | Modeling of a three-layer piezoelectric bimorph beam with hysteresis[END_REF]. So, by combining the first equation of the set in (1) with the control function in (4), the generalized Bouc-Wen model is:

y = kU -z ż = U A -|z| ψ U, U , z (5)
with

ψ U, U , z =β 1 sgn U z + β 2 sgn U U + β 3 sgn (U z) + β 4 sgn U + β 5 sgn (z) + β 6 sgn (U ) (6)

A NEW MULTIVARIABLE GENERALIZED BOUC-WEN: MODELING AND IDENTIFICATION

Multivariable Generalized Bouc-Wen modeling

In this section, we propose to extend the monovariable model in equation ( 5) into a multivariable hysteresis model with n inputs and n outputs.

Let the following expression be the extension of the monovariable model represented in [START_REF] Tuma | Optimal scan trajectories for high-speed scanning probe microscopy[END_REF]:

    y 1 y 2 . . . y n     y =     k 11 k 12 • • • k 1n k 21 k 22 • • • k 2n . . . . . . . . . k n1 k n2 • • • k nn     K     U 1 U 2 . . . U n     U -     z 1 z 2 . . . z n     z (7)     ż1 ż2 . . . ż n     ż =     A 11 A 12 • • • A 1n A 21 A 22 • • • A 2n . . . . . . . . . A n1 A n2 • • • A nn     A      U1 U2 . . . Un      U -     ψ 11 ψ 12 • • • ψ 1n ψ 21 ψ 22 • • • ψ 2n . . . . . . . . . ψ n1 ψ n2 • • • ψ nn     Ψ(U, U ,z)       U1 0 U2 . . . 0 Un       diag( U )     |z 1 | |z 2 | . . . |z n |     |z| (8)
i.e.

y = KU -z ż = A U -Ψdiag( U )|z| (9)
where |z| ∈ R n is a state vector and A ∈ R n×n , the hysteresis scale control matrix. Ψ is represented by equation ( 10) and corresponds to the multivariable hysteresis shape control matrix, extended from equation (4).

Ψ(U, U , z) = Bdiag(sgn( U z)) + Cdiag(sgn(U U )) +Ddiag(sgn(U z)) + Ediag(sgn( U )) +F diag(sgn(z)) + Gdiag(sgn(U )) (10)
B, C, D, E, F and G are R n×n matrices controlling the shape of the multivariable hysteresis.

Identification of Multivariable Generalized Bouc-Wen

By developping equations ( 7) and ( 8), [START_REF] Kuhnen | Inverse feedforward controller for complex hysteretic nonlinear-ities in smart-material systems[END_REF] can be derived:

                     y i =   n j=1 k ij U j   -z i żi = n j=1 A ij Uj - n j=1 ψ ij Uj |z j | ψ ij = B ij sgn( Uj z j ) + C ij sgn(U j Uj ) + D ij sgn(U j z j ) +E ij sgn( Uj ) + F ij sgn(z j ) + G ij sgn(U j ) (11)
In the above equation, the parameters to be identified are

k ij , A ij , B ij • • • G ij with 1 ≤ i ≤ n and 1 ≤ j ≤ n
The subscript i denotes the considered output inside the vector y, while j the considered input inside the vector U . Furthermore, for i = j, we have the direct hysteresis and for i = j, we have the coupling hysteresis.

To identify the parameters of the new multivariable Generalized Bouc-Wen model hysteresis, the following identification procedure is proposed:

Step 1 (experimental characterization): Apply a periodic (sinusoidal or triangular) input voltage U j (starting by j = 1) to the system and leave all the remaining inputs (i.e. all U k with k = j and 1 ≤ k ≤ n) equal to zero. Capture the output displacement y j due to this voltage U j (direct axis), and capture also the coupling caused by this voltage to the other axis (y k , k = j). In the model, these will correspond to (U j , y i ) (with i = 1, ..., j, ...n) which, from equation [START_REF] Kuhnen | Inverse feedforward controller for complex hysteretic nonlinear-ities in smart-material systems[END_REF], are derived as follows:

       y i = k ij U j -z i żi = A ij Uj -ψ ij Uj |z j | ψ ij = B ij sgn( Uj z j ) + C ij sgn(U j Uj ) + D ij sgn(U j z j ) +E ij sgn( Uj ) + F ij sgn(z j ) + G ij sgn(U j ) (12)
Step 2 (identification of the direct hysteresis parameters): Identify, first, the parameters k jj , A jj • • • G jj from the direct hysteresis (U j , y i ) (with i = j). This is done by using the following equations:

       y j = k jj U j -z j żj = A jj Uj -ψ jj Uj |z j | ψ jj = B jj sgn( Uj z j ) + C jj sgn(U j Uj ) + D jj sgn(U j z j ) +E jj sgn( Uj ) + F jj sgn(z j ) + G jj sgn(U j ) ( 13
)
Step 3 (identification of the coupling hysteresis parameters): Then, identify the parameters k ji , A ji • • • G ji from the coupling hysteresis ((U j , y i ) with i = j). This is done using the following equations:

             y i = k ij U j -z i żi = A ij Uj -ψ ij Uj |z j | ψ ij = B ij sgn( Uj z j ) + C ij sgn(U j Uj ) + D ij sgn(U j z j ) +E ij sgn( Uj ) + F ij sgn(z j ) + G ij sgn(U j ) żj = A jj Uj -ψ jj Uj |z j | (14)
After that, repeat the previous three steps for j = 2, 3, ..., n to find the remaining parameters.

Parameters identification at the individual steps

Step 2 and Step 3 can be performed with classical technique of identification. This individual identification is similar to the monovariable hysteresis model identification. For instance, in [START_REF] Sireteanu | A Genetic Algorithms Method for Fitting the Generalized Bouc-Wen Model to Experimental Asymmetric Hysteretic Loops[END_REF], the monovariable Generalized Bouc-Wen parameters were identified by using techniques such as Genetic Algorithms. In this paper, we have performed the identification in the individual steps (Step 2 and Step 3) by using the non-linear least square Matlab optimization toolbox (lsqnonlin) mentioned in [START_REF] Song | Generalized Bouc-Wen model for highly asymmetric hysteresis[END_REF] to minimize the sum of squared errors over the concerned entire hysteresis loop.

FEEDFORWARD CONTROL OF THE MULTIVARIABLE GENERALIZED BOUC-WEN

Assuming now that the parameters of the multivariable Bouc-Wen hysteresis model in [START_REF] Merry | Modeling and compensation of asymmetric hysteresis in a piezo actuated metrological AFM[END_REF] are identified, we present in this section a (multivariable) compensator for the hysteresis. The compensator is built in a way such that the general feedforward controlled system as depicted in Fig. 1a satisfies the following condition:

y ref = y (15)
where y ref (with

y ref = (y ref 1 y ref 2 • • • y ref n ) T
) is the desired (reference) input vector. The first equation of the multivariable model in (9) becomes therefore:

y ref = KU -z (16
) From ( 16), the compensator is directly yielded as follows which has an inverse multiplicative structure:

U = K -1 (y ref + z) (17)
where the control voltage U is the output of the compensator and the desired displacement y ref is its input. It is worth to notice that this compensator is implementable because K is invertible (non singular). Indeed, when designing a multi-DOF actuator, the first natural expectation is that each input control U i affects a corresponding output displacement and axis y i : in other words, the direct transfers K ii are non-null. The only case where K is singular is when applying a voltage U i , couplings y j (j = i) exist in a such a way that their sum compensates for the the output y i . This case is very rare and in general, designers of actuators will never allow systems having such characteristics.

In details, the compensator of equation ( 17) is written as in ( 18) while the simplified diagram for the multivariable hysteresis compensator with the proposed approach is represented in Fig. 1-b.

    U 1 U 2 . . . U n     U =     K 11 K 12 • • • K 1n K 21 K 22 • • • K 2n . . . . . . . . . K n1 K n2 • • • K nn     -1 K -1      y ref 1 + z 1 y ref 2 + z 2 . . . y ref n + z n      y ref +z (18)
The proposed approach is devoted to multi-DOF piezo- electric actuators. Another advantage is that there is no additional calculation to derive the compensator, except the inversion of K. As soon as the model is identified, the compensator is derived thanks to the inverse multiplicative structure.

EXPERIMENTAL RESULTS AND DISCUSSIONS

Experimental setup

The actuator used for the experiments is a PT230.94 piezotube scanner (from PIceramic company) often used in atomic force microscopy (AFM). This is an actuator designed to provide deflections along with three directions (X, Y and Z). The actuator is made up of PZT material (lead zirconate titanate), four external electrodes (+x, -x, +y and -y), and an inside electrode for ground. +U and -U voltages can be applied on +x and -x electrodes (+y and -y respectively) to obtain tube deflections along the X direction (Y direction). Extension (elongation) along the Z direction is obtained by applying simultaneously +U on the four external electrodes.

The tube deflections are measured by using LC2420 optical displacement sensors (from Keyence Company). These sensors are set to have a resolution of 50nm, an accuracy of 200nm and a bandwidth of 1500kHz.

As PT230.94 operating voltage range is ±250V , a high voltage amplifier is used to amplify the control voltage.

The sensors and the amplifier are connected to the computer through a 1103-dSPACE board (Fig. 2).

The tubular shape of the actuator does not allow linear measurement with the optical displacement sensors. For this reason, a small cube (with faces pointing towards the sensors) is placed on the top of the piezotube.

System characterization, modeling and identification

This section is devoted to the characterization and modeling of the piezotube and the parameters identification.

• Multivariable Hysteresis Characterization:

We apply first a sine voltage U x with an amplitude of 200V and a frequency of 0.1Hz to the x electrodes (+x and -x) and we set U y equal to zero. The resulting displacements X and Y are measured thanks to the optical sensors. Then we plot X versus U x (direct transfer) and Y versus U x (coupling). They are pictured in Fig. 3a-(1) and in Fig. 3a-(3) respectively.

Then, U x is set equal to zero and we apply the same sine voltage but to the y electrodes (+y and -y). The resulting displacements X and Y are once again measured and we afterwards plot X versus U y (coupling) and Y versus U y (direct transfer). The plots are pictured in Fig. 3a-(2) and in Fig. 3a-(4) respectively. Notice that the chosen frequency (0.1Hz) used for the characterization and identification corresponds to the operational frequency when controlling the system.

• Generalized Bouc-Wen parameters Identification and Validation: This is done by following the steps discussed in section III.B. From the experimental hysteresis curves, we have identified the parameters of the multivariable Generalized Bouc-Wen model that fits with the experimental data by following these steps and by applying a nonlinear least square optimization with Matlab.

With the identified parameters, the Generalized Bouc-Wen set of equations ( 7), ( 8) and ( 10) can now be written as:

X Y = 0.1873 0.0432 -0.0193 0.2047 K U x U y - z x z y (19)
We mention that the model parameters are identified only from the external hysteresis loop, i.e from the hysteresis loops captured when the input voltage amplitude is in its maximal operating range (200V ). Fig. 3b-( 1), ( 2), ( 3) and ( 4) picture the comparison between the experimental hysteresis and the simulation of the identified model. In the figures, hysteresis with the external loop and hysteresis obtained with lower voltages are plotted. This comparison shows the good conveniency of the identified model. The efficiency of the Generalized Bouc-Wen modeling is very remarkable in Fig. 3b-( 2) and ( 3) where the asymmetry of the hysteresis is important. This asymmetry cannot be tracked by the classical Bouc-Wen model.

Hysteresis compensation and Comparison between compensated and uncompensated systems

After the model identification and validation, we have implemented in the Matlab-Simulink software the hysteresis compensator of Fig. 1-b. The compensator is built by using the parameters identified in section V-B. When the piezotube actuator is controlled, the output deflections X and Y versus the input references X ref and Y ref are represented in Fig. 4.

Comparing the initial hysteresis in Fig. 3a with the results in Fig. 4 obtained with the compensator, we can clearly see that the proposed compensator has considerably suppressed the hysteresis along X and along Y axis, and reduced the coupling amplitudes. Table .1 summarizes the performances of the piezotube with and without the proposed hysteresis compensator. The hysteresis which were initially 18.8% and 18.2% for the X and Y displacements respectively were reduced to less than 3% and than 1.5% (direct transfers). We also remark that the input-output map of these direct transfers is with unitary gain. Concerning the couplings hysteresis, they were reduced from 31.2% and 34.7% for the X and Y displacements respectively to almost 0%. The amplitudes of the couplings were also reduced. The remaining couplings are due to the fact that the coupling hysteresis of the piezotube actuator were not regular and then the generalized Bouc-Wen technique could not exactly be modelled. This non-regularity involves other kinds of nonlinearities to be combined with the hysteresis, which is not the objective of this paper.

Lissajous curves trajectories tracking

Lissajous trajectories are used in [START_REF] Tuma | Optimal scan trajectories for high-speed scanning probe microscopy[END_REF] as a novel reference inputs for high speed scanning probe microscopy. In this

(4) (3) expected X coupling X coupling Y expected Y h Uy[V] Y[µm] (2)
(a) Hysteresis Characterization: (1) and ( 3), X and Y displacements when only Ux is applied (Uy = 0), and ( 2) and ( 4) when only Uy is applied (Ux = 0). A multivariable generalized Bouc-Wen model for multi-DoF systems with hysteresis has been presented in this paper. The proposed model is applicable to multi-DoF piezoelectric actuators exhibiting asymmetric hysteresis both in the expected (direct) transfers and in the couplings transfers. Then, a procedure to identify the parameters of the multivariable model has been proposed. Afterwards, we have proposed a compensator devoted to the multivariable generalized Bouc-Wen model. Based on the inverse multiplicative structure, there is no additional calculation of the compensator. The ease of implementation of the compensator was also discussed. Finally, the proposed modeling, identification procedure and compensation technique have been applied to a piezoelectric tube scanner often used in atomic force microscopy and working in 2-DoF.

Experimental tests with classical signal references and more complex references trajectories (Lissajous curves) demonstrated the efficiency of the proposed approach.

Future work will be focused on the adaptation/modification of the proposed approach in order to take into account the modeling and compensation of multivariable hysteresis with rate-dependent feature.

Fig. 1 .

 1 Fig. 1. (a): block diagram of a feedforward control of multivariable hysteresis. (b): Generalized Bouc-Wen based multivariable hysteresis compensator. Z(.) denotes a non-linear operator characterizing the second equation of (9)

Fig. 2 .

 2 Fig. 2. Experimental setup description.

  Experimental curves and the simulation of the identified 2-DoF Generalized Bouc-Wen model with different hysteresis loops.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Hysteresis characterization and Multivariable Generalized Bouc-Wen model representation.

Fig. 5 .

 5 Fig. 5. Lissajous curves tracking before and after compensation. Subscripts a and b stand respectively for af ter and bef ore compensation

Table 1 .

 1 Comparison between system responses, before and after compensation

		Before	After
		Compensation	Compensation
	Hysteresis
	(h/H)
	Expected X	12.8/68=18.8%	1.8/60=3.0%
	Expected Y	11.5/63=18.2%	0.9/60=1.5%
	Coupling X	1/3.2=31.2%	0.0/1.2=0%
	Coupling Y	0.8/2.3=34.7%	0.0/1.1=0%
	Couplings
	amplitudes
	Coupling X	3.2 µm	1.1 µm
	Coupling Y	2.3 µm	1.2 µm
	Coupl X/ Exp X	3.2/68=4.7%	1.1/60=1.8%
	Coupl Y / Exp Y	2.3/63=3.6%	1.2/60=2.0%