Mélodie Monteil
email: melodie.monteil@ensta-paristech.fr

Cyril Touzé
email: cyril.touze@ensta-paristech.fr

Olivier Thomas
email: olivier.thomas@ensam.eu

Nonlinear vibrations of steelpans: analysis of mode coupling in view of modal sound synthesis

Steelpans are musical percussions made from steel barrels. During the manufacturing, the metal is stretched and bended, to produce a set of thin shells that are the differents notes of the instrument. In normal playing, each note is struck, and the sound reveals some nonlinear characteristics which give its peculiar tone to the instrument. In this paper, an experimental approach is first presented in order to show the complex dynamics existing in steelpan's vibrations. Then two models, based on typical modal interactions, are proposed to quantify these nonlinearities. Finally, one of them is observed in free oscillations simulations, in order to compare the internal resonance model to the steelpan vibrations behaviour in normal playing. The aim is to identify the important modes participating in the vibrations in view of building reduced-order models for modal sound synthesis.

INTRODUTION

Steelpans belong to a musical instruments family coming from the island of Trinidad and Tobago. They are usually played in steelbands (see Fig. 1(a)), that are orchestras composed of steelpans covering a range of several octaves.

A steelpan is a tuned percussion, built from cylindrical steel barrels subjected to several stages of metal forming that stretch and bend the structure. The top of the barrel is pressed, hammered, punched and burnt in order to obtain a sort of main thin bowl within which convex substructures are formed. Each convex dome corresponds to a musical note, which natural frequency is precisely tuned according to harmonic relationships ( f , 2f , 3 (or 4) f , ...). Usually, this instrument is played by striking each note with a stick covered by a piece of rubber, as shown on the Fig. 1(b). When a note is stroke with a stick, vibration amplitudes are such that geometric nonlinearities cannot be neglected, and are recognized as a key feature for explaining the peculiar tone of the steelpan [START_REF] Achong | Mode locking on the non-linear notes of the steelpan[END_REF]. This nonlinearity combined with harmonic relationships, due to steelpan tuning, can activate internal resonances creating strong energy transfers between eigenmodes. Vibrations modelling have been proposed in a series of papers by Achong et al. [START_REF]The steelpan as a system of non-linear modelocalized oscillators, I: Theory, simulations, experiments and bifurcations[END_REF]. In these works, the steelpan is considered as a nonlinear system of oscillators, and energy transfers between normal modes of vibration are highlighted. Rossing et al. [START_REF] Rossing | Vibrational mode shapes in Caribbean steelpans. I. Tenor and double second[END_REF][START_REF] Rossing | Vibrational mode shapes in Caribbean steelpans. II. Cello and bass[END_REF] have performed modal analyses by holographic techniques to observe modal interactions between harmonically tuned notes. More recently, numerical modal analyses with the finite element methods have been proposed [START_REF] Gay | Finite Element Modelling of Steelpan Acoustics[END_REF]. The steelpans sound radiation has been recently addressed experimentally by Copeland et al. [START_REF] Muddeen | Sound radiation from Caribbean steelpans using nearfield acoustical holography[END_REF]. Finally, some metallurgical issues during the steelpan making have been considered in [START_REF] Murr | Fabrication of aluminum, Caribbeanstyle, musical pan: Metallurgical and acoustical characterization[END_REF].

In this contribution, we propose a refined analysis of modal couplings and energy exchanges occurring in nonlinear vibrations of steelpan. A modal analysis first reveals the appearance of pairs of modes from the second harmonic, for each note of the pan. These pairs of modes, having nearly equal eigenfrequencies, are interpreted as a consequence of the localization of the vibration into the notes area. Secondly, experiments in forced vibrations allow to reveal the complex nature of the energy exchanges between the modes, that are excited for very small levels of vibration amplitudes. Simple original models including a 1:2:2 and a 1:2:4 resonance are then fitted to the experiments, showing that: (i) mode pairing substantially complexifies the dynamical behaviours and favours the appearance of unstable regimes, (ii) few details of the resonance curves are missed by those 3-dofs models, indicating that even at very small amplitudes of vibrations, a complex dynamics involving more than 3 modes is at hand. Those findings are used to derive oscillator models for sound synthesis. The first results with three dofs show that the main features (energy transfers and enveloppe modulation) are recovered.

STEELPAN TUNING AND LINEAR ANALYSIS

In steelpan making, a fine tuning of the notes already formed is performed at the end of the building process. Most of the time, the steelpan maker tunes the first three harmonics of each note [START_REF] Kronman | Steel pan tuning: a handbook for steel pan making and tuning[END_REF]. He begins with the fundamental frequency (the pitch) by modelling the center of the note. Then, he focuses on the overtones. He adjusts the frequency of two upper harmonics (partials) by modelling the periphery of the note area, as it is shown of Fig. 2(a). This procedure can be easily understood by considering that each harmonic of a given note is associated to a particular vibration mode shape, as it will be explained in the following.

Modal analysis is usually used to characterize the linear behaviour of a structure by identifying eigenfrequencies, mode shapes and modal damping coefficients. A homemade non-contact coil/magnet exciter is used to excite the steelpan at a given point. The equivalent point force is estimated by recording the current intensity in the coil [START_REF] Thomas | Non-linear vibrations of free-edge thin spherical shells: modal interaction rules and 1:1:2 internal resonance[END_REF]. The steelpan vibratory response, in velocity, is measured with a laser vibrometer.

The steelpan, shown on the figure 2(b), is a right barrel of a double second (middle-high frequency steelpan). It is composed of 19 precisely tuned notes, distributed on three concentric circles, the lower notes being on the outer circle. Previous studies have shown that the vibration remains confined around the solicitation area [START_REF] Rossing | Vibrational mode shapes in Caribbean steelpans. I. Tenor and double second[END_REF]. Hence, the scan is more particularly focused on G3 (of fundamental frequency f 1 ) and its harmonically tuned neighbours G4 (2f 1 ) and G5 (4f 1 ). Fig. 3 shows the transfer function, measured on the excitation point, in the frequency range [0, 1700] Hz, and the associated mode shapes of the structure. One can see that the first three modes are perfectly tuned like f 1 , 2f 1 and 4f 1 , while the fourth and the fifth departs a little from the perfect harmonic relationship, and are slighlty shifted from the exact 6f 1 and 8f 1 relation. More precisely, around f 1 , the modal shape of the structure is focused on the excited note only. Around 2f 1 , a double peak is clearly visible indicating that the mode is degenerate with two mode shapes around the same frequency. The first one has for eigenfrequency f 2 = 390 Hz and is composed of the second vibration pattern of the G3 note together with the fundamental vibration mode of the G4 note. The second one has its eigenfrequency at f 3 = 397.8 Hz and its mode shape is similar except the fact that the pattern on the G4 note is out of phase. Finally the measurement reveals also that at 4f 1 , two degenerate modes are also at hand, with eigenfrequencies f 4 = 789.5 Hz and f 5 = 799.3 Hz, and companion mode shapes. This modal analysis reveals that in steelpan vibrations, modes appear by parts from the second partial. This is a consequence of the fact that vibrations are strongly localized into notes areas, and is a classical feature in mode localization that has been observed for instance in simple beam systems [START_REF] Luongo | Mode Localization in Dynamics and Buckling of linear Imperfect Continuous Structures[END_REF].

Forced vibrations at higher force amplitudes will now be detailed to depict how energy is transferred between these modes.

FORCED OSCILLATIONS

The previous linear analysis shows that the first frequencies of the steelpan note are organized as a 1:2:2:4:4 sequency. In a nonlinear regime, and in the case of thin curved structures, harmonic relationships create strong modal interactions via internal resonances [START_REF] Nayfeh | Nonlinear interactions[END_REF]. In order to investigate the nonlinear behaviour of the steelpan, the same home-made coil/magnet exciter is used, but higher force amplitudes are applied. The coil/magnet exciter has been thoroughly analyzed for calibration in [START_REF] Thomas | Asymmetric non-linear forced vibrations of free-edge cicular plates. Part 2: experiments[END_REF], where it has been shown that for the vibration amplitudes encountered, harmonic distortion is less than 1%, ensuring a clean and reproducible harmonic excitation. An external sinusoidal current is applied with an external excitation frequency around the first eigenfrequency (f dr f 1 ). The vibratory response of one point of the excited note is measured. The first three harmonics (at f 1 , 2f 1 and 4f 1 ) of the response are precisely measured. Then, two analytical models of internal resonance relationships are proposed to fit the experiment and identify some modal interactions in steelpans vibrations. More details can be found in [START_REF] Monteil | Complicated dynamics exhibited by thin shells displaying numerous internal resonances : application to the steelpan[END_REF][START_REF] Monteil | Vibrations non linéaires de steelpans : couplages modaux via la résonance interne 1:2:2[END_REF].

Experiments

In Fig. 4 the nonlinear response of the steelpan excited around 197.5 Hz with I = 5A is shown. The harmonics 1, 2 and 4 of the recorded displacement are shown, they are denoted respectively by w H1 , w H2 and w H4 . A strong coupling between these harmonics, oscillating at f 1 , f 2 2f 1 and f 3 4f 1 , is revealed. Markers are inserted into the figures to precisely locate, in frequency, the different eigenfrequencies of the system. First, a 1:2:2 internal resonance is observed through the complex behaviour of the measured amplitude at twice the excitation frequency. The energy transfered gives a large amplitude of vibration to the second (390 Hz) and the third eigenmode (789.5 Hz). Secondly, coupling with the fourth harmonic is also observed where w H4 seems to be slaved to w H2 due to the similar global shapes of the two corresponding amplitudes. One can assume that a 1:2:2:4 resonance is here activated. A strong peak of amplitude exists around f dr 192 Hz. At this frequency, the first mode, directly excited, reaches 0.06mm, the second harmonic 0.08mm and the fourth one 0.003mm. Considering the thickness of a typical steelpan initial barrel around 1 mm, one can conclude that geometric nonlinear effects are exhibited for vibration amplitudes of 1/16 times the thickness. Then, a quasiperiodic regim, around 193 Hz, reveals a complex dynamics of the system. This implies that there is a difference between the forward and the backward excitation. A jump phenomenom is observed at 192.6 Hz. Compared to the theory [START_REF] Monteil | Vibrations non linéaires de steelpans : couplages modaux via la résonance interne 1:2:2[END_REF][START_REF] Nayfeh | Nonlinear Oscillations[END_REF], the shapes of the solutions are slightly different than a 1:2 or a 1:2:2 internal resonance, probably because the cubic nonlinearity is not anymore negligible. An other remark on the effects of cubic nonlinearities can be the shift of the curves to the left part of the graph, compared to the linear initial frequencies. The maximum of each harmonic response is observed around f dr 192 Hz.

3-dofs internal resonance models fitting

The complicated dynamics exhibited by the forced vibrations is now modeled by two simple systems involving three internally resonant modes, so as to highlight the most salient features of the dynamics of the steelpan. The two models involve either a 1:2:2, or a 1:2:4 eigenfrequency relationship. For these 3-dofs models, analytical solutions are accessible via multiple scales analysis [START_REF] Monteil | Vibrations non linéaires de steelpans : couplages modaux via la résonance interne 1:2:2[END_REF]. Model fitting to experimental measurements will thus shed new light on the identification of nonlinear coupling coefficients as well as energy transfers.

Considering the steelpan as a curved thin structure, the transverse displacement w can be discretised by expanding w on the linear modes basis {Φ k } as

w(x, t) = N k=1 Φ(x)q k (t), (1)
where q k is the modal coordinate and Φ k is the mode shape amplitude value at the spatial point x.

Now two 3-dofs simple models are proposed, considering the most important coupling found in the modal analysis, and thus displaying respectively a 1:2:2 and the 1:2:4 internal resonance. The 1:2:2 internal resonance considers three modes which eigenfrequencies are tuned such that ω 3 ω 2 2ω 1 . The associated normal form of the amplitude of the three corresponding modes reads :

     q1 + ω 2 1 q 1 = ε [-2µ 1 q1 -α 1 q 1 q 2 -α 2 q 1 q 3 + F 1 cos Ωt] , q2 + ω 2 2 q 2 = ε -2µ 2 q2 -α 3 q 2 1 , q3 + ω 2 3 q 2 = ε -2µ 3 q3 -α 4 q 2 1 , (2)
where ω k = 2πf k denotes the angular frequency of the mode k and µ k its damping coefficient. Ω = 2πf dr is the external frequency and F 1 the amplitude of the external force. In first approximation, only quadratic nonlinear coupling terms α i are kept. Only four of them are present. They correspond to the resonant terms. All other possible nonlinear quadratic terms have no importance for the global dynamics and can be cancelled. Nonlinear coefficients, damping terms and external forcing are assumed to be small as compared to the linear oscillatory part, and thus are scaled by 1. According to the multiple scales development, the solutions of the dynamical system Eqs. (2) are obtained as:

     q 1 = a 1 cos (Ωt + ϕ 1 ) , q 2 = a 2 cos (2Ωt + ϕ 2 ) , q 3 = a 3 cos (2Ωt + ϕ 3 ) , (3)
where a k and ϕ k are amplitudes and phases of the solution q k , respectively. In term of transverse displacement, the equation (1) leads to w(x, t) = w 1 cos(Ωt + ϕ 1 ) + w 2 cos(2Ωt + γ 2 ), where w 1 = Φ 1 a 1 and where w 2 and γ 2 are combinations of Φ 2 , Φ 3 , a 2 , a 3 , ϕ 2 and ϕ 3 .

The 1:2:4 internal resonance, is a 3-dofs nonlinear dynamical system constructed with the following eigenfrequency relationship: ω 3 2ω 2 4ω 1 . It reads:

     q1 + ω 2 1 q 1 = ε [-2µ 1 q1 -α 5 q 1 q 2 + F 1 cos Ωt] , q2 + ω 2 2 q 2 = ε -2µ 2 q2 -α 6 q 2 1 -α 7 q 2 q 3 , q3 + ω 2 3 q 3 = ε -2µ 3 q3 -α 8 q 2 2 , (4
) where α 5,6,7,8 are the four new nonlinear coupling coefficients associated to the 1:2:4 internal resonance system. Solutions of Eqs. ( 4) are:

     q 1 = a 1 cos (Ωt + ϕ 1 ) , q 2 = a 2 cos (2Ωt + ϕ 2 ) , q 3 = a 3 cos (4Ωt + ϕ 3 ) . (5)
In that case, two upper-harmonics, oscillating at 2Ω and 4Ω, simultaneously appear. The transverse displacement is w(x, t) = w 1 cos(Ωt + ϕ 1 ) + w 2 cos(2Ωt + φ 2 ) + w 3 cos(4Ωt + φ 3 ), where w 1 = Φ 1 a 1 , w 2 = Φ 2 a 2 and w 3 = Φ 3 a 3 . Fig. 5 presents the 1:2:2 model fitting (in red color), first compared to the experiment, then compared to a 1:2:4 internal resonance model fitting. Experimental fitting, presented in Fig. 5, determines the value of each nonlinear coefficient while the linear parameters (ω k , µ k ) are identifyed by modal analysis. The amplitude of the external current is I = 2A. The two 3-dofs models allow to locally identify some modal interactions in steelpan vibrations. Each of them locally fits parts of the experimental curves. Markers of linear frequencies indicate the mode activated by nonlinear coupling. This fitting shows that simple models can be used to enhance the comprehension of the complicated dynamics experimentally observed. More complicated scenarii invoking also the presence of the modes at 6f 1 and 8f 1 may also be activated in certain vibratory regimes.

The steelpan is harmonically tuned (both between notes and between harmonics in a note). Therefore, the thinness of the curved structure and the amplitude of sollicitation of the note lead to nonlinear behaviour with modal interactions and energy transfers, consequently. Internal resonance models displaying a 1:2:2 and a 1:2:4 frequency relationships allow to recover the main features of the FRFs,

F k = ω k /2π [Hz] 197.1 389.2 788.6 ξ k = µ k /ω k [%] 0.0045 0.001 0.001 Φ k (x)
0.295 0.339 0.06 while non-modelled effect appears to be easily interpreted. The most complete model for that case should be a 1:2:2:4:4 one, unfortunately analytical solutions for that problem are not tractable.

NONLINEAR MODEL FOR SOUND SYNTHESIS 4.1 Free oscillations measurement

In Fig. 6, two time-frequency representations of free oscillations of the G3 note in normal playing are shown. A weak stroke (see Fig. 6(a)) reveals that most of energy is stored within the first two harmonics. Comparatively, a strong stroke (see Fig. 6(b)) has much more energy on the higher frequencies. It can be noticed that, on these two spectrograms, there is no energy at F = 3f 1 . For the low frequency range, only the tuned harmonics are excited when the steelpan is played. harmonic. First, on the Fig. 7(a), it can be observed that the second harmonic reaches its maximum later than the first one does. Then, the amplitude of the maximum is a little bit higher than the first one. The fourth harmonic presents a very small amplitude. For a strong stroke (see Fig. 7(b)), the amplitude of the second harmonic is much larger than the first one. It can be deduced that energy is transfered from the first mode to the second one, and the intensity of the transfer depends on the intensity of the excitation. The fourth one has much more energy than in the first case. Its amplitude is very close to the first harmonic. Thus, it is shown that initial conditions are very important in the steelpan response.

Free oscillations simulation

The 1:2:4 adjusted model is now analyzed in free oscillations behaviour. The ordinary differential equation solver (ODE45) of Matlab is used. The linear parameters and the nonlinear coupling coefficients values are deduced from the experimental fitting in forced oscillations (see Fig. 5). An initial condition in displacement is given so as to mimic a delta dirac temporal excitation. Fig. 8 shows the result of two simulations. The first one (Fig. 8(a)) is for a weak initial condition (w 0 = 0.15 mm) and the second one (Fig. 8(b)) is for a larger one (w 0 = 0.65 mm). The acceleration of the three transverse displacements (w 1 , w 2 and w 3 ) are represented versus time. Compared to the filtered measurement (Fig. 7), the global time evolution of a steelpan sound is recovered for the two initial conditions. The right amplitude ratios between modes are observed. Also, one can remark that the period of the oscillations is qualitatively the same in the experi- These time simulations confirm that with regard to the global dynamics carried by the envelope modulation of the main harmonics, the 3-dofs model allows to recover the most salient nonlinear features, and are thus identified as a key component in steelpans's sound and vibrations.

CONCLUSION

Nonlinear vibrations of a steelpan presenting a 1:2:4 modal tuning (f , 2f , 4f ) have been investigated thanks to a refined modal analysis, measurements of frequency-response curves in forced oscillations and time domain simulations in free, impacted vibrations, using 3-dofs models identified from the forced response.

The main outcomes of the present study are the following • The modal analysis has clearly evidenced the fact that modes appear by pairs from the second harmonic of each note, a feature that is only scarcely mentioned in the literature. This pairing is interpreted as a consequence of the strong localization of the vibrations into the notes area. Indeed, simpler systems consisting of beams with stiffness imperfections reveals the same degeneracy which is thus generic for structural systems having at least two minima of stiffness [START_REF] Luongo | Mode Localization in Dynamics and Buckling of linear Imperfect Continuous Structures[END_REF].

• A consequence of this mode pairing is that the dynamical equations are complexified. Analytically, one can also show that the appearance of 1:2:2 (instead of simpler 1:2) or e.g. 1:2:2:4:4 (instead of simpler 1:2:4) internal resonances in a nonlinear systems favour instabilities [START_REF] Monteil | Vibrations non linéaires de steelpans : couplages modaux via la résonance interne 1:2:2[END_REF].

• The forced response analysis reveals the complex dynamics of steelpans with appearance of energy exchange and quasiperiodic regimes for very small amplitudes of vibration.

• 3-dofs models that are analytically solvable have been fitted to measured resonance curves, showing that the main features are recovered by considering 1:2:2 and 1:2:4 resonances. However small details in the resonance curves are not fitted, advocating for the fact that the precise dynamics is more complex with energy exchange between pairs of modes and thus solutions comprising e.g. 1:2:2:4:4 resonances and even for some amplitudes 1:2:2:4:4:6 and 1:2:2:4:4:8.

• The temporal simulations allow recovering the main features of the energy exchange, once again confirming that the most important part of the dynamics is carried by the identified resonances. However, small details are still missing to the ear for a sound synthesis approach, stating clearly that the higher modes (at 6f and 8f) have to be included in a simple modal model for sound synthesis.

At the conference, sound synthesis will be presented in order to more precisely quantify the importance of the higher modes for the sound production, as well as the impact model.

  Copyright: c 2013 Mélodie Monteil et al. This is an open-access article distributed under the terms of the Creative Commons Attribution 3.0 Unported License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Figure 1 .

 1 Figure 1. Pictures of steelpans : (a) Steelpans of a typical steelband. (b) Steelpan in normal playing.

Figure 2 .

 2 Figure 2. (a) Steelpan tuning. (b) Modal analysis of the steelpan used for the experiment (right barrel of a double second) excited on G3.

Figure 3 .

 3 Figure 3. Modal analysis of the steelpan excited on the note G3. FRF measurement and associated mode shapes.

Figure 4 .

 4 Figure 4. Frequency response curves I = 5A, excitation frequency f dr in the vicinity of the first mode. Perlines: harmonics 1, 2 and 4 of the measured displacement w. Forward (black line: ) and backward (grey line: ) frequency sweeps. Linear eigenfrequency markers f 1 (--) ; 2f 1 (-• -) and 4f 1 (• • • ).

Figure 5 .

 5 Figure 5. Frequency response curves and experimental fitting (I = 2A) (black line). Nonlinear coupling coefficients identified: 1:2:2 internal resonance model (blue line) -α 1 = 0.42, α 2 = 0.46, α 3 = 0.46, α 4 = 0.39 and 1:2:4 internal resonance model (red line) -α 5 = 0.24, α 6 = 0.34, α 7 = 0.88, α 8 = 0.53.

Fig. 7

 7 Fig. 7 represents the measurements of the Fig. 6 filtered on the harmonics 1, 2 and 4 (denoted H 1 , H 2 , H 4 ), showing the evolution with the time of the amplitude of each

Figure 6 .

 6 Figure 6. Experimental spectrograms of two free oscillation measurement with two stick-impact intensities. (a) Weak stroke. (b) Strong stroke.

4 Figure 7 .Figure 8 .

 478 Figure 7. Same measurements than Fig. 6, filtered on harmonics 1, 2 and 4. Evolution of each harmonics with time.


X You

G G Perez-Sanchez

C Gosset

Continuously-tunable, from MHz to GHz range, short pulse train generation and its dynamic properties

Keywords: short pulse train generation, large tunability, SOA nonlinear amplification, duty cycle tailoring

We propose a new method to generate a continuously tunable repetition rate pulse train. It is based on the filtering of the negative nonlinear phase shift induced by the amplification of a rectangular optical signal in a semiconductor optical amplifier (SOA). A constant pulse width of about 40 ps with a huge continuous tunability of the repetition rate, ranging from 1 MHz to 5 GHz, was experimentally demonstrated. The dynamic properties of the proposed versatile pulse source are studied.

INTRODUCTION

In the field of instrumentation and characterization techniques, narrow pulse sources play an important role. For instance, the asynchronous sampling pulsed coherent detection provides a blind and asynchronous characterization method of optical signal and data flows, with THz bandwidth potentiality [START_REF] Westlund | Software-synchronized all-optical sampling for fiber communication systems[END_REF][START_REF] Gallion | Bandwidth and dynamic range of a pulsed local oscillator coherent optical detection. Application to linear optical sampling[END_REF]. In the field of optical communication, it allows extracting the intrinsic statistical properties as constellation, eye diagram, BER and SNR. This technique requires a local oscillator as a train of brief optical pulse (typically in the 1 ps range or below) with repetition rate typically in the 500 MHz -1 GHz range corresponding to a trade-off between phase noise impairments and cost effective solution. On another hand, application such as pump-probe experiment requires repetition rate in the 10 MHz range to allow relaxation of high-speed phenomena in materials.

Passive Mode-locked lasers are interesting candidates to generate short optical pulse train signal. Fiber laser are adapted to the generation of low repetition rate signal, typically bellow 100 MHz, thanks to long cavities. For higher frequency range, semiconductor laser exhibits impressing performances in term of repetition rate, typically above 10 GHz [START_REF] Gosset | Phaseamplitude characterization of a high repetition rate laser quantum dash passively mode-locked laser[END_REF][START_REF] Akrout | Generation of 10 GHz optical pulses with very low timing jitter using one section passively mode locked quantum dash based lasers operating at 1.55 µm[END_REF]. In the range 100 MHz-10 GHz, it doesn't exist solution based on mode-locking process, due to not suited particle lifetime of the active medium. Furthermore, mode-locked lasers, whatever technology used, suffer from a lack of tunability of the repetition rate, and semiconductor based lasers also exhibits high phase noise of spectral components that can appear as a bottleneck for application such as asynchronous coherent sampling.

Cavity-less pulse sources allow the generation of high repetition rate signal with very good phase noise properties. But the proposed techniques don't allow the generation of periodic signal below 1 GHz [START_REF] Wiberg | Cavity-Less 40 GHz Pulse Source Tunable Over 95 nm[END_REF][START_REF] Wiberg | Photonics processor for anlog-to-digital converter using a cavity-less pulse source[END_REF].

In this paper, we present a simple setup for the generation of versatile pulse train based on the nonlinear amplification of a square signal with a semiconductor optical amplifier. It allows a large continuous tunability of the repetition rate, from 1 MHz to 5 GHz, with a constant pulse width around 40 ps. Furthermore, by using a tunable continuous wave (CW) laser, we obtain also a continuous wavelength tunability. This simple setup is a preliminary step toward a sub-ps versatile pulse train source by combining additional compression stages based on fiber nonlinearities. In the next section, we will first describe the principle of operation of the proposed scheme. Then, in section III, the pulse source is characterized by evaluating it over the characteristics of the driven electrical signal, as the repetition rate and the duty cycles, and of the bias point of the Mach-Zehnder Modulator (MZM). In the last section, the experimental results will be discussed before to conclude.

PRINCIPLE OF OPERATION AND EXPERIMENTAL SETUP

The proposed pulse source is based on the nonlinear amplification using a semiconductor optical amplifier (SOA) of a square waveform signal. SOA non-linearities lead to the front edge distortion, such that the output signal can be seen as the superposition of the linearly amplified signal, and a pulse located at the front edge. Due to the non-zero phaseamplitude coupling factor, the SOA induces optical chirp during the transient amplification, leading to a dissymmetric optical spectrum. As the SOA induces a strictly negative nonlinear phase shift during the dynamic process of signal amplification, the part of the spectrum located toward short wavelength is a replica of the input spectrum, while the one located toward long wavelength is related to the transient regime occurring during the amplification of the front edge. Hence, by filtering out the long wavelength spectrum, the transient signal of the induced-distortion signal train is isolated from the rectangular wave.

Figure 1 shows the experimental setup. The intensity of a tunable CW laser is modulated with a MZM by a pulse pattern generator (PPG). The role of the PPG is to easily generate a periodic rectangular waveform driven signal, with tunable repetition rate and duty cycle. An attenuator at the output of the MZM is used to adjust the power of the modulated signal at the input of the SOA, in order to control the nonlinear phase shift induced by the amplification of the SOA. Due to the polarization dependence of the SOA, a polarization controller (PC) is used at the input of the SOA to maximize the amplification, and thus the induced nonlinear phase shift of the amplified signal. The carrier frequency of the tunable laser is set to 1540.162 nm, and injection current of the SOA 200 mA. The centre wavelength of the flat top optical filter is set at 1540.52 nm (bandwidth: 650 pm).

RESULTS AND ANALYSIS

In this section, we investigate the influence of the driven electrical signal and the bias point of the SOA. We will show their strong influence on the amplification dynamic, and thus on the characteristic of the generated pulses after filtering out the single-sideband (SSB) signal corresponding to the transient amplified signal.

Influence of the repetition rate of the driven electrical signal

As an illustration, we show on figure 3 the spectra of a square signal at the output of the SOA for several repetition rates and a 50% duty cyle. The shapes of the long wavelength spectra for 1GHz and 100 MHz repetition rate signal are identical, with a power ratio of 10 dB. The power ratio corresponds to the repetition rate ratio. It means that for both repetition rates, the energy contains in the transient signal is identical. For the 10 GHz repetition rate, the shape of the long wavelength signal is different. Actually, the energy contains in the transient signal is smaller than for 1GHz and 100 MHz repetition rates. This is due to the dynamic of the SOA. For 1 GHz and 100 MHz repetition rate of the input square signal, the time between two square pulses is long enough (respectively 500 ps and 5 ns) to obtain a complete gain recovery. This is not the case for 10 GHz repetition rate of the input signal. As a consequence, the front edge distortion is reduced and thus, for the latter case, the SPM effect is reduced, leading to a narrower spectrum compared to the 1 GHz and 100 MHz repetition rate. Hence, the chirp is reduced in this configuration.

Influence of the duty cycle of the driven electrical signal

For a given repetition rate, the duty cycle of the driven electrical signal plays an important role in the pulse train generation. To control the period of the signal, we set the length of the sequence generated by the PPG to obtain sequence duration equal to the required period. The PPG is clocked with a 10 GHz clock signal. The duty cycle is then fixed with the number of consecutive marks running within the sequence, while other bits are set to zero. The generated electrical signal is then applied to the driver of the MZM. The rise time of the front edge of the optical signal is measured to be 25 ps. The sequence length is set to 16 bits with a bit duration of 100 ps, leading to a repetition rate of the signal of 625 MHz. The influence of the duty cycles on the pulse generation process is investigated from 1/16 to 15/16 for 8 values totally.

We show in figure 5 the amplified signal at the output of the SOA. The peak intensity of the transient signal is constant for duty cycle lower than 9/16, measured to be 0.6 mW. The duration of the off part of the input signal (number of consecutive zeros) is large enough to allow the full recovery of the gain. For larger duty cycle, the peak intensity decreases, down to 0.2 mW for 15/16 duty cycle due to the lack of recovery time of the gain. As a consequence, the saturated gain is reduced, and the power of the amplified mark decrease with respect to the duty cycle.

For low duty cycle, below 3/16, we observe that the time required to reach the saturation is not enough. Hence, the energy contains in the transient part of the signal is reduced, and thus the energy in the pulse after filtering. We show on figure 6 the impact of the duty cycle on the optical spectrum. The nonlinear part of the spectra, corresponding to the long wavelength, is almost similar for duty cycle between 3/16 and 9/16. As told above, it indicates that the dynamics of the SOA is fast enough to allow a complete gain recovery before the next front edge amplification. For duty cycle from 11/16 to 15/16, we measure a decreasing of the optical power related to long wavelength spectrum, telling us the lack of gain recovery time of the SOA. We observe also that the spectrum related to the 1/16 duty cycle is reduced due to the too short duration of the input signal. Finally, thanks to Parseval theorem, we observe on figure 7 that the energy contains in the pulse after filtering is reduced for 1/16 and 11/16 to 15/16 duty cycles. To confirm the dynamic behaviour of the proposed pulse source illustrated on figures 5 to 7, we measure the mean power of the signal after filtering for several duty cycle, shown on Figure 8. Measurements were done with 100 MHz repetition rate in order to obtain a better resolution. We observe the constant behaviour for intermediate duty cycle range, between 6% and 90%. For duty cycle upper than 90%, the dynamics of the SOA limits the pulse conversion efficiency. For duty cycle lower than 6%, the efficiency is limited by the duration of the input signal.

CONCLUSION

We proposed a new method to generate optical pulse train at several repetition rates, from 1MHz to 5 GHz, with constant pulse width measured to be 40 ps. The dynamical properties of the pulse generator have been studied. They indicate the potentiality of the method to generate pulse train with very low repetition rate, up to a few GHz. The repetition rate upper limitation is induced by the gain recovery time of the SOA, while no limitation exists for very low repetition rate. It is even possible to generate a single pulse with this technique. The use of faster SOA based on quantum dots materials could be interesting to increase further the repetition rate. Nevertheless, the SOA needs to have a highly nonlinear behaviour to not decrease the pulse conversion efficiency. Hence a trade-off between dynamic and conversion efficiency has to be considered for the design of the SOA. This technique is a well promising technique to generate short optical pulses with versatile repetition rate, as it could be associated with compression stage.

Figure 1 .

 1 Figure 1. Experimental setup and the symbolic pulse train at different measurement pointsIn the lower right of figure1, the schematic intensity of the optical waveform is depicted. Figure2shows measured temporal waveforms at point b (fig.2a), and point c (fig.2b) after filtering out the transient regime. The repetition rate is set to 1 GHz with a 50% duty cycle. Figure2c and 2dshow the related optical spectrum.

  Figure 2 shows measured temporal waveforms at point b (fig. 2a), and point c (fig. 2b) after filtering out the transient regime. The repetition rate is set to 1 GHz with a 50% duty cycle. Figure 2c and 2d show the related optical spectrum.

Figure 2 .

 2 Figure 2. a) Temporal waveform of the amplified signal at the output of the SOA. b) Temporal waveform after filtering with a square flat top optical filter. c) Spectrum of a). d) Spectrum of b).

Figure 3 .

 3 Figure 3. Spectra of 100 MHz, 1 GHz and 10 GHz signals amplified by the SOAThe repetition rate tunability is evaluated for 1 MHz, 10 MHz, 100 MHz, 500 MHz and 5 GHz. The different pulse profiles are measured with a 30 GHz bandwidth sampling oscilloscope and depicted in figure4. We observe a constant amplitude for repetition rates up to 1 GHz, corresponding to a constant energy in the transient regime. A constant pulse width of 40 ps is obtained. At 5 GHz, we observe a strong diminution of the amplitude of the front edge distortion, correlated with an increase of the pulse width, due to the impact of the SOA dynamics. Hence, it induces a decrease of the efficiency conversion.

Figure 4 .

 4 Figure 4. Pulse profiles for several repetition rate

Figure 5 .

 5 Figure 5. Pulse profiles for several duty cycles from 1/16 to 15/16

Figure 6 .

 6 Figure 6. Spectra of the signal amplified by SOA measured for several duty cycles

Figure 8 .

 8 Figure 8. Pulse power (a.u.) for several duty cycles

ACKOWLEDGEMENT

This work is supported by the "OCELOT" project, from the French National Research Agency (ANR).


H Rifaï

J.-F Guerrero-Castellanos

N Marchand
email: nicolas.marchand@gipsa-lab.grenoble-inp.fr

G Poulin-Vittrant
email: guylaine.poulin-vittrant@univ-tours.fr

Biomimetic-Based Output Feedback for Attitude Stabilization of a Flapping-wing Micro Aerial Vehicle †

The paper deals with the development of a bounded control law for Flappingwing Micro Aerial Vehicles (FMAVs) that mimics a strategy adopted by animal flapping flyers to stabilize their orientation. The control consists on generating torques about the body's principal axes by means of a modulation of the wing angle amplitudes. It is known that flapping flyers orient their body without any numerical computation or estimation of their current attitude. Therefore, the proposed control law is computed using the direct measurements of on-board sensors mimicking animal sensitive organs, more specifically the halteres, legs sensilla and magnetic sense. The technological equivalents of these biological sensors are respectively three rate gyros, a tri-axis accelerometer and a tri-axis magnetometer. Besides, the control signal is bounded to keep the wing angle amplitudes below the maximal values. Owing to its simplicity, this control law is suitable for applications where on-board computational resources are limited. The stability of the closed loop system is proved based on Lyapunov analysis and averaging theory. The effectiveness of the proposed control law is shown in simulations. The robustness with respect to external disturbances is also shown emphasizing the importance and need of the bounded control.

Introduction

Animal flapping flyers depict several techniques to achieve flight and maneuvers moving their wings, body and legs [START_REF] Dudley | The biomechanics of insect flight: form, function, evolution[END_REF]. Principally, these techniques get benefit of the wing morphology, aerodynamic effects, sensory and actuating systems, flight control and obstacle avoidance mechanisms. For example, insects are able to change, very quickly, their speed and direction of flight, within almost 100 ms especially during predatory phases [START_REF] Dudley | The biomechanics of insect flight: form, function, evolution[END_REF]. They are also able to accomplish a transitory lateral or backwards flight. These maneuvers are performed by a displacement of the abdomen or an asymmetrical evolution of the wings using an amplitude modulation [START_REF] Nachtigall | Neuro-muscular control of dipteran flight[END_REF][START_REF] Alexander | Wind tunnel studies of turns by flying dragonflies[END_REF] or phase modulation [START_REF] Mountcastle | Vortexlet models of flapping flexible wings show tuning for force production and control[END_REF], frequency modulation occurs very rarely [START_REF] Nachtigall | Neuro-muscular control of dipteran flight[END_REF].

Besides, the flapping flyer is endowed with information, issued from multiple sensitive organs, about its state, interaction with the environment as well as local destination to determine the way to reach it. For insects, one can cite the ocelli, the compound eyes and other biological sensors detailed in the following [START_REF] Dudley | The biomechanics of insect flight: form, function, evolution[END_REF][START_REF] Alexander | Nature's flyers: Birds, insects and the biomechanics of flight[END_REF][START_REF] Campolo | Attitude stabilization of a biologically inspired robotic housefly via dynamic multimodal attitude estimation[END_REF]) (and references therein):

• Halteres: are gyroscopic biological sensors, present at the wing bases, that detect the rotational movement of the body and allow to determine its angular velocity along the three axes,

• Sensilla: are cuticular sense hair detecting chemical or mechanical stimuli. They are present on the antenna, wings and legs. Legs sensilla, for example, allow to determine the direction of the gravity field with respect to the insect's body,

• Magnetic sense allows to determine the direction of the earth magnetic field with respect to the insect's body. Using a fusion of information issued from these organs, the insect determines its trajectory and adapts its body's velocity and orientation to track it. One should emphasize that these actions are executed by the insect without any instantaneous numerical knowledge of its position, orientation or velocity.

Insects performance has encouraged the design of Flapping-wing Micro Aerial Vehicles (FMAVs). The micro aerial vehicle is a small size aircraft, having a maximal dimension of 6 inches †, and intended to perform an autonomous flight, thanks to an onboard control system, comprising a set of sensors and a dedicated integrated circuits. Due to their small size and flapping movement, FMAVs fly in zones characterized by low Reynolds numbers (10 2 -10 4 ). They develop, therefore, extra non-stationary aerodynamic forces, which help them reduce the forces generated by the actuators to accomplish the flight and, consequently, reduce the energetic consumption. FMAVs produce low noise. They are able to accomplish vertical taking-off and landing as well as stationary flight in hovering mode. However, their major drawback is still the difficulty of identifying and implementing the complex mechanisms carried out by insects to perform maneuvers [START_REF] Dudley | The biomechanics of insect flight: form, function, evolution[END_REF][START_REF] Hedrick | Flight control in the hawkmoth Manduca sexta: the inverse problem of hovering[END_REF]. The FMAVs are intended for use in areas inaccessible for people or that require high accuracy of intervention.

The present paper deals with the modeling and attitude stabilization of a FMAV for scenarios necessitating to maintain a stable orientation in front of a scene in order to monitor or investigate it. A simple model mimicking the flapping flight is presented, including the wings degrees of freedom, the sensorimotor system as well as the developed aerodynamic forces. The major contribution of this work is the development of a biomimetically inspired control strategy aiming to stabilize the orientation of the FMAV. The inputs of the control law are the direct measurements of on-board sensors, equivalent to those which an insect is equipped with, without the need of an explicit attitude reconstruction (represented by the Euler angles in R 3 , quaternion in S 3 ⊂ R 4 or rotation matrix in SO(3)). Hence, unlike conventional approaches, the computational cost used for the attitude estimation/reconstruction is avoided. Note that attitude estimation is frequently carried out by means of Extended Kalman Filters or Nonlinear observers which represent a high computational cost. The control law takes into consideration the wing angle amplitude bounds, characterizing the animal species. It has the easiness of a PD controller, the derivative term is obtained by means of the halteres (angular velocity) and the proportional term by means of legs sensilla and magnetic sense (direct measurement of attitude error). Moreover, the control is very simple and therefore is suitable to be implemented in real time. It is also independent of the FMAV body's inertia, modeling or aerodynamic errors and robust with respect to external disturbances. Unlike fuzzy controllers or fractional PID, the bounds of the control signal are taken into account explicitly and the stability properties are well established.

Models of FMAVs existing in literature consider them as rigid bodies subject to external forces and torques generated by the flapping wings whose amplitudes may be modulated (Deng et al. 2006a;[START_REF] Rakotomamonjy | Longitudinal modelling and control of a flapping-wing micro aerial vehicle[END_REF], phase or even both of them [START_REF] Chung | Neurobiologically inspired control of engineered flapping flight[END_REF][START_REF] Oppenheimer | Dynamics and control of a biomimetic vehicle using biased wingbeat forcing functions: Part I-aerodynamic model[END_REF][START_REF] Couceiro | Modeling and control of biologically inspired flying robots[END_REF]. Few works have treated the problem of controlling the orientation of flapping aerial vehicles based on sensor measurements. A proportional derivative output feedback is proposed in [START_REF] Schenato | Attitude control for a micromechanical flying insect via sensor output feedback[END_REF]) based on ocelli and halteres measurements. A linear quadratic optimal control is proposed in [START_REF] Deng | Model identification and attitude control for a micromechanical flying insect including thorax and sensor models[END_REF](Deng et al. , 2006b) ) based on ocelli, magnetic compass and halteres measurements. In [START_REF] Epstein | Biologically inspired feedback design for drosophila flight[END_REF], halteres measurements are used to estimate the pitch rate and stabilize the corresponding angle with a proportional derivative controller with pole placement. Rate gyros, accelerometer and magnetometer sensors, representing respectively the halteres, legs sensilla and magnetic compass, are used in [START_REF] Campolo | Attitude stabilization of a biologically inspired robotic housefly via dynamic multimodal attitude estimation[END_REF] to estimate the insect's attitude (rotation matrix), used after, along with angular velocity measurement, in a state feedback control law. One should emphasize that, in the aforementioned works, the proposed linear control laws cannot be sufficiently robust with respect to external disturbances representing wind for example. Note also that building a control law using estimated attitude state is not a biomimetic approach. In fact, insects have no numerical instantaneous determined orientation knowledge (angles or rotation matrix) but only information issued directly from their biological sensors [START_REF] Chapman | The insects: Structure and function[END_REF].

The paper is organized as follows. In section 2, the model of the FMAV is proposed in agreement with insect models, including the sensorimotor systems, aerodynamic forces and body's dynamics. A biologically inspired bounded control law, based directly on the measurements of some embedded sensors and aiming to stabilize the orientation of the FMAV, is presented in section 3. Simulation results are addressed in section 4 as well as some robustness tests. Finally, conclusions are given in section 5 and future works are introduced.

Insect flight versus biomimetic robot flight

The present section deals with the design and modeling of a flapping-wing micro aerial vehicle. Design concerns mainly the choice of the sensors and actuators to embark on the FMAV. Therefore, the insect's sensorimotor system is shortly presented and correlated to the technological corresponding equipments. Modeling concerns principally the establishment of a simple mathematical model allowing to represent the degrees of freedom of the flapping flight as well as aerodynamic forces and mechanisms deployed by nature's flapping flyers to perform maneuvers. One should note that a beating wing is in interaction with the surrounding air flow, which creates aerodynamic forces perpendicular to its surface. These forces generate the linear and rotational movements of the FMAV's body.

Wing degrees of freedom

The movement of a wing is a combination of several elementary actions [START_REF] Thomson | Experiment-based optimization of flapping wing kinematics[END_REF][START_REF] Chapman | The insects: Structure and function[END_REF]): flapping, feathering, lagging or elevation besides flexion and torsion (Figure 1 and 2). Flapping is an up-and-down movement of the wing represented by a rotation of the wing about its tangential axis t of a flapping angle φ. Feathering is a rotation of the wing about its span-wise axis r of a rotation angle ψ. Lagging or elevation is a forward-backward movement of the wing parallel to the body, modeled by a rotation about a normal axis n of a deviation or elevation angle θ. Flexibility of the wing allows it to resist to turbulence, provides a gentler flight and increases the aerodynamic force relative to a same size rigid wing [START_REF] Shyy | Recent progress in flapping wing aerodynamics and aeroelasticity[END_REF]. Torsion is a twist movement of the wing, providing an aerodynamic stability [START_REF] Senda | Effects of flexibly torsional wings in flapping-of-wings flight of butterfly[END_REF]. The directions of the wing rotational axes are chosen such that r is oriented from the wing base to its tip along the wingspan, t is parallel to the wing chord, oriented from trailing to leading edge and n is perpendicular to the wing plane oriented so that the three-sided frame F w ( r, t, n) is direct. Note that the frame F w should be indexed, left F w l and right F w r , relative to the left and right wings respectively. Maneuvers of animal flapping flyers are performed primarily by asymmetrical movements of the wings, a displacement of the abdomen and legs. The movement of the wings generates aerodynamic forces, which resultant is perpendicular to the wing surface. Therefore, the flapping/rotation movements of the wings generate the vertical/longitudinal movements of the insect, respectively. A difference in the amplitudes of the left and right flapping/rotation angles allows to generate a roll/yaw movements, respectively (Figure 2). The pitch movement can be created by controlling the wing's elevation degree of freedom, or by changing the center of gravity of the body tilting it upward or downward. The second solution allows to reduce the number of embarked actuators and therefore is adopted in the sequel. Finally, the lateral movement is generated by coupling the body vertical lift force with the roll angle, maneuver performed by many birds and insect species [START_REF] Dudley | The biomechanics of insect flight: form, function, evolution[END_REF].

The wing angles are characterized by their maximal amplitudes and their wingbeat frequency within their predefined trajectory. These values are specific to each species. Generally, flapping-wing flying creatures move their wings according to a sine function with higher harmonics, the fundamental frequency is equal to the wingbeat frequency [START_REF] Dudley | The biomechanics of insect flight: form, function, evolution[END_REF]. A wingbeat period is divided into two phases: the downward phase of the flapping movement or downstroke and, the upward phase or upstroke [START_REF] Chapman | The insects: Structure and function[END_REF][START_REF] Dudley | The biomechanics of insect flight: form, function, evolution[END_REF]. During downstroke, the flapping flyer orients the dorsal side of the wing to the flow. It orients the ventral side during the upstroke. The model proposed in this work considers the wing as a rigid body, flapping in the mean stroke plane, defined by taking the deviation angle θ to zero. This choice is motivated by the possibility to control insects using only two wing angles [START_REF] Lentink | Nature-inspired flight -beyond the leap[END_REF]. Flapping and rotation angles, φ and ψ, are assumed to vary according to saw tooth and pulse functions respectively, such that the wing changes its orientation at the end of each half stroke (Figure 3). The time variation of the wing angles is given by (2.1):

φ(t) = φ 0 (1 -2t κT ) 0 ≤ t ≤ κT φ 0 (2 t-κT (1-κ)T -1) κT < t ≤ T ψ(t) = ψ 0 sign(κT -t) 0 ≤ t ≤ T θ(t) = 0 0 ≤ t ≤ T (2.1)
where sign designates the classical sign function, T is the wingbeat period, κ is the ratio of downstroke duration to the wingbeat period and should verify 0 < κ < 0.5 to produce a positive aerodynamic lift force over a wingbeat period, φ 0 and ψ 0 are respectively the amplitudes of flapping and rotation angles. φ 0 and ψ 0 , considered for left and right wings, will define the control inputs in the following. This choice is motivated by the fact that wingbeats have generally fixed frequency, except during some maneuver phases. Turns can be created by asymmetrically changing the left and right wing angle amplitudes [START_REF] Nachtigall | Neuro-muscular control of dipteran flight[END_REF][START_REF] Alexander | Wind tunnel studies of turns by flying dragonflies[END_REF]). Note that the wingbeat frequency has been taken equal to 100 Hz considering in the following the model of a diptera insect [START_REF] Dudley | The biomechanics of insect flight: form, function, evolution[END_REF]. Note also that the given angle parametrization does not represent the real movement of the wings, but the desired trajectory to be achieved by the wing actuators.

Remark 1. The wing angle parametrization adopted in the present work is not unique. Other parameterizations can be used as well. The only condition is that they allow the generation of an aerodynamic body lift force that, averaged over a wingbeat period, should be positive and able to balance the FMAV's weight. The control strategy and control law proposed in the following remain always applicable.

Other non-symmetric wing-parametrizations have been proposed in the literature. In [START_REF] Chung | Neurobiologically inspired control of engineered flapping flight[END_REF], all three angles have been modulated. However, embarking many actuators is weight costly. Therefore, strategies have been oriented to decrease their number. In [START_REF] Oppenheimer | Dynamics and control of a minimally actuated biomimetic vehicle: Part II-control[END_REF], two actuators controlling the wings have been used besides a bobweight actuator in the body. In [START_REF] Oppenheimer | Dynamics and control of a biomimetic vehicle using biased wingbeat forcing functions: Part I-aerodynamic model[END_REF], the bobweight has been removed and the wings parametrization adapted to use only one actuator per wing. This configuration necessitates the development of complex control laws to stabilize the translational and rotational movement of the FMAV. In [START_REF] Finio | Asymmetric flapping for a robotic fly using a hybrid power-control actuator[END_REF]), three actuators have been introduced, one principal and two secondary to accomplish an asymmetric control of the wing angle amplitudes. Another work has considered one active and one passive degrees of freedom allowing to control four wing angles, the aim is to ensure the vertical flight, the orientation, notably the roll and yaw angles, is passively regulated [START_REF] Sreetharan | Passive torque regulation in an underactuated flapping wing robotic insect[END_REF].

Actuation

The periodic wingbeats are accomplished by the muscles, present at the wing bases or in the thorax and, vibrating with a fixed frequency. This movement can be identified as a vibration at the undamped natural frequency. Therefore, piezoelectric actuators are the most suitable for this application. The reverse effect of piezoelectricity, consisting on applying voltage and retrieving mechanical movement, is of interest. The alternative voltage is delivered by an electronic converter designed specifically for piezoelectric actuators. These actuators behave as reactive loads [START_REF] Janocha | New approach to a switching amplifier for piezoelectric actuators[END_REF][START_REF] Campolo | Efficient charge recovery method for driving piezoelectric actuators with quasi-square waves[END_REF] introducing some non-linearities (hysteresis, creep) that can be avoided using a local control [START_REF] Kuhnen | A new drive concept for high-speed positioning of piezoelectric actuators[END_REF]. The controller inputs are the reference angles defined in (2.1). The local closed loop of controller-actuator is regulated to behave as a first-order filter having fast dynamics so that the influence of the actuator on the global system dynamics is despised (Figure 3):

Ä = Är -λ 1 ( Ȧ -Ȧr ) -λ 2 (A -A r ) (2.2)
with A the amplitude of the flapping or rotation angles at the actuator's output and A r is the reference amplitude at the actuator's input. λ 1 and λ 2 are computed using pole placement such that the time constant of the local closed loop is equal to τ = 0.1 ms: λ 1 = 2 τ and λ 2 = 1 τ 2 . Remark 2. Note that the piezoelectric actuators available on the market have very fast dynamics with a time response reaching the microsecond range (PZT piezoceramic actuators for example). The conservative hypothesis considered for τ = 0.1 ms is therefore absolutely realistic even when the actuator is coupled to the wing. Note that if the time response of the wing-actuator system is less than 0.1 ms, it will be more advantageous for the transparency of the actuator and its influence on the FMAV's movement.

Sensory system

In attitude stabilization theory, the FMAV's body should reach a desired orientation and maintain it, which means that, at stationary hovering flight, the body axes should be aligned to some reference axes with a null angular velocity. Within insect bio-sensors, the halteres, legs sensilla and magnetic sense contribute to control the body's orientation. These sensitive organs have some technological equivalents: respectively, the rate gyros, accelerometers and magnetometers.

Define a fixed frame in the space F f ( x f , y f , z f ) and a mobile frame attached to the FMAV's body at its center of gravity F b ( x b , y b , z b ) (indexes f and b stand for fixed and body, respectively) (Figure 2). Define also a rotation matrix R ∈ SO(3) = {R ∈ R 3×3 : R T R = I 3 , det R = 1} allowing the transformation from the fixed frame F f to the body frame F b .

Rate gyros

Three rate gyros are mounted orthogonally on the FMAV's body such that their axes coincide with the body's axes. The sensors deliver the angular velocity measurements about the body's axes. The measured angular velocity, ω G , is given by:

ω G = ω + ξ G (2.3)
where ω ∈ R 3 is the angular velocity of the body, ξ G ∈ R 3 is a white gaussian noise of the rate gyros.

Magnetometer

A tri-axis magnetometer is mounted on the FMAV such that its axes coincide with the body's axes. The sensor gives the measurement of the fixed magnetic field in F b :

s b M = R s f M + ξ M (2.4)
where s f M ∈ R 3 is the magnetic field in the fixed frame F f , R is the rotation matrix from the fixed frame F f to the body's frame F b and ξ M ∈ R 3 is a vector of gaussian white noise. Note that the magnetic field in the fixed frame is considered constant and equal to

s f M = [ 1 2 , 0, - √ 3
2 ] T in the geographic zone where the simulations are performed.

Accelerometer

A tri-axis accelerometer is mounted on the FMAV such that its axes coincide with the body's axes. The sensor measures the FMAV's acceleration in F b . It is given by:

s b A = R s f A + ξ A
where s f A = (a -ge 3 ) with e 3 = [0, 0, 1] T , g ∈ R and a ∈ R 3 are respectively the gravity and body's acceleration vectors expressed in the fixed frame F f . R ∈ R 3×3 is the rotation matrix from the fixed frame F f to the body frame F b . ξ A ∈ R 3 is a gaussian white noise. Considering that the FMAV moves at low accelerations such that a g and normalizing, the accelerometer measurement is given by:

s b A = R g n + ξ A (2.5)
with g n = s f A = [0, 0, -1] T the normalized gravity vector.

Remark 3. The FMAV experiments acceleration along the body-fixed direction z b for the taking-off, landing and translational horizontal movements. Because in the present work, only the attitude stabilization is considered, the assumption a g holds and the tri-axis accelerometer can be used as a reference vector sensor. The gravity vector is the reference vector and the acceleration is considered as a disturbance. If the translational movement has to be considered, a compensation of the vertical acceleration should be envisaged.

Aerodynamics

The flapping movement of the wings within a surrounding steady air flow generates quasisteady aerodynamic forces: a lift force, perpendicular to the wing, enhancing the flight and a drag force, parallel and opposing to the air flow [START_REF] Dudley | The biomechanics of insect flight: form, function, evolution[END_REF]. Drag forces generated by the wings are neglected in the present work considering that the wings are made of materials having a sufficiently small friction coefficient and that the drag is generated only by the body [START_REF] Dudley | The biomechanics of insect flight: form, function, evolution[END_REF]. In addition to the quasi-steady aerodynamic forces, the flapping flyer is subject to unsteady aerodynamic forces generated by the rotation of the wing about its radial axis r. This rotation creates vortices, adding air flow to the following wingbeats, and creating additional forces. These forces allow the flapping-wing animal to ensure maneuvers such as quick variations of velocity or flight direction and instantaneous turns. All forces are considered applied at the wing's center of pressure, located at l r equal to 65% of the wing length L from the wing base (l r = 0.65L) and at l t equal to 25% of the chord length C h from the leading edge (l t = 0.25C h ) [START_REF] Schenato | Controllability issues in flapping flight for biomimetic micro aerial vehicles (MAVs)[END_REF]. The center of pressure is considered belonging to the radial axis r, given in the wing frame F w by p w = [l r , 0, 0] T .

Quasi-steady aerodynamic force

The quasi-steady force is generated by the pressure of the air flow exerted on the wing surface. This force is perpendicular to the wing and is applied at the wing's center of pressure. It is oriented to the opposite direction of the wing's velocity. The magnitude of this force is given by:

f qs = - 1 2 ρC w S w v w |v w | (2.6)
with ρ the air density, S w the wing surface and v w the wing velocity. C w is the aerodynamic coefficient of the wing.

C w = C(1 + C f ) during downstroke and C w = C(1 -C f )
during upstroke, with C ≈ 3.5 is the force coefficient derived empirically in [START_REF] Dickinson | Wing rotation and the aerodynamic basis of insect flight[END_REF][START_REF] Schenato | Controllability issues in flapping flight for biomimetic micro aerial vehicles (MAVs)[END_REF]) and C f is a coefficient chosen so that the aerodynamic force is 20% greater during downstroke relative to upstroke. This dissymmetry emphasizes the fact that the convex dorsal side of the wing is oriented to the flow during the downstroke, while the concave ventral side of the wing is opposed to the flow during upstroke. The wing camber alteration is due to the stroke reversal of the air circulation around the wing, reducing the effective area of the wing [START_REF] Dudley | The biomechanics of insect flight: form, function, evolution[END_REF]. Therefore, the downstroke force is presumably higher than the upstroke one.

Rotational force

The rotation of the wing about its radial axis deviates the surrounding air flow. The wing reacts to this phenomenon by creating additional rotational circulation [START_REF] Sane | Review The aerodynamics of insect flight[END_REF] and consequently a rotational force modeled by [START_REF] Rakotomamonjy | Longitudinal modelling and control of a flapping-wing micro aerial vehicle[END_REF]:

f r = πρl r C 2 h ( 3 4 - l t C h )v w ψ (2.7)
with ψ the first derivative of the rotation angle.

Added mass force

The added mass phenomenon is created by the acceleration of the additional fluid mass surrounding the wing when it accelerates and rotates. It can be modeled by [START_REF] Rakotomamonjy | Longitudinal modelling and control of a flapping-wing micro aerial vehicle[END_REF]:

f m = π 4 ρLl r C 2 h φ (2.8)
with φ the second derivative of the flapping angle.

In addition to these forces, the wing is subject to other phenomena like the wake capture, delayed stall, wagner effect, etc. that have a minor contribution to the total wing lift force and are difficult to model [START_REF] Sane | Review The aerodynamics of insect flight[END_REF].

The aerodynamic force generated by a wing is applied at its center of pressure, has the direction of the normal vector of the wing n. Its expression in the wing's frame F w is given by:

f w = f qs + f r + f m (2.9)
The total aerodynamic force f b ∈ R 3 , generated by the left and right wings, expressed in the body frame F b is given by

f b = R b l f w l + R b r f w r (2.10)
with R b l,r are the rotation matrices from the left or right wing frames to the body frame. The aerodynamic torque expressed in the body frame is given by

τ b = p b l × f b l + p b r × f b r (2.

Body's dynamics

The body can be considered as a whole entity and is modeled as a rigid body, to which are attached two wings. The wing inertia is neglected in the present work because their mass is less than 5% of the body's mass [START_REF] Schenato | Controllability issues in flapping flight for biomimetic micro aerial vehicles (MAVs)[END_REF]. The effect of the wings' inertia due to the high flapping frequency is considered beyond the scope of this paper. Therefore, the FMAV is modeled as a rigid body subject to aerodynamic forces and torques that generate its movement. It is subject also to viscous and gravitational forces.

The rotational kinematics and dynamics of the FMAV are given by (2.12,2.13). The translational dynamics are not presented here for sake of simplicity. The aerodynamic lift force balances the gravity effect in order to stabilize the FMAV in hovering mode. For more details, readers can refer to previous works [START_REF] Rifai | Bounded control of a flapping wing micro drone in three dimensions[END_REF].

Ṙ = Rω × (2.12) ω = J -1 (τ b -ω × Jω) (2.13)
ω ∈ R 3 is the angular velocity with respect to the mobile frame F b . τ b ∈ R 3 is the aerodynamic torque vector defined in F b . J ∈ R 3×3 is the inertia matrix of the body relative to F b . R is the rotation matrix from the fixed frame F f to the mobile frame F b . ω × is the skew symmetric matrix associated to the vector ω and related to the cross product × such that ω × a = ω × a, with a ∈ R 3 .

2.6. Control strategy Aerodynamic forces are function of the wing angles, their derivatives and some geometric and aerodynamic parameters. Controlling the orientation of the FMAV amounts from controlling the amplitude or the frequency of the wing angles. Considering that the nature's flapping flyers use generally a fixed characteristic frequency except during instantaneous maneuvers execution [START_REF] Dudley | The biomechanics of insect flight: form, function, evolution[END_REF][START_REF] Nachtigall | Neuro-muscular control of dipteran flight[END_REF], the wing angle amplitude control is considered.

On the other hand, the wingbeat frequency adopted in the present work is relatively high (100 Hz). Aerodynamic forces and torques, which are generated at the wingbeat frequency, affect the body's movement only by their average values, computed over a wingbeat period. This is proved based on the averaging theory [START_REF] Khalil | Nonlinear systems[END_REF][START_REF] Bullo | Averaging and vibrational control of mechanical systems[END_REF][START_REF] Vela | Averaging and control of nonlinear systems[END_REF]) which is a strategy usually used in FMAVs control [START_REF] Schenato | Controllability issues in flapping flight for biomimetic micro aerial vehicles (MAVs)[END_REF]. Denote by u = (φ l (t), φ r (t), ψ l (t), ψ r (t)) the flapping and rotation angles for left and right wings, v = (φ l 0 , φ r 0 , ψ l 0 , ψ r 0 ) the amplitudes of these angles (2.1), then u = v f 2 (t). Denote also by x = (R, ω) the FMAV's rotational state, the model given by (2.12-2.13) can be written in a compact form as:

ẋ = f 1 (x, u) (2.14)
Let x = ( R, ω) denote the averaged state over a wingbeat period T . The compact form of the FMAV's model, averaged over a wingbeat period, is given by:

ẋ = f1 (x, v) (2.15)
Averaging theory shows that an exponential stable equilibrium state for the averaged dynamics of a high frequency oscillating system (x = x e ) is also a stable equilibrium state for the oscillating (time variant) system: there exists k > 0 such that x(t) ⊗ x(t) < kT for all t ∈ [0, ∞), with ⊗ designates a comparison operator and T the wingbeat period which is small in the present case.

As mentioned previously, the FMAV is controlled indirectly by means of the wing angle amplitudes v that will be computed using a feedback h(•) of direct sensors measurements ω and s k , k ∈ {1, . . . , n} with n the number of on-board reference sensors:

v = h(ω, s k ),
(2.16)

Using the inverse reasoning, the control strategy is defined as follows. On one hand, based on the desired orientation and the current one measured by on-board sensors, a control torque is computed. This control is actually equal to the average torque τ over a wingbeat period: τ = U(ω, s k ). On the other hand, based on the expression of the aerodynamic torque, a relation between the mean torque τ and the wing angle amplitudes is determined. Recall that the flapping and rotation angle amplitudes (φ 0 , ψ 0 ) are the control inputs. Therefore: τ = Λ(φ 0 , ψ 0 ). The inverse relation can be used to compute the wing angle amplitudes, that should be applied at the beginning of a wingbeat period, function of the control torque:

(φ 0 , ψ 0 ) = Λ -1 (τ ) = Λ -1 (U(ω, s k )). Therefore, h(•) = Λ -1 (U(•)).
Two actuators are used for each wing to generate the flapping and rotation angles, creating then the roll and yaw movements of the FMAV. The pitch movement is generated by moving a small mass inside the body using the ElectroWetting On Dielectric (EWOD) technology [START_REF] Renaudin | Droplet manipulation using SAW actuation for integrated microfluidics[END_REF]. A similar strategy using a bobweight actuator inside the FMAV has also been adopted in [START_REF] Oppenheimer | Dynamics and control of a minimally actuated biomimetic vehicle: Part II-control[END_REF]).

Remark 4. Analyzing only the rotational movement of the FMAV, the system can be considered as fully actuated. However, when studying the translational movement, the system will become underactuated and a coupling of the rotational and translational dynamics should be performed to execute a three dimensional movement. Hence arises the importance of the attitude control.

Computing the averaged roll and yaw torques, function of the flapping and rotation angle amplitudes, one has the explicit form of the trigonometric function Λ(•):

τ1 = βl r φ r 2 0 cos ψ r 0 -φ l 2 0 cos ψ l 0 τ3 = αl r φ r 2 0 sin ψ r 0 -φ l 2 0 sin ψ l 0 (2.17) with α = 2 T 2 1+(1-2κ)C f κ(1-κ) ρCS w l 2 r β = 2 T 2 1-2κ+C f κ(1-κ) ρCS w l 2 r
Each species is characterized by maximum flapping and rotation angles. In a practical point of view, this represents also a technological constraint because the actuators inputs should be bounded to avoid their saturation. Therefore:

0 ≤ φ 0 ≤ φ 0max 0 ≤ ψ 0 ≤ ψ 0max (2.
18)

The roll and yaw torques are therefore bounded. Note also that the movement of the mass inside the body is limited by the body length, the mass and velocity used within the EWOD technology, therefore, the pitch control is bounded too. Resuming, τ ∈ [-τ max , τ max ].

Biologically inspired attitude stabilization

The attitude stabilization problem consists in reaching a desired orientation and maintaining it all over the time. The error between the desired orientation defined by R d and the current orientation defined by R is denoted by R e = RR T d . Therefore, the attitude stability condition is expressed as:

ω → 0 R e → I 3
with I 3 the three dimensional identity matrix. The angular velocity can be directly accessed based on the measurements of the rate gyros. However, no measurement of the rotation matrix can be obtained without using estimation strategies. Therefore, a biologically inspired technique is used in this work. It consists on defining the attitude error by means of the FMAV's on-board sensors which are equivalent to the animal flapping flyers' sensitive organs as it was described earlier. It is actually equal to the error between the on-board sensor measurements and the desired measurements without any computation of the body's orientation. Let s b k and s f k , k ∈ {1, . . . , n}, denote the measurement of the on-board sensor k in the body frame F b and fixed frame F f , respectively, with n the total number of the embarked reference sensors. The projections of the unit vector s k in the body and fixed frames are linked through the rotation matrix R such as s b k = Rs f k . Note that s b k is called vector observation and s f k reference vector. The attitude error γ is defined as:

γ = 1 n n k=1 (Rs f k ) × (R d s f k ) = 1 n n k=1 s b k × R d s f k (3.1)
Reaching the desired attitude is expressed by a rotation matrix R equal to the desired one R d , i.e. relative to the accessible information γ = 0. Reversely, γ = 0 means that s b k and R d s f k are collinear (3.1). Since s k is a unit vector, then

s b k = ±R d s f k . Two cases can be identified: i) s b k = R d s f k i.e
. the vectors observation have the same orientation of the reference vectors desired direction in F b which defines the final goal to reach, and ii)

s b k = -R d s f k i.e
. the vectors observation have the opposite direction of the desired reference vectors direction F b . Therefore, they are defined by a rotation of 180 • (π rad) with respect to a plane, axis or point. The symmetry surface depends actually on the number of reference sensors or more precisely on the number of non-collinear reference vectors n such that k ∈ {1, . . . , n}. Three cases are analyzed: a) n = 1: If only one vector observation is afforded, it is linked to the reference vector by means of the rotation matrix R as s b 1 = Rs f 1 providing a two-dimensional constraint and unveiling any rotation about the axes s b 1 or s f 1 , in the body and fixed frames respectively. Therefore, at least two different sensor measurements should be provided instantaneously to have complete information about the attitude error [START_REF] Markley | Nonlinear attitude filtering methods[END_REF]: n ≥ 2. The case n = 1 is therefore excluded. b) n = 2: The symmetry surface is limited to an axis n perpendicular to the plane containing the two vectors R d s f k , k ∈ {1, 2}, at the vectors intersection. Therefore,

s f k = -R d s f k = R s R d s f k
with R s the rotation matrix defining the symmetry, i.e. the rotation matrix about the axis n of 180 • . In the specific case where R d s f k belong to one of the basis planes (e f 1 , e f 2 ), (e f 1 , e f 3 ) or (e f 2 , e f 3 ), then the rotation is performed respectively with respect to the axes e f 3 , e f 2 and e f 1 . These special cases chiefly facilitate the problem and therefore, a choice of the reference sensors should be done this way. This technique is adopted in the present paper. The matrix R s is given by:

R s ∈ {diag(-1, 1, 1), diag(1, -1, 1), diag(1, 1, -1)} (3.2) c) n ≥ 3: s b k = -R d s f k = SR d s f
k with S the reflexion matrix defined by diag(-1, -1, -1). Note that S is not a rotation matrix because det S = -1 and Tr S = -3. No rotation matrix can therefore be identified and the configuration s b k = -R d s f k cannot be reached. Resuming, the set of couples (ω, R) for which the attitude stability is reached, (ω = 0, γ = 0), is given by:

E = {(ω, R) ∈ T SO(3) : ω = 0, R = PR d } (3.3) with P = {diag(1, 1, 1), diag(-1, 1, 1), diag(1, -1, 1), diag(1, 1, -1)} (3.4)
The block diagram of the closed loop attitude control is given in Figure 4 where the different blocks have been defined previously. In the following, the control torques will be addressed. The control forces are set such that they balance the gravity and ensure no movement along the three axes of the fixed frame F f .

Bounded attitude control

A bounded control torque computed using the direct measurements of on-board sensors is proposed. This control takes into consideration the saturation of the actuators and is simple to implement. First, the asymptotical stability of the rigid body is shown, then the proof is extended to show the stability of the FMAV. It is assumed that there exist at least two non-collinear vector observations embarked on the FMAV such that n ≥ 2.

Proposition 1. Consider the flapping micro-aerial vehicle's rotational dynamics described by (2.12-2.13). The attitude error γ between current and desired body's orientations is defined by (3.1) with n ≥ 2. The FMAV body's angular velocity ω = [ω 1 , ω 2 , ω 3 ] T is measured by three rate gyros. Define the bounded control τ = [τ 1 , τ2 , τ3 ] T by:

τj = -sat Nj λω j ρ j + λγ j , j ∈ {1, 2, 3} (3.5)
with ω, γ are respectively the averaged angular velocity and attitude error computed over a wingbeat period T , N j is the bound of the control torque component τj , sat Nj (•) is a classical saturation function, λ is a positive real parameter such that 0 < λ ≤ min(N j /2) and ρ j , j ∈ {1, 2, 3}, are positive scaling parameters. The control torque (3.5) stabilizes the FMAV at (ω, R) = (0, R d ) with a domain of attraction equal to

T SO(3)\{(0, R s R d )}, R s is given in (3.2).
Proof. Let's first prove the stability of a rigid body whose state vector is the averaged state vector of a flapping FMAV. The dynamics of the rigid body are then given by (2.12-2.13) where the angular velocity is denoted by ω, the rotation matrix by R, the control torque by τ and the attitude error by γ. Note that the derivative of an averaged vector observation sb k is linked to the vector by means of

ṡb k = sb k × ω = -ω × sb k (3.6)
Note also that the reference vectors are well known and fixed vectors; the average over a wingbeat period is equal to the vector. Consider first the following definite positive Lyapunov function:

V = 1 2 ωT J ω (3.7)
J is the inertia matrix of the body. V is trivially positive definite and radially unbounded. The derivative of (3.7) along the trajectories of the closed-loop system is given by Let Φ = {ω : |ω| ≤ ρ+ε} for some ε > 0 and ρ = [ρ 1 , ρ 2 , ρ 3 ] T . Outside Φ, i.e. |ω j | > ρ j +ε, using the unitary condition |γ j | ≤ 1, it follows that |λ ωj ρj + λγ j | ≥ λ ε ρj and that λ ωj ρj + λγ j and ωj have the same sign. Therefore,

V = ωT J ω = ωT (-ω × J ω) =0 +ω T τ = ωT τ =
V = - 3 j=1 ωj sat Nj λω j ρ j + λγ j ≤ - 3 j=1 λε ρ j |ω j | < - 3 j=1 λε(ρ j + ε) ρ j < -3λε < 0
Consequently, ω enters in Φ. During this time, γj cannot diverge since it is structurally unitary and therefore bounded.

Once in Φ, one has:

λω j ρ j + λγ j ≤ 2λ + λε ρ j
Taking ε sufficiently small and using the assumption that 2λ < N j , then:

λω j ρ j + λγ j ≤ N j
Consequently, sat Nj operates in a linear region, and the control torque becomes:

τj = λ ρ j (ω j + ρ γj ) (3.9)
In Φ, consider the Lyapunov function W defined by:

W = 1 2 ωT J ω + λ n n k=1 [1 -( Rs f k )(R d s f k )] (3.10)
W is a continuous and positive definite function on TSO(3) because W > 0 and W (0, Rd) = 0. Since sb k = Rs f k and s f k constant, the derivative of (3.10) is given by:

Ẇ = ωT J ω - λ n n k=1 ( ṡb k R d s f k ) (3.11)
By means of (3.1, 3.6) and the fact that ωT (-ω × J ω) = 0, (3.11) becomes:

Ẇ = ωT τ + λ n n k=1 ωT (s b k × R d s f k ) = ωT τ + λω T γ = 3 j=1 (ω j τj + λω j γj ) = 3 j=1 Ẇj (3.12)
Analyzing Ẇj for j ∈ {1, 2, 3}, one gets from τj in (3.9) and (3.12) that:

Ẇj = -ω j λω j ρ j + λγ j + λω j γj = - λ ρ j ω2 j and Ẇ = 3 j=1 Ẇj = -λ 3 j=1 ω2 j ρ j ≤ 0 (3.13)
The derivative of the Lyapunov function W is negative semidefinite. Recall that TSO(3) is compact. Therefore, for any initial condition (ω(0), R(0)) T SO(3), the set:

Ω = {(ω, R) ∈ T SO(3) : W (ω, R) ≤ W (ω(0), R(0))}
is a compact, positively invariant set of the closed-loop. From LaSalle Invariance Principle, it follows that all solutions that start in Ω converge to the largest invariant subset of Ω belonging to Ω. Ẇ (ω, R) ≡ 0 implies ω ≡ 0, then, substituting this last identity into the closed-loop system (2.12,2.13,3.5), one has

Ω = {(ω, R) ∈ T SO(3) : ω ≡ 0, γ ≡ 0}
with γ defined in (3.1) over the averaged measurements. Then, the largest invariant subset of Ω is given by E (3.3, 3.4).

Therefore, E defines the set of equilibrium points for which the Lyapunov function W (3.10) represents a minimum (W = 0) and a local maximum (W = 2λ) corresponding respectively to the rotation matrices R = R d and R = R s R d with R s given in (3.2). If the system is at one of these points at t 0 = 0, it will remain there for all t ≥ t 0 . The control law acts then to ensure the convergence of the closed-loop solutions, whose initial conditions do not verify (ω, R) = (0, R s R d ), to the stable equilibrium point given by (ω, R) = (0, R d ) and corresponding to W = Ẇ = 0. Therefore, the rigid body is asymptotically stable with a domain of attraction equal to

T SO(3) \ {(0, R s R d )}.
Once the asymptotic stability of the rigid body proven, the stability of the FMAV subject to the control torque (3.5) and computed over the averaged dynamics of the system (ω, R) will be presented. For that, the averaging theory recalled in §2.6 is used. It states that for a high frequency oscillating system, the averaged and time varying dynamics are very close and therefore, a stable equilibrium of the averaged dynamics is also a stable equilibrium of the time varying dynamics: ω -ω < k 1 T and R RT < k 2 T for k 1,2 > 0 and T the wingbeat period. In other words, the state vector (ω, R) of the FMAV, subject to the control torque (3.5), converges to the equilibrium point (0, R d ) with a domain of attraction equal to T SO(3) \ {(0, R s R d )}.

Simulations and robustness tests

Two strategies have been adopted in the literature to test the control laws developed for animal flapping flyers. The first is based on models of the flapping flight including aerodynamics, body and wing dynamics (Deng et al. 2006a;[START_REF] Rakotomamonjy | Longitudinal modelling and control of a flapping-wing micro aerial vehicle[END_REF][START_REF] Oppenheimer | Dynamics and control of a biomimetic vehicle using biased wingbeat forcing functions: Part I-aerodynamic model[END_REF][START_REF] Couceiro | Modeling and control of biologically inspired flying robots[END_REF]. Note that well developed models of flapping flight do not exist because the aerodynamics at low Reynolds numbers have not been yet perfectly identified [START_REF] Dudley | The biomechanics of insect flight: form, function, evolution[END_REF]. The second strategy consists on testing the control in real time on a prototype [START_REF] Finio | System identification and linear time-invariant modeling of an insect-sized flapping-wing micro air vehicle[END_REF][START_REF] Perez-Arancibia | First controlled vertical flight of a biologically inspired microrobot[END_REF] where simple models are proposed and identified. One should emphasize that at microscopic scale, only the vertical movement has been performed on a flapping flyer prototype [START_REF] Perez-Arancibia | First controlled vertical flight of a biologically inspired microrobot[END_REF] where the pitch movement is avoided by attaching the FMAV using some wires, the laser sensor measuring the vertical altitude is external as well as the power system. Vertical movement control of a 101 mg FMAV has been also presented in [START_REF] Duhamel | Altitude feedback control of a flapping-wing microrobot using an on-board biologically inspired optical flow sensor[END_REF]. Note that very challeging mass and size of FMAVs can be reached nowadays with the progress of microelectronics. The on-board electronic circuit developed within the OVMI/EVA project weighs 100 mg and is composed of a microprocessor, inertial and optic flow sensors. Piezoelectric actuators can reach 10 mg and are used within FMAVs [START_REF] Steltz | Dynamometer power output measurements of piezoelectric actuators[END_REF]. Other low size and weight microelectronic devices exist on the market for such applications. The power supply remains very difficult to embark.

In the present work, the first strategy is adopted. For simulations, the physical data of a dipteran insect are used [START_REF] Dudley | The biomechanics of insect flight: form, function, evolution[END_REF]. The mass is of 200 mg and the wingbeat frequency of 100 Hz. The maximum amplitudes for flapping and rotation angles are taken respectively to φ 0max = 60 • and ψ 0max = 90 • . The wingspan and wings surface are assumed respectively to 2L = 3 cm and 2S w = 1.14 cm 2 , allowing to generate a vertical ascendant movement using admissible flapping angles amplitudes.

Based on the defined numerical values, the admissible set for control torques Ω τ1,τ3 can be determined (2.17) and plotted (Figure 5). Ω τ1,τ3 is approximated to the largest ellipse E r that fits inside Ω τ1,τ3 (Figure 5) for computation simplification reasons. The control torques τ1 and τ3 respect then an ellipsoidal admissible set defined by

[τ 1 τ3 ]Q[τ 1 τ3 ] T ≤ 1 (4.1)
where Q is a diagonal definite positive matrix defining the ellipse's semi-axes denoted a r and b r . Practically, if τ1 ≥ a r (3.5), τ1 could be saturated to a r , and τ3 will be equal to zero. A privilege is given to the roll control in order to bring the FMAV to the horizontal plane. To avoid having a null yaw control torque and to give preference to the roll movement, 70% of a r is attributed to N 1 , τ3 will be calculated by (4.1) defining then a set Ω r (Figure 5).

The control torque (3.5) is applied to the FMAV to be validated in simulation. Three sensors are embarked on the FMAV to mimic a part of the sensory system of an insect (halteres, legs sensilla and magnetic sense): three rate gyros, a tri-axis accelerometer and a tri-axis magnetometer whose models are given by (2.3, 2.5, 2.4), respectively. The sensors measurements are not perfect, additive white gaussian noise is considered having a standard deviation of σ G = 10 -3 rad/s for the rate gyros, σ A = 10 -4 m/s 2 for the accelerometer and σ M = 0.01 mgauss for the magnetometer. Because of the low inertia of the FMAV, the developed control torque has also very small values of almost 10 -5 N • m (Figure 5). Since the control law is computed through a feedback of direct sensor measurements, the sensors should have a remarkable precision in order to guarantee a reasonable signal to noise ratio. Note that sensors with suitable precision or sensibility have already been developed [START_REF] Cheinet | Measurement of the sensitivity function in time-domain atomic interferometer[END_REF]. Even if the averaging of the angular velocity and attitude error reduce the noise influence, it can not be totally eliminated. The tuning parameters considered in simulations are given by (N 1 , N 2 , N 3 ) = (0.7 a r , 10 -5 , br ar

a 2 r -τ 2 1 ) with (a r , b r ) = (1.859 • 10 -5 , 5.843 • 10 -5 ), λ = 5 • 10 -6 , (ρ 1 , ρ 2 , ρ 3 ) = ( 2.2•10 -3 ar , 15.5, 5•10 -4 ar br √ a 2 r -τ 2 1
).

The evolution of the roll, pitch and yaw angles as well as the angular velocities measured by the rate gyros and the control torques are plotted in Figure 6. Note that the body's angles are not used in the control torque, they are plotted only to show the convergence of the FMAV's orientation. The initial orientation is (-40 • , -25 • , 50 • ) for the roll, pitch and yaw angles, respectively. The initial angular velocity is null. The sensor measurements are presented in Figure 7 and the wing angle amplitude envelops in Figure 8. The stability is reached in a sufficiently fast time which makes the control law suitable for real-time implementation on a FMAV. Moreover, it presents comparable values to those observed in true insects [START_REF] Dudley | The biomechanics of insect flight: form, function, evolution[END_REF]. The rate gyros noise is detectable in Figure 6. The accelerometer's noise is more detectable than that of the magnetometer. Their measurements converge to s f A = [0, 0, -1] T and s

f M = [ 1 2 , 0, - √ 3
2 ] at the equilibrium. Recall that the FMAV is stabilized in hovering mode. Therefore, the flapping angles converge to the value that generates an aerodynamic lift balancing its weight. The rotation angles converge to zero so that no translational movement is generated.

The robustness of the control law is tested with respect to external disturbances. These disturbances simulate wind or rain drops affecting the FMAV and creating a body torque of (1.2 10 -5 , 2 10 -5 , 1.2 10 -5 ) N • m, applied at t = 1.5 s during 10 wingbeat periods. Noting that a rain drop weighs about 5 • 10 -6 N and the FMAV body's inertia is of 10 -8 Kg • m 2 , the disturbance applied has a great effect on the body's movement. It will destabilize the FMAV all over the disturbance taking it far away from the equilibrium. The control law acts later on to bring the FMAV back to stability. The evolution of the roll, pitch and yaw angles, angular velocities and control torques is plotted in Figure 9 zoomed to the disturbance zone in Figure 10. Note that the bound of the yaw torque depends on the value of the roll torque as explained previously. The reference sensor measurements are given in Figure 11 and the wing angle amplitudes in Figure 12. Notice that the flapping angle of the left wing is bounded during the disturbance at φ max = 60 • avoiding the saturation of the actuator and emphasizing the development of a bounded control law. One should also note that an insect subject to such a high disturbance will lose all control of its position and orientation and regain it only when the disturbance is over. The measurement along the roll axis is plotted with dot-dashed blue line, the pitch axis with dashed red line and yaw axis with continuous black line.

Conclusions and future works

Even if the gap between the flapping-wing animal flight and its technological reproduction is still very large, the robotic and control communities are exerting big effort to develop flapping micro aerial vehicles that mimic the nature's flight the most accurately possible. The main contribution of this work is the development of a control law stabilizing the FMAV's attitude. It is based directly on the measurements of some embedded sensors without the need of computing the orientation. This strategy reproduces the insects' one to stabilize their orientation based on their halteres, legs sensilla and magnetic sense. The control torque is bounded, allowing to take into account the amplitude bounds of the flapping wings in order to avoid the saturation of the actuators and ensure the stability of the body. Note that it allows different saturation bounds along the three axes. It is also not restricted to symmetric bodies and independent of the inertia matrix. Moreover, it is simple to compute and is adaptable for real time implementation. As shown in simulations, the control law is robust with respect to external disturbances. The boundedness of the control torque helps preserving the piezoelectric actuators not saturated and guarantee then their linear behavior even in hard conditions. Future works will consider the development of bounded control force aiming to control the FMAV's trajectory based on sensor's measurements. In fact, insects can determine the sun direction using their ocelli. Moreover, the light polarization direction can be determined using the compound eyes. Based on polarized light compasses, for example, one can determine the direction of flight and couple it with the attitude control in order to ensure a movement in the three-dimensional space.

Figure 1 .

 1 Figure 1. Degrees of freedom of a wing: flapping, lagging (elevation) and feathering

Figure 2 .

 2 Figure 2. Coordinate frames: inertial fixed frame F f ( x f , y f , z f ), body attached frame F b ( x b , y b , z b ) and left wing attached frame F w l ( r l , t l , n l )

Figure 3 .

 3 Figure 3. Shape of the wing angle evolution through an illustrative example: The flapping (top) and rotation (bottom) wing angles, the reference angles in dashed red line and the real angles (delivered by the actuators) in continuous blue line.

  11) with p b l,r is the position of the left or right wing aerodynamic center in the body frame computed by p b l,r = R b l,r p w . It is considered, in the present work, that the aerodynamic and pressure centers coincide. The left and right wing velocities, expressed in the body frame, are computed by deriving the wing center of pressure positions v b l,r = ṗb l,r . Their projections in the wing frames are obtained by a simple rotation v w l,r = R l,r b v b l,r .

Figure 4 .

 4 Figure 4. The block diagram of the FMAV's attitude control based on sensors measurements.

Figure 5 .

 5 Figure 5. Yaw torque versus roll torque, defining the saturation set Ωτ 1 ,τ 3 approximated to an ellipse Er then to a set Ωr.

Figure 6 .Figure 7 .

 67 Figure 6. The attitude (left), angular velocity (middle) of the FMAV going from initial roll, pitch and yaw angles (-40 • , -25 • , 50 • ) and null angular velocity and the control torques (right) applied to the FMAV.

Figure 8 .

 8 Figure 8. The envelops of the left and right wing angles.

Figure 9 .Figure 10 .Figure 11 .Figure 12 .

 9101112 Figure 9. Robustness with respect to external disturbances: The (left), angular velocity (middle) and the control torques (right) applied to the FMAV.

† The definition is given by the DARPA (Defense Advanced Research Projects Agency).

† The paper lies within the scope of the project OVMI/EVA sponsored by the French National Research Agency.