B Bergeot

S Bellizzi

B Cochelin

Passive suppression of helicopter ground resonance using nonlinear energy sinks attached on the helicopter blades

Keywords: Helicopter ground resonance, Passive control, Nonlinear energy sink, Complex multi-blade coordinate transformation, Relaxation oscillations, Strongly modulated response

This paper investigates the passive control of a rotor instability named helicopter Ground Resonance (GR). The passive device consists of a set of essential cubic nonlinear absorbers named Nonlinear Energy Sinks (NES) each of them positioned on a blade. A dynamic model reproducing helicopter GR instability is presented and transformed to a time-invariant nonlinear system using a multi-blade coordinate transformation based on Fourier transform mapping the dynamic state variables into a non-rotating reference frame. Combining complexification, slow/fast partition of the dynamics and averaging procedure, a reduced model is obtained allowed us to use the so-called geometric singular perturbation analysis to characterize the steady state response regimes. As in the case of a NES attached to the fuselage, it is shown that under suitable conditions, GR instability can be completely suppressed, partially suppressed through periodic response or strongly modulated response. Relevant analytical results are compared, for validation purposes, to direct integration of the reference and reduced models.

Introcuction

A Nonlinear Energy Sink (NES) is a passive nonlinear absorber which introduces non linear coupling between the main structure and the absorber for controlling the vibratory response of the structure under external excitation. A NES consists on a nonlinear oscillator with an essentially nonlinear restoring force and generally a small mass compared to the main structure. The objective of the NES is to capture the energy of the main structure and to dissipate it. This is done through the Targeted Energy Transfer (TET) concept. TET is based on the capacity of the essential nonlinear oscillator to be resonant at any frequencies giving possible tuning with the main structure to be controlled. TET has been widely studied in the literature [START_REF] Vakakis | Energy pumping in nonlinear mechanical oscillators: Part II -Resonance capture[END_REF][START_REF] Gourdon | Nonlinear energy pumping under transient forcing with strongly nonlinear coupling: Theoretical and experimental results[END_REF][START_REF] Vakatis | Nonlinear Targeted Energy Transfer in Mechanical and Structural Systems[END_REF][START_REF] Starosvetsky | Strongly modulated response in forced 2dof oscillatory system with essential mass and potential asymmetry[END_REF].

NES also represents an alternative device for control systems with self-excitation. The possible suppression of the limit cycle oscillations of a van der Pol oscillator coupled to a NES is demonstrated numerically in [START_REF] Lee | Suppression of limit cycle oscillations in the van der Pol oscillator by means of passive non-linear energy sinks[END_REF] and the self-excitation response regimes are investigated theoretically in [START_REF] Gendelman | Bifurcations of self-excitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF] where an asymptotic analysis of the system related to slow/super-slow decomposition of the averaged flow reveals periodic responses, global bifurcations of different types and basins of attraction of various self-excitation regimes. A series of papers [START_REF] Lee | Suppression aeroelastic instability using broadband passive targeted energy transfers, part 1: Theory[END_REF][START_REF] Lee | Suppression aeroelastic instability using broadband passive targeted energy transfers, part 2: Experiments[END_REF][START_REF] Lee | Enhancing the robustness of aeroelastic instability suppression using multi-degree-of-freedom nonlinear energy sinks[END_REF][START_REF] Gendelman | Asymptotic analysis of passive nonlinear suppression of aeroelastic instabilities of a rigid wing in subsonic flow[END_REF] are dedicated to aeroelastic instability. It is demonstrated that a NES coupled to a rigid wing in subsonic flow can partially or even completely suppress aeroelastic instability. The suppression mechanisms are investigated numerically in [START_REF] Lee | Suppression aeroelastic instability using broadband passive targeted energy transfers, part 1: Theory[END_REF] and [START_REF] Lee | Enhancing the robustness of aeroelastic instability suppression using multi-degree-of-freedom nonlinear energy sinks[END_REF] for one and multi-DOF NES. Several aspects of the suppression mechanisms are validated expermentally in [START_REF] Lee | Suppression aeroelastic instability using broadband passive targeted energy transfers, part 2: Experiments[END_REF]. Finally an asymptotic analysis is reported in [START_REF] Gendelman | Asymptotic analysis of passive nonlinear suppression of aeroelastic instabilities of a rigid wing in subsonic flow[END_REF] demonstrating the existence of the three passive suppression mechanisms based on TET. Suppression of aeroelastic instability of a long span bridge model is also considered in [START_REF] Vaurigaud | Passive control of aeroelastic instability in a long span bridge model prone to coupled flutter using targeted energy transfer[END_REF] whereas [START_REF] Luongo | Aeroelastic instability analysis of nes-controlled systems via a mixed multiple scale/harmonic balance method[END_REF] is dedicated to a general class of nonlinear multi degree of freedom system including aeroelastic forces. In [START_REF] Nankali | Suppression of regenerative instabilities by means of targeted energy transfers[END_REF][START_REF] Gourc | Quenching chatter instability in turning process with a vibro-impact nonlinear energy sink[END_REF] it is shown that a NES is able to quench chatter instability in turning process whereas in [START_REF] Viguié | Using passive nonlinear targeted energy transfer to stabilize drill-string systems[END_REF] the stabilization of drill-string systems is considered. Note that in [START_REF] Habib | Suppression of limit cycle oscillations using the nonlinear tuned vibration absorber[END_REF], it is demonstrated that a nonlinear tuned vibration absorber possessing a linear spring and a nonlinear spring whose mathematical form is determined according the nonlinearity in the host system can suppress limit cycle oscillations.

In this paper, we focus on the instability in helicopter rotor named Ground Resonance (GR) problem. GR is related to the coupling of the rotor blade in-plane motion with the airframe motion on its landing gear. The standard reference of the GR analysis is the paper by Coleman and Feingold [START_REF] Coleman | Theory of self excited mechanical oscillations of helicopter rotor with hinged blades[END_REF] where it is established, considering an isotropic rotor, that GR is due to a frequency coalescence between a lag mode and the fuselage mode. The range of rotor speeds for which this frequency coalescence occurs is predicted analytically. More references can be found in [START_REF] Johnson | Helicopter theory[END_REF][START_REF] Bramwell | Bramwell's helicopter dynamics[END_REF][START_REF] Kunz | Nonlinear analysis of helicopter ground resonance[END_REF][START_REF] Krysinski | Instabilité mécanique: contrôle actif et passif[END_REF] and a recent analysis of helicopter GR with asymmetric blades can be found in [START_REF] Sanches | Parametrically excited helicopter ground resonance dynamics with high blade asymmetries[END_REF]. Traditionally, GR instability is prevented by increasing damping in the landing gear and blades [START_REF] Done | A simplified approach to helicopter ground resonance[END_REF][START_REF] Tang | Influence of nonlinear blade damping on helicopter ground resonance instability[END_REF]. A robustness analysis is discussed in [START_REF] Quaranta | Lead-lag damper robustness analysis for helicopter ground resonance[END_REF]. Active control of GR has been also studied in [START_REF] Krysinski | Instabilité mécanique: contrôle actif et passif[END_REF].

The initial application of a NES to control GR instability was considered in [START_REF] Bergeot | Analysis of steady-state response regimes of a helicopter ground resonance model including a non-linear energy sink attachment[END_REF][START_REF] Bergeot | Passive suppression of helicopter ground resonance instability by means of a strongly nonlinear absorber[END_REF]. A theoretical/numerical analysis of the steady-state responses of a helicopter model with a minimum number of degrees of freedom that can reproduce Helicopter Ground Resonance instability when a NES is attached on the fuselage in an ungrounded configuration was performed. The system was first simplified using successively Coleman transformation [START_REF] Coleman | Theory of self excited mechanical oscillations of helicopter rotor with hinged blades[END_REF] and binormal transformation [START_REF] Done | A simplified approach to helicopter ground resonance[END_REF] followed by a complexification-averaging method together with geometric singular perturbation. Four steady-state responses are highlighted and explained analytically: complete suppression, partial suppression through strongly modulated response, partial suppression through periodic response and no suppression of the Helicopter Ground Resonance.

In this study, the same problem is considered but now the NES on the fuselage is removed and it is replaced by a set of NES, each of them been positioned on a blade. A similar configuration has been considered in [START_REF] Bab | Lateral vibration attenuation of a rotor under mass eccentricity force using non-linear energy sink[END_REF][START_REF] Bab | Vibration mitigation of a rotating beam under external periodic force using a nonlinear energy sink (nes)[END_REF] where the efficiency of NES on rotating part of a rotor has been analyzed under mass eccentricity force [START_REF] Bab | Lateral vibration attenuation of a rotor under mass eccentricity force using non-linear energy sink[END_REF] and under an external forcing [START_REF] Bab | Vibration mitigation of a rotating beam under external periodic force using a nonlinear energy sink (nes)[END_REF]. One of the advantages of this configuration is to reduce the weight due to the additional absorbers. Here, response regimes are investigated extending the theoretical/numerical analysis developed in [START_REF] Bergeot | Analysis of steady-state response regimes of a helicopter ground resonance model including a non-linear energy sink attachment[END_REF] using a Fourier coordinate transformation as in [START_REF] Ghazavi | Bifurcation and nonlinear analysis of nonconservative interaction between rotor and blade row[END_REF]. In contrast to the Colman transformation, this multi-blade coordinate transformation can be easily combined with the complexification-averaging method together with geometric singular perturbation.

The next section introduces the helicopter model including a fuselage and four blades with a NES attached on each blade. Only one direction of the fuselage motion, the lag motion of each blade and the lag motion of each associated NES are considered. Based on Fourier coordinate transform, the model is reduced to a time-invariant nonlinear system involving complex variables. In Sect. 3, linear stability analysis is performed on this reduced model showing that this model is able to reproduce GR phenomenon. In Sect. 4, the analytical procedure [START_REF] Bergeot | Analysis of steady-state response regimes of a helicopter ground resonance model including a non-linear energy sink attachment[END_REF] based on averaging method together with geometric singular perturbation theory is adapted to the reduced system. It is shown that it is able to characterize the situations for which trivial solution is unstable. In Sect. 5, numerical analysis is performed and discussed considering analytical results validated from direct integration of the reference model and the reduced model. Finally in Sect. 6, few words about the comparison between blade and fuselage NES attachments are expounded.

The model

Initial equations of motion

The helicopter model studied is shown Fig. 1. It consists of a fuselage and a 4-blades rotor rotating at a constant speed Ω. This model is very similar to that described for example in [START_REF] Johnson | Helicopter theory[END_REF][START_REF] Bramwell | Bramwell's helicopter dynamics[END_REF][START_REF] Krysinski | Instabilité mécanique: contrôle actif et passif[END_REF].

The fuselage is a simple mass-spring damped system with mass m y , spring constant k y , viscous damper of damping coefficient c y and the translational DOF y. It is assumed that the center of inertia G f of the fuselage at rest coincides with the origin O of the earth-fixed system of coordinates (O, x 0 , y 0 , z 0 ). Each blade is assumed to be a point mass G δ,i (with i ∈ [START_REF] Vakakis | Energy pumping in nonlinear mechanical oscillators: Part II -Resonance capture[END_REF][START_REF] Starosvetsky | Strongly modulated response in forced 2dof oscillatory system with essential mass and potential asymmetry[END_REF]) with mass m δ linked to the axis (O, z 0 ) with a bar without mass of length L δ and an articulation with torsional spring k δ and viscous damper of damping coefficient c δ . Only lag motions of the blades, characterized by the lag angle δ i (with i ∈ [START_REF] Vakakis | Energy pumping in nonlinear mechanical oscillators: Part II -Resonance capture[END_REF][START_REF] Starosvetsky | Strongly modulated response in forced 2dof oscillatory system with essential mass and potential asymmetry[END_REF]), are taken into account.

This simple model is used to study the effect of attaching a NES in ungrounded configuration on each blade of the helicopter (as shown Fig. 2) in order to control lagging motion and consequently ground resonance instability.

Each NES is a nonlinear-spring mass damped system (see Fig. 2) with mass m h at the mass point G ϕ,i (with i ∈ [START_REF] Vakakis | Energy pumping in nonlinear mechanical oscillators: Part II -Resonance capture[END_REF][START_REF] Starosvetsky | Strongly modulated response in forced 2dof oscillatory system with essential mass and potential asymmetry[END_REF]), linear spring constant k 1h , cubic spring constant k 3h , and viscous damper of damping coefficient c h . Each NES is placed at distance L ϕ from the axis (O, z 0 ) and it is characterized with only the lag angle DOF ϕ i which is related to the relative displacement h i by

h i = L ϕ sin(ϕ i -δ i ).
Starting from the position of the center of inertia and the corresponding velocity of the i-th blade (respectively the i-th NES) in the plane (O, x 0 , y 0 ), defined as

OG δ,i = L δ cos (Ωt + ψ i + δ i ) y + L δ sin (Ωt + ψ i + δ i ) , (1)
and,

v δ,i =   -L δ Ω + δi sin (Ωt + ψ i + δ i ) ẏ + L δ Ω + δi cos (Ωt + ψ i + δ i )   , (2)
and,

OG ϕ,i = L ϕ cos (Ωt + ψ i + ϕ i ) y + L ϕ sin (Ωt + ψ i + ϕ i ) , ( 3
)
v ϕ,i = -L ϕ (Ω + φi ) sin (Ωt + ψ i + ϕ i ) ẏ + L ϕ (Ω + φi ) cos (Ωt + ψ i + ϕ i ) (4)
with,

ψ i = - π 2 (i -1) (5)
the equations of motion of the whole system is derived using Lagrange equations from the kinetic energy T , the potential energy V and the Rayleigh dissipation function D written as

T = 1 2 m y ẏ2 + 4 i=1 1 2 m δ v 2 δ,i + 4 i=1 1 2 m ϕ v 2 ϕ,i , (6)
V = 1 2 k y y 2 + 4 i=1 1 2 k δ δ 2 i + 4 i=1 1 2 k 1 sin 2 (ϕ i -δ i ) + 4 i=1 1 4 k 3 sin 4 (ϕ i -δ i ), (7)
D = 1 2 c y ẏ2 + 4 i=1 1 2 c δ δ2 + 4 i=1 1 2 c ϕ sin 2 ( φi -δi ), (8)
with

k 1 = L 2 ϕ k 1h , k 3 = L 4 ϕ k 3h and c ϕ = L 2 ϕ c h . Assuming small angles (|ϕ i | 1 and |δ i | 1)
and also small angle differences (|ϕ i -δ i | 1, the equations of motion reduce to a nonlinear differential system with periodic coefficient

                                     (m y + 4(m δ + m ϕ )) ÿ + c y ẏ + k y y + M δ 4 i=1 δi -Ω 2 δ i cos (Ωt + ψ i ) -2Ω δi sin (Ωt + ψ i ) + M ϕ 4 i=1 φi -Ω 2 ϕ i cos (Ωt + ψ i ) -2Ω φi sin (Ωt + ψ i ) = 0 (9a) I δ δi + c δ δi + c ϕ ( δi -φi ) + k δ δ i + k 1 (δ i -ϕ i ) + k 3 (δ i -ϕ i ) 3 + M δ ÿ cos (Ωt + ψ i ) = 0, i = 1, 4 (9b)
I ϕ φi -c ϕ ( δi -φi ) -k 1 (δ i -ϕ i ) -k 3 (δ i -ϕ i ) 3 + M δ ÿ cos (Ωt + ψ i ) = 0, i = 1, 4. (9c)
where " ˙" denotes the derivative with respect to time t and M δ = m δ L δ and

I δ = m δ L 2 δ (respectively M ϕ = m ϕ L ϕ and I ϕ = m ϕ L 2
ϕ ) are the static moment and the moment of inertia of each blade (respectively NES).

It is important to note that the cosine and sine functions do not now depend on the angle variables δ i and ϕ i (see for the record Eqs. ( 1)-( 4)) as the result of the small angle assumption (all the cosine and sine terms have been expanded in a first-order Taylor series around the trivial equilibrium position).

Scaled equations of motion: the Reference Model (RefM)

Our objective is to show that the lightweight ungrounded NES can be effective passive absorbers and local energy dissipators to mitigate instabilities.

Defining the lightweight of a NES with respect to a blade as the ratio of the moments of inertia

= I ϕ I δ . ( 10
)
and introducing the barycentric coordinates v i (t) and w i (t) as

v i = δ i + ϕ i and w i = δ i -ϕ i , for i = 1, • • • , 4, (11)
Eqs. ( 9) become

                                                       ÿ + λ y ẏ + ω 2 y y + 1 2(1 + ) 4 i=1 (12a) Ω sin (Ωt + ψ i ) -(S d + T d )(Ω + Ω + 2 vi ) + 2 (T d -S d ) ẇi + cos (Ωt + ψ i ) (S d + T d )(v i -Ω 2 v i ) -(T d -S d )( ẅi -Ω 2 w i ) = 0 (12b) vi + λ δ vi + ẇi + 1 + ω 2 δ v i + w i + 1 + (S c + T c )ÿ cos (Ωt + ψ i ) = 0, i = 1, 4 (12c) ẅi + λ δ vi + ẇi + 1 + ω 2 δ v i + w i + 1 -(T c -S c )ÿ cos (Ωt + ψ i ) +µ(1 + ) ẇ + α 1 (1 + )w i + α 3 (1 + )w 3 i = 0, i = 1, 4 (12d)
where the following notations defining rescaled parameters

ω 2 y = k y / (m y + 4(m δ + m ϕ )) , ( 13a
)
ω 2 δ = k δ /I δ , ( 13b
)
λ y = c y / (m y + 4(m δ + m ϕ )) , ( 13c
)
λ δ = c δ /I δ , ( 13d
) µ = c ϕ /I δ , ( 13e
)
α 1 = k 1 /I δ , ( 13f
)
α 3 = k 3 /I δ , ( 13g
)
S c = M δ /I δ = 1/L δ , ( 13h
) 2 T c = M ϕ /I δ , ( 13i
)
S d = 2M δ / (m y + 4(m δ + m ϕ )) , (13j) 2 T d = 2M ϕ / (m y + 4(m δ + m ϕ )) , (13k)
have been used.

In the sequel, we will be assumed that

1 with λ y , λ δ , µ, α 1 , α 3 , S d , T d ,S c , T c ∼ O(1)
. We point that expression (13i) defining M ϕ /I δ = /L ϕ is already of the order of . In the same way, if we reason with scaled parameters Sd = 2M δ / (m y + 4(m δ + m ϕ )) and Td

= 2M ϕ / (m y + 4(m δ + m ϕ )), therefore, Td = Sd (M ϕ /M δ ) = Sd (L δ /L ϕ ) = S d 2 (L δ /L ϕ ).
Therefore, 2 appears in Eqs. (13i) and (13k) resulting of the assumption that the lengths L δ and L ϕ are the same order of magnitude. Eqs. [START_REF] Luongo | Aeroelastic instability analysis of nes-controlled systems via a mixed multiple scale/harmonic balance method[END_REF] defines the Reference Model (RefM).

The time-invariant system

To develop an analytic method for characterizing the motions of the helicopter when Ground Resonance (GR) can occur, RefM has to be reduced (and/or simplified) to a time-invariant system. Classically, a Coleman transform [START_REF] Coleman | Theory of self excited mechanical oscillations of helicopter rotor with hinged blades[END_REF] is used. Coleman transform consists in a change of variables which transforms individual motions of the blades (described by the lagging angles) into collective motions described by the so-called Coleman coordinates. Coleman transform applied to a linear system with periodic coefficients gives a linear system with constant coefficients.

Here a complex transformation, first proposed in [START_REF] Ghazavi | Bifurcation and nonlinear analysis of nonconservative interaction between rotor and blade row[END_REF] and used for example in [START_REF] Bab | Lateral vibration attenuation of a rotor under mass eccentricity force using non-linear energy sink[END_REF], is preferred to Coleman transform because this complex transformation deals easily with nonlinearities. Indeed, we can first work in the rotating reference frame applying the following complex transformation of the variables v i and w i (for i = 1 to 4)

v i (t) = 4-1 k=0 a k (t)e -jkψ i , with i = 1, • • • , 4 (14)
given equivalently

a k (t) = 1 4 4 i=1 v i (t)e jkψ i , with k = 0, • • • , 3 (15)
and,

w i (t) = 4-1 k=0 b k (t)e -jkψ i , with i = 1, • • • , 4 (16)
given equivalently

b k (t) = 1 4 4 i=1 w i (t)e jkψ i , with k = 0, • • • , 3 (17)
with j 2 = -1. The new variable a k and b k (for k = 0 to 3), called modal coordinates, are now complex and satisfy the following relations a 4-k = a * k and b 4-k = b * k where " * " is the usual notation for the complex conjugate.

Substituting Eqs. ( 14) and ( 16) into Eqs. [START_REF] Luongo | Aeroelastic instability analysis of nes-controlled systems via a mixed multiple scale/harmonic balance method[END_REF], it can be observed that the modal coordinates a 0 and a 2 are uncoupled from a 1 , a

λ y = 0.2, ω y = 1, S d = 1, T d = 0.1, S c = 3, T c = 1, λ = 0.2, ω δ = 2, µ = 0.2, α 1 = 0.5 and = 0.1.
Finally, equations of motion are written in the inertia reference frame using the variables A and B defined as follow

A(t) = a(t)e jΩt and B(t) = b(t)e jΩt . ( 18
)
where A and B are also complex variables. The equations of motion associated to y, A and B reduce to

                                 ÿ + λ y ẏ + ω 2 y y + S d + T d 1 + Ä + Ä * + 2 S d -T d 1 + B + B * = 0 (19a) Ä + λ δ 1 + -2jΩ Ȧ + -Ω(Ω + (iλ δ + Ω)) + ω 2 δ 1 + A + 1 2 (S c + T c ) ÿ + 2 λ δ 1 + Ḃ + -j λ δ Ω + ω 2 δ 1 + B = 0 (19b) B + µ -2jΩ + µ + λ δ + 1 Ḃ + α 1 (1 + ) 2 + Ω -j 2 λ δ + (1 + ) 2 µ -(1 + )Ω + ω 2 δ 1 + B + 1 2 (S c -T c ) ÿ + λ δ 1 + Ȧ + -j λ δ Ω + ω 2 δ 1 + A + 3α 3 (1 + ) |B| 2 B = 0. ( 19c
)
High harmonic terms in e 4jtΩ appearing in Eqs. (19c) have been ignored. The resulting equations are now time-invariant system of equations. They will be used to perform linear stability analysis of the trivial solution and also as starting point of the nonlinear analysis.

Preliminary results

Ground resonance analysis

Helicopter ground resonance is a consequence of a dynamic instability which appears in the RefM for given values of the rotor speed Ω. Considering the time-variant (or parametrically excited) systems ( 9) or [START_REF] Luongo | Aeroelastic instability analysis of nes-controlled systems via a mixed multiple scale/harmonic balance method[END_REF] the instability is a parametric combination resonance of the summed type, see e.g. [START_REF] Nayfeh | Nonlinear Oscillations[END_REF]. This king of instability can be analyzed using classical methods for linear timevariant ordinary differential equations as Floquet theory, see e.g. [START_REF] Nayfeh | Nonlinear Oscillations[END_REF][START_REF] Nayfeh | Applied nonlinear dynamics[END_REF]. However, is it easier (and equivalent) to perform the stability analysis on the time-invariant system [START_REF] Bramwell | Bramwell's helicopter dynamics[END_REF] in which the instability becomes a mode coupling instability. Usually, the ground resonance phenomenon is analyzed on a time-invariant system resulting of Coleman transformation (see [START_REF] Johnson | Helicopter theory[END_REF] for more details) instead of the used complex transformation defined by Eqs. [START_REF] Gourc | Quenching chatter instability in turning process with a vibro-impact nonlinear energy sink[END_REF][START_REF] Viguié | Using passive nonlinear targeted energy transfer to stabilize drill-string systems[END_REF][START_REF] Habib | Suppression of limit cycle oscillations using the nonlinear tuned vibration absorber[END_REF][START_REF] Coleman | Theory of self excited mechanical oscillations of helicopter rotor with hinged blades[END_REF].

It is easy to show that the only fixed point of Eqs. ( 19) is the trivial solution y = A = A * = B = B * = 0. As usual, stability of the fixed point is found by looking the sign of the eigenvalues real parts of the Jacobian matrix of the vector function F 1 evaluated at the trivial equilibrium. The vector function F 1 characterizes the system of Eqs. [START_REF] Bramwell | Bramwell's helicopter dynamics[END_REF] when it is formally written in state-space form

U = F 1 (U), with U = y A A * B B * ẏ Ȧ Ȧ * Ḃ Ḃ * t . ( 20
)
The evolution, with respect to the rotor speed Ω, of the 10 complex conjugates eigenvalues of the Jacobian matrix of the vector function F 1 are plotted in Fig. 3. Near Ω = ω y + ω δ , a phenomenon of frequency coalescence is observed, the real part of one of the eigenvalues becomes positive and a dynamic instability occurs; this is the helicopter ground resonance. Due to the presence of the NES attachments, a second frequency coalescence appears near Ω = ω y + √ α 1 ≈ 1.7. This can be a problem if one wants to build such a device. However, the assumption that the lengths L δ and L ϕ are the same order of magnitude leading to Eq. (13k) implies a weak coupling between the NES and the fuselage. As a result, in this situation, this frequency coalescence does not produce a dynamic instability. Indeed, no real part of the eigenvalues become positive (see Fig. 3(b)).

Some steady-state response regimes

Using numerical integration of the RefM, Eqs. [START_REF] Luongo | Aeroelastic instability analysis of nes-controlled systems via a mixed multiple scale/harmonic balance method[END_REF], four different types of response regimes which may be highlighted. They are classified into two categories depending on the fact that the trivial solution of the RefM, calculated in Sect. 3.1, is stable or not:

• The trivial solution of the RefM is stable:

-Complete suppression. In this case, the additional damping due to the NES attachment stabilizes the system and the GR instability is completely suppressed.

• The trivial solution of RefM is unstable:

-Partial suppression through Periodic Response (PR). In this case, the steady-state response regime is periodic with frequency close to ω y1 .

-Partial suppression through Strongly Modulated Response (SMR). In this case, the steadystate response regime is a quasiperiodic regime which exhibits a "fast" component with frequency close to ω y and a "slow" component corresponding to the envelope of the signal. The term "Strongly modulated response" has been introduced by Starosvetsky and Gendelman [START_REF] Starosvetsky | Strongly modulated response in forced 2dof oscillatory system with essential mass and potential asymmetry[END_REF] for the study of a harmonic forced linear system coupled to a NES.

-No suppression of GR. The NES is not able to maintain stable steady-state regimes. We observe exponential growth of the system.

These four responses are also observed by Bergeot et al. [START_REF] Bergeot | Analysis of steady-state response regimes of a helicopter ground resonance model including a non-linear energy sink attachment[END_REF] when one NES is attached on the fuselage of the helicopter.

Theoretical study: asymptotic analysis

Our objective is to develop a similar procedure as in Bergeot et al. [START_REF] Bergeot | Analysis of steady-state response regimes of a helicopter ground resonance model including a non-linear energy sink attachment[END_REF] to characterize the periodic solution with period ω y for Ω in the neighborhood of ω y +ω δ . The proposed method combines reduction by complexfication-averaging (CA-X) procedure [START_REF] Manevitch | Complex representation of dynamics of coupled nonlinear oscillators[END_REF][START_REF] Vakatis | Nonlinear Targeted Energy Transfer in Mechanical and Structural Systems[END_REF] and analysis using geometric singular perturbation theory (GSPT) [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF][START_REF] Jones | Geometric singular perturbation theory[END_REF][START_REF] Desroches | Mixed-mode oscillations with multiple time scales[END_REF].

The slow-flow

Starting from Eqs. [START_REF] Bramwell | Bramwell's helicopter dynamics[END_REF], complexification2 can be achieved re-writing the variable y as

ξ = ẏ + jω y y, ( 21
)
giving reciprocally,

y = ξ -ξ * 2jω y , ẏ = ξ + ξ * 2 and ÿ = ξ - jω y 2 (ξ + ξ * ) . ( 22
)
Then, averaging method is developed in the way of the complexification-averaging method first introduced by Manevitch [START_REF] Manevitch | Complex representation of dynamics of coupled nonlinear oscillators[END_REF] and discussed in detail by Vakakis et al. [START_REF] Vakatis | Nonlinear Targeted Energy Transfer in Mechanical and Structural Systems[END_REF]. Only the main steps are given here after.

We assume that the variable ξ, A and B may be broken down into fast and slow components introducing the following representation

ξ 1 = φ 1 (τ )e jωyt , A = φ 2 (τ )e jωyt , B = φ 3 (τ )e jωyt , ( 23
)
where τ = t (τ defines the super-slow time scale, see hereafter for more details about this designation) and φ i (with i ∈ [1, 3]) is the complex slow modulated amplitude of the fast component e jωyt .

The following steps are performed to obtain the slow flow of Eqs. ( 19):

1. Eqs. ( 22) and ( 23) are substituted into Eqs. (19);

2. An averaging over one period of the frequency ω y is performed;

3. The resulting equations are written with respect to the super-slow time scale, taking into account the fact that

dφ i (τ ) dt = dφ i (τ ) dτ and d 2 φ i (τ ) dt 2 = 2 dφ i (τ ) dτ 2 , ( 24
)
4. The equations are rewritten with respect to the slow time scale t;

5. A first-order Taylor series around = 0 is performed.

Finally, Eqs. ( 19) are reduced as

       φ1 = f 1 (φ 1 , φ 2 , φ 3 ) (25a) φ2 = f 2 (φ 1 , φ 2 , φ 3 ) (25b) φ3 = f 3 (φ 1 , φ 2 , φ 3 , ) (25c)
where the functions f 1 , f 2 and f 3 are given by

                                 f 1 (φ 1 , φ 2 , φ 3 ) = - λ y 2 φ 1 + ω 2 y S d φ 2 (26a) f 2 (φ 1 , φ 2 , φ 3 ) = - 1 2j - jω y S c 4ω δ φ 1 + (ω δ + 2σ + jλ δ )φ 2 -ω δ φ 3 (26b) f 3 (φ 1 , φ 2 , φ 3 , ) = - ω 2 δ µ -2jω δ φ 2 + ω 2 δ -α 1 + jµω δ µ -2jω δ φ 3 - 3α 3 µ -2jω δ φ 3 |φ 3 | 2 - µ -2jω δ jω y (T c + S c ) 4 φ 1 + ω 2 δ -jλ δ ω δ φ 2 + ω 2 δ + α 1 -2σω δ -jµ(ω δ + σ) φ 3 + 3α 3 φ 3 |φ 3 | 2 . ( 26c
)
Note that Ω has been replaced by

Ω = ω y + ω δ + σ, (27)
introducing the detuning term σ (with σ ∼ O( 1)) to analyze the system for Ω around ω y + ω δ . Eqs. ( 25) describe the slow-flow of the RefM [START_REF] Luongo | Aeroelastic instability analysis of nes-controlled systems via a mixed multiple scale/harmonic balance method[END_REF] and fixed points of this slow-flow trivially characterize periodic responses of the RefM model. It will be called the Reduced Model (RedM).

Geometric singular perturbation analysis

Eqs. ( 25) can be transformed by switching from the slow time scale t to the super-slow time scale τ = t as

     φ 1 = f 1 (φ 1 , φ 2 , φ 3 ) (28a) φ 2 = f 2 (φ 1 , φ 2 , φ 3 ) (28b) φ 3 = f 3 (φ 1 , φ 2 , φ 3 , ) , (28c)
where " " denotes the derivative with respect to time τ . Eqs. ( 25) and ( 28) highlight the "slow/fast" nature of the system (recalling that is small). In this paper, we prefer to use the terminology introduced by Gendelman and Bar [START_REF] Gendelman | Bifurcations of self-excitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF] for which the terms fast and slow are replaced by slow and super-slow respectively, whereby the term fast is reserved for the fast component e jωyt . Therefore, [START_REF] Quaranta | Lead-lag damper robustness analysis for helicopter ground resonance[END_REF] and Eqs. (28) consist of one slow complex variable φ 3 and two super-slow complex variables φ 1 and φ 2 . As already mentioned, Eq. ( 25) is the slow-flow of the system written at the time scale t.

Some of steps that have been developed to obtain equivalently Eqs. (25) (or [START_REF] Bab | Lateral vibration attenuation of a rotor under mass eccentricity force using non-linear energy sink[END_REF]) are very hard to justify mathematically. So it will be important to numerically verify that the behaviors of the solutions captured by RedM is also observed with the complete model RefM.

Definition of the Critical Manifold

Solutions of Eqs. (25) (or equivalently [START_REF] Bab | Lateral vibration attenuation of a rotor under mass eccentricity force using non-linear energy sink[END_REF]), can exhibit slow and super-slow epochs characterized by the speed at which the solution advances. The Geometric Singular Perturbation Theory (GSPT) [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF][START_REF] Jones | Geometric singular perturbation theory[END_REF][START_REF] Desroches | Mixed-mode oscillations with multiple time scales[END_REF] is used to describe these dynamics for 0 < 1. To achieve that, stating = 0, Eqs. ( 25) and ( 28) are reduced respectively to the slow subsystem

       φ1 = 0 (29a) φ2 = 0 (29b) φ3 = f 3 (φ 1 , φ 2 , φ 3 , 0) , ( 29c
)
and to the super-slow subsystem

     φ 1 = f 1 (φ 1 , φ 2 , φ 3 ) (30a) φ 2 = f 2 (φ 1 , φ 2 , φ 3 ) (30b) 0 = f 3 (φ 1 , φ 2 , φ 3 , 0) . ( 30c
)
Then, we use the following result of the GSPT: if 0 < 1, the dynamics of the system during slow (resp. super-slow) epoch is given by the dynamic of the slow subsystem (29) (resp. super-slow subsystem [START_REF] Ghazavi | Bifurcation and nonlinear analysis of nonconservative interaction between rotor and blade row[END_REF]).

The so-called Critical Manifold (CM) [START_REF] Jones | Geometric singular perturbation theory[END_REF] associated to Eqs. ( 29) and ( 30) is defined by the algebraic equation (30c) as

CM := (z 1 , z 2 , z 3 ) ∈ C 3 f 3 (z 1 , z 2 , z 3 , 0) = 0 . ( 31
)
The CM is the set of the fixed points of the slow subsystem Eqs. [START_REF] Bab | Vibration mitigation of a rotating beam under external periodic force using a nonlinear energy sink (nes)[END_REF] and it is also the subspace where the trajectories of Eqs. (30) take place.

From the expression of the function f 3 , the CM takes the form

φ 2 = φ 3 F (|φ 3 |) , ∀φ 1 , ( 32
)
where the complex-valued function F of a real positive variable is given by

F (x) = F R (x) + jF I (x) = 1 - α 1 ω 2 δ - 3α 3 ω 2 δ x 2 + j µ ω δ . ( 33
)
Note that the CM only depends on the resonance frequency of the blades and the damping and the cubic spring constant of the NES.

The structure of the CM as defined by Eqs. ( 32) and ( 33) has the same form as the structure of the CM obtained in [START_REF] Bergeot | Analysis of steady-state response regimes of a helicopter ground resonance model including a non-linear energy sink attachment[END_REF] considering only one NES on the fuselage replacing the resonance frequency ω y associated to the fuselage motion by the resonance frequency ω δ associated to the lag motion. Moreover, this structure is similar to the structure of the CM obtained by Starosvetsky and Gendelman [START_REF] Starosvetsky | Strongly modulated response in forced 2dof oscillatory system with essential mass and potential asymmetry[END_REF] 3 studying harmonic forced vibration of one DOF linear system coupled to a NES and also in [START_REF] Gendelman | Asymptotic analysis of passive nonlinear suppression of aeroelastic instabilities of a rigid wing in subsonic flow[END_REF][START_REF] Gendelman | Bifurcations of self-excitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF] studying a nonlinear self-excitated system coupled to a NES. In [START_REF] Starosvetsky | Strongly modulated response in forced 2dof oscillatory system with essential mass and potential asymmetry[END_REF][START_REF] Gendelman | Asymptotic analysis of passive nonlinear suppression of aeroelastic instabilities of a rigid wing in subsonic flow[END_REF][START_REF] Gendelman | Bifurcations of self-excitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF] the studied systems involve only two DOF therefore the CM is defined in C 2 . The distinctive feature of this work is the fact that the system involves three DOF and hence the CM is defined in C 3 , keeping an equivalent form in the (φ 3 , φ 2 )-space.

The CM can be analyzed in the real domain. To achieve this, polar coordinates are introduced

φ 1 = N 1 e jθ 1 , φ 2 = N 2 e jθ 2 m and φ 3 = N 3 e jθ 3 . ( 34
)
and we compute successively the module and the argument of [START_REF] Nayfeh | Applied nonlinear dynamics[END_REF], that lies to

     N 2 2 = H(N 3 ), (35a)
θ 2 = θ 3 + arctan F I (N 3 ) F R (N 3 ) , ( 35b
)
where

H(x) = x 2 F R (x) 2 + F I (x) 2 , ( 36
)
is a real-valued function.

The shape of the CM in the (N 3 , N 2 )-plane depends on the existence of local extrema of the function H. The local extrema of the function H(x) are given by the positive roots of its derivative H (x). An easy calculus shows that two local extrema occur at

N 3,M = 1 3 2 ω 2 δ -α 1 + ω 2 δ -α 1 2 -3µ 2 ω 2 δ α 3 ( 37
)
N 3,m = 1 3 2 ω 2 δ -α 1 - ω 2 δ -α 1 2 -3µ 2 ω 2 δ α 3 , ( 38
)
if the following relation holds

µ < 1 √ 3 ω 2 δ -α 1 ω δ , ( 39
)
and in this case N 3,M < N 3,m . This situation is illustrated Fig. 4(a). Including the third variable φ 1 as N 1 , the CM is a surface in the 3D-space plotted Fig. 4(b).

Stability analysis on Critical Manifold

Stability analysis of the CM can be performed in terms of fixed points of the slow subsystem [START_REF] Bab | Vibration mitigation of a rotating beam under external periodic force using a nonlinear energy sink (nes)[END_REF].

As in [START_REF] Starosvetsky | Strongly modulated response in forced 2dof oscillatory system with essential mass and potential asymmetry[END_REF], it can be shown that the stability condition of the CM is equivalent to

H (N 3 ) > 0. ( 40
)
Hence the stability range of the CM is characterized by the points (N 3,m , N 2,m ) and (N 3,M , N 2,M ) where

N 2,M = H (N 3,M ) and N 2,m = H (N 3,m ) (41)
which are therefore called fold points [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF]. At each fold point it is associated a landing point, (N 3,d , N 2,m ) and (N 3,u , N 2,M ) respectively, where N 3,d and N 3,u are solutions of

H (N 3,d ) = H (N 3,m ) and H (N 3,u ) = H (N 3,M ) (42)
respectively.

A typical Critical Manifold and it stability range are depicted in (N 3 , N 2 )-plane (see Fig. 4(a)) and in the (N 3 , N 2 , N 1 )-space (see Fig. 4(b)). In the (N 3 , N 2 , N 1 )-space, each fold point defines a folded line (L M and L m ) co-linear to the N 1 -axis. The points N 3,d and N 3,u define respectively the line L d and L u which are also co-linear to the N 1 -axis.

Steady-state regimes

The shape and the stability property of the CM (i.e. the existence of folded lines on which the stability of the CM changes) shown in Fig. 4 allows to define at least three steady-state regimes of the RedM model [START_REF] Quaranta | Lead-lag damper robustness analysis for helicopter ground resonance[END_REF] in coherence with the steady-state response regimes described in Sect. 3.2.

To describe the steady regimes of the slow flow we consider situations on which, after an transient regime, the trajectory of the system arrives at a point M 0 on the CM. After that, three steady-state regimes may be considered:

A fixed point of the RedM is reached. These situations corresponds to a periodic solutions of the RefM. Note that to have equivalence between fixed point of the RedM and Periodic Solutions (PR) of the RefM, the RedM must be written using polar coordinates and argument differences must be considered, see Appendix A of [START_REF] Bergeot | Analysis of steady-state response regimes of a helicopter ground resonance model including a non-linear energy sink attachment[END_REF] for more details.

Relaxation oscillations. The S-shape of the CM in the (N 3 , N 2 )-plane suggests also the possible existence of relaxation oscillations [START_REF] Grasman | Asymptotic Methods for Relaxation Oscillations and Applications[END_REF]: from M 0 the system reaches a point S 0 ∈ L M , jumps to a point S 1 ∈ L u and undergoes a super-slow evolution (in the stable domain of the CM) until it reaches a point S 2 ∈ L m . After another jump from S 2 to S 3 ∈ L u and a super-slow evolution (in the stable domain of the CM), the trajectory of the system returns to S 0 ∈ L M . Such scenario of relaxation oscillations for the slow-flow can explain the existence of Strongly Modulated Responses [START_REF] Gendelman | Asymptotic analysis of passive nonlinear suppression of aeroelastic instabilities of a rigid wing in subsonic flow[END_REF][START_REF] Gendelman | Bifurcations of self-excitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF][START_REF] Starosvetsky | Strongly modulated response in forced 2dof oscillatory system with essential mass and potential asymmetry[END_REF] (SMR) for the RefM. Note that if µ > 1

√ 3 ω 2 δ -α 1
ω δ , the S-shape nature of the CM is lost and therefore relaxation oscillations are note possible.

Explosion. Until a first jump the slow-flow evolves the same way as for relaxation oscillations mechanism. However, instead of moving toward a stable fixed point or the folded line L m , the trajectory of the system follows the CM to the infinity. This scenario explains no suppression regime for the RefM.

The existence of one of the three steady-state regimes described above (or an other) depends of the position and the stability of the fixed points of the RedM [START_REF] Quaranta | Lead-lag damper robustness analysis for helicopter ground resonance[END_REF]. Indeed, a stable fixed point of the RedM placed on the stable part of the CM is a necessary condition to obtain PR of the RefM [START_REF] Luongo | Aeroelastic instability analysis of nes-controlled systems via a mixed multiple scale/harmonic balance method[END_REF]. On the other hand, the relaxation oscillations of the RedM (i.e SMR for the RefM) can exist if both folded lines L M and L m have attractive parts. Position and stability of the fixed points of Eq. ( 25) and attractivity (or repulsively) of the folded lines are determined through the study of the super-slow subsystem (30).

Super-slow subsystem analysis

Introducing the CM Eq. ( 32) in Eqs. (30a) and (30b), the super-slow subsystem is written only with respect to the variables φ 1 and φ 3

       dφ 1 dτ = f 1 (φ 1 , φ 3 F (|φ 3 |) , φ 3 ) (43a) d [φ 3 F (|φ 3 |)] dτ = f 2 (φ 1 , φ 3 F (|φ 3 |) , φ 3 ) . ( 43b
)
Using the polar coordinates [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF] and separating real and imaginary parts, Eqs. (43) can be finally reduced (after some calculation steps) to the following form

               g(N 3 ) dN 3 dτ = f N 3 (N 1 , N 3 , δ 31 ) (44a) g(N 3 ) dδ 31 dτ = f δ 31 (N 1 , N 3 , δ 31 ) (44b) dN 1 dτ = f N 1 (N 1 , N 3 , δ 31 ) , (44c)
where

g(x) = H (x) 2 . ( 45
)
From Eqs. (44), it is possible to detect fixed points and folded singularities. The (regular) fixed points of Eqs. (44), {N e 1 , N e 3 , δ e 31 }, are defined by as

f N 3 (N e 1 , N e 3 , δ e 31 ) = 0, (46a)
f δ 31 (N e 1 , N e 3 , δ e 31 ) = 0, (46b)
f N 1 (N e 1 , N e 3 , δ e 31 ) = 0 (46c) g(N e 3 ) = 0. (46d)
If 1, fixed points computed from Eqs. (46) corresponds to fixed points of the RedM model [START_REF] Quaranta | Lead-lag damper robustness analysis for helicopter ground resonance[END_REF]. As usual, stability of the fixed points are found by looking the sign of the eigenvalues real parts of the Jacobian matrix of the vector function

F 2 = (f N 3 /g, f δ 31 /g, f N 1 ) evaluated at {N e 1 , N e 3 , ∆ e 31 }.
From this, a theoretical bifurcation diagram can be obtained showing amplitudes and stability range of the fixed point (resp. periodic regimes) of the RedM (25) (resp. RefM [START_REF] Luongo | Aeroelastic instability analysis of nes-controlled systems via a mixed multiple scale/harmonic balance method[END_REF]) with respect of the detuning term σ.

As a remark, Eqs. (44) have the same form that the equivalent equations obtained when one NES is attached on fuselage (see Eq. ( 68) of [START_REF] Bergeot | Analysis of steady-state response regimes of a helicopter ground resonance model including a non-linear energy sink attachment[END_REF]). In the previous cited paper, detailed analysis Eqs. (44) is proposed. In particular, folded singularities, which are signs of SMR in the system and gives indications about attractivity (or repulsively) of the folded lines are highlighted and computed. This analysis is not recalled here.

Benchmark of theoretical results and discussion

In this section, a bifurcation analysis deduced from the theoretical results discussed in previous section is compared with a bifurcation analysis obtained from numerical integration of the RefM model [START_REF] Luongo | Aeroelastic instability analysis of nes-controlled systems via a mixed multiple scale/harmonic balance method[END_REF] and of the RedM model [START_REF] Quaranta | Lead-lag damper robustness analysis for helicopter ground resonance[END_REF]. The comparison between numerical simulations of RefM and RedM (i.e. the slow-flow) is important because the capacity of the RedM to reproduce the behavior of RefM reflect the validity and quality of all mathematical developments which derive from (even those not presented in this paper). It is important to note that the NES is not subject to an optimization. The parameters used in this section are chosen to illustrate with only two sets of parameters the potential of the NES to modify the response regime after the bifurcation of the trivial solution.

The detuning parameter σ is used as the bifurcation parameter. For each given σ, the solver NDSolve available in Mathematica is used to solve the ordinary differential equations (( 12) and ( 25)) T ). Numerical bifurcation diagrams are then built starting from the smaller value of σ solving the ordinary differential equations with initial conditions as a small perturbation of the trivial solution and next increasing σ and solving the ordinary differential equations with initial conditions corresponding to the steady-state regimes obtained for the previous value of σ.

A first comparison is shown in Fig. 5 where the bifurcation diagrams associated to the three variables N 1 (Fig. 5 First of all, it interesting to note that the responses obtained by numerical integrations of Eqs. [START_REF] Luongo | Aeroelastic instability analysis of nes-controlled systems via a mixed multiple scale/harmonic balance method[END_REF] and Eqs. ( 25) are very close. Moreover, these responses are also very close to the responses predicted by the analytical method when the stability criterion is satisfied i.e for -1.5 ≤ σ ≤ -0.13 and 0.64 ≤ σ ≤ 1.5 (Fig. 5, gray dot markers). We can observe that, starting from σ = -1.5, the transition from complete suppression of GR instability to partial suppression through PR occurs when the trivial solution loses its stability (at σ ≈ -0.6, numerical value given by the theoretical approach). In the same way, the transition from partial suppression through PR to partial suppression through SMR, manifested as a jump of the steady-state amplitudes of the numerical simulations, occurs when the stable fixed points of the RedM lose their stability (at σ ≈ -0.13, numerical value given by the theoretical approach). For -0.13 < σ < 0.64, fixed points of the RedM are unstable and SMR are observed on numerical simulations. Finally, for σ > 0.64 (numerical value given by the theoretical approach), the trivial solution get its stability back and complete suppression is again observed on numerical simulations for both RefM and RedM.

To complete this first comparative study, a trajectory of the RedM (25) when SMR occurs is plotted, Fig. 6, in the (N 3 , N 2 )-plane and compared to the Critical Manifold given by Eq. (35a). The numerical value σ = 0 is used. This simulation illustrates the relaxation oscillations scenario.

These results illustrate how the analytical approach is able to predict different types of steady state regimes and explain the suppression mechanism. However, two comments have to be made on the limitation of the proposed procedure. The first point concerns the existence of the SMR. The proposed procedure allows finding only necessary conditions for the SMR and does not give any sufficient condition as it time-invariant done for example in [START_REF] Starosvetsky | Strongly modulated response in forced 2dof oscillatory system with essential mass and potential asymmetry[END_REF] considering external forced systems problem or for example in [START_REF] Gendelman | Bifurcations of self-excitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF] considering self-excited systems. Indeed, in works in which the superslow flow subsystem can be reduced to a one-dimensional system [START_REF] Domany | Dynamic responses and mitigation of limit cycle oscillations in Van der Pol-Duffing oscillator with nonlinear energy sink[END_REF][START_REF] Gendelman | Asymptotic analysis of passive nonlinear suppression of aeroelastic instabilities of a rigid wing in subsonic flow[END_REF][START_REF] Gendelman | Bifurcations of self-excitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF], global structure of possible response regimes can be deduced from the linear stability analysis of the slow-flow (i.e. the RedM). In our study, the super-slow subsystem (44) is a three-dimensional system. In this case, the linear stability analysis is not sufficient to predict unequivocally the nature of the steady-state response regimes. Indeed, depending of the initial conditions, a stable fixed point of the RedM in the (δ 31 , N 3 , N 1 )-space can never be reached. The second point concerns the co-existence of distinct stable steady-state responses. The nature of the steady-state response regimes (PR or SMR) depends of the stability of the fixed points of the RedM but also of their basin of attraction. In the example of Fig. 5 the basin of attraction of the stable fixed points reduces when σ approaches the transition value σ ≈ -0.13. Consequently the transition from partial suppression through PR to partial suppression through SMR may occur at a smaller value than that which corresponds to the loss of stability of the fixed point of the slow-flow. This is illustrated in Fig. 7. Indeed, in Fig. 7(a), the RefM and RedM are initialized with initial conditions close to the stable steady-regime and the expected responses (periodic regime for the RefM and fixed point for the RedM) are actually observed. In Fig. 7(b), other initial conditions are chosen (see figure caption the exact value) and SMR is observed even if the fixed point of the RedM is stable. These previous observations highlight that situations in which a competition between a stable fixed point and a stable relaxation oscillations cycle of the slow-flow can occur. Similar problems due to the relatively high dimension of the system was already encountered in previous study by the authors [START_REF] Bergeot | Analysis of steady-state response regimes of a helicopter ground resonance model including a non-linear energy sink attachment[END_REF] and also when a NES is used as a passive control device for a nonlinear elastic string, in internal resonance conditions, excited by an external harmonic force [START_REF] Luongo | Nonlinear energy sink to control elastic strings: the internal resonance case[END_REF].

A second comparison is shown in Fig. 8 considering the same set of numerical values except for λ y = 0.3, S d = 1, α 1 = 0.05 and α 3 = 7. For this set of values, the suppression of the GR instability is not satisfied over the whole σ-interval [-3, 3]. As in Fig. 5, starting from σ = -3., the transition from complete suppression of GR instability to partial suppression through PR occurs when the trivial solution loses its stability (at σ ≈ -1.8, numerical value given by the theoretical approach). Partial suppression through PR is satisfied up to σ ≈ -1.45, where the stable fixed points of the RedM lose their stability. For σ > 1.8 (numerical value given by the theoretical approach), the trivial solution get its stability back and complete suppression is again observed. Finally, for -1.45 ≤ σ ≤ 1.8, numerical simulations show responses with partial suppression through SMR and responses with no suppression of GR. Unfortunately, the theoretical bifurcation diagram does not give elements to discriminate SMR, the reason is also the fact that the super-slow subsystem (44) is a three-dimensional system.

As a final remark, to conclude unequivocally about the nature of the steady-state response of the system for a given set of parameters and a given set of initial conditions two informations are missing: the stability range of the relaxation oscillations and the size of the basins of attraction of the stable solutions. This may be subject of future work.

About the comparison between blade and fuselage NES attachments

The configuration using NES on the fuselage (as presented in [START_REF] Bergeot | Analysis of steady-state response regimes of a helicopter ground resonance model including a non-linear energy sink attachment[END_REF]) and configuration using NES attached directly on the blades (as presented here) have the same potential to mitigate GR instability.

This section does not present a rigorous comparison between the two configurations but only few words about it are expounded. The reason is the fact that the efficiency of a NES depends of its design. Therefore, the comparison between NES attached on the fuselage and NES attached on the blades depends more, in our view, on the technical implementation of each solution, which is not the purpose of the paper. However, the blade NES attachments leave more possibilities for engineers for a possible implementation because in this case the small dimensionless parameter depends on both the mass ratio and the squared length ratio between the NES and the blades. Regarding to the fuselage NES attachment configuration in which the small dimensionless parameter depends only on the mass ratio between the NES and the fuselage.

Conclusion

We studied the response of a ground resonance helicopter model including ungrounded NES attached on the blades of the engine. An helicopter model involving fuselage, blades and the cubic NES dynamics was reduced applying a complex multi-blade transformation to the blades and to the NES, defining the Reference Model (RefM).

As usual, in the framework of NES properties exploration and particularly in the context of dynamic instabilities mitigation, four steady-state response regimes were highlighted: complete suppression, partial suppression through periodic response, partial suppression through strongly modulated response and no suppression of the GR instability.

The first regime corresponds to the local stability of the trivial solution of the RefM. To explain all the regimes the slow-flow of the system is determined using complexification-averaging (CA-X). The presence of a small dimensionless parameter related to the moments of inertia of the NES in the slow-flow system implies that it involves one "slow" complex variable and two "super-slow" complex variables. The "super-slow/slow" nature of the system allowed us to use the so-called geometric singular perturbation theory to analyze it.

In particular, the Critical Manifold of the slow-flow was determined as well as its theoretical bifurcation diagram which predicts the possible periodic responses of the RefM et their stability. The S-shape of the Critical Manifold, involving two folded lines, and the associated stability properties provide an analytical tool to explain the existence of three regimes: periodic response regimes, strongly modulated response regimes and no suppression regimes which appear when the trivial solution is unstable. Moreover, the theoretical bifurcation diagram gives indications about the nature the steady-state responses of the system. However, because the super-slow flow subsystem cannot be reduced to a one-dimensional system, the nature the steady-state responses cannot be predicted unequivocally. Indeed, two informations are missing: the stability range of the SMR and the size of the basins of attraction of the stable solutions. This may be subject of future work.

Figure 1 :

 1 Figure 1: Descriptive diagram of the used helicopter system. (a) Overview of the system. (b) View from the top.

Figure 2 :

 2 Figure 2: Descriptive diagram of the NES attachment on the i-th blade. The length [OG ϕ,i ] is called L ϕ .

  3 , b 1 , b 3 and the rotor lateral vibration y. Moreover the modal coordinates b 0 and b 2 are linearly uncoupled from a 1 , a 3 , b 1 , b 3 and the rotor lateral vibration y. In a linear approximation, the components a 0 , a 2 , b 0 and b 2 are neglected. Recalling that a 3 = a * (a) (b)

Figure 3 :

 3 Figure 3: Evolution of imaginary (a) and real (b) parts of the eigenvalues of the vector function F 1 evaluated at the trivial equilibrium. Parameters used: λ y = 0.2, ω y = 1, S d = 1, T d = 0.1, S c = 3, T c = 1, λ = 0.2, ω δ = 2, µ = 0.2, α 1 = 0.5 and = 0.1.

Figure 4 :

 4 Figure 4: Critical Manifold (CM). Following parameters are used: ω δ = 2, α 1 = 0.5, α 3 = 2 and µ = 0.1. (a) In the (N 3 , N 2 )-plane and (b) In the (N 3 , N 2 , N 1 )-space.

Figure 5 :

 5 Figure 5: Comparison between theoretical bifurcation diagram, obtained from Eq. (46) (gray dots stable fixed points and black dots for unstable fixed points) and maximum steady-state amplitudes obtained from numerical simulations of the RefM (12) (blue squares) and from numerical simulations of the RedM model (25) (red circles). Bifurcation diagram and maximum steady-state amplitudes are plotted for the variables (a) N 1 , (b) N 2 and (c) N 3 as a function of σ. When the trivial solution is unstable, parameter range in which PR or SMR are observed are specified on the figure. The parameters used are: λ y = 0.3, ω y = 1, λ δ = 0.2, ω δ = 2, S d = 1, T d = 1, S c = 2, T c = 2, µ = 0.1, α 1 = 0.05, α 3 = 7 and = 0.01.

Figure 6 :

 6 Figure 6: Comparison between the trajectory of the simulated RedM (25) in the plane (N 3 , N 2 ) and the Critical Manifold (35a). Parameters and initial conditions are the same as for Fig. 5 with σ = 0.

  (a)), N 2 (Fig. 5(b)) and N 3 (Fig. 5(c)) are plotted for σ ∈ [-1.5, 1.5] (see Fig. 5 caption to know the used parameters). The chosen set of parameters (see Fig. 5 caption) has been used to illustrate the suppression of the GR instability over a whole σ-interval [-1.5, 1.5]. The corresponding numerical values of the parameters characterizing the RefM model (12) are deduced from (13), they are not reported here.

Figure 7 :

 7 Figure 7: Comparison between numerical integration of the RefM (12) (gray line) and of the RedM model (25) (black line). The parameters are the same as for Fig. 5 with σ = -0.3 and two set of initial condition are used: (a) the systems are initialized close to its stable steady-state regime (fixed point for the RedM and PR for the RefM), (b) the initial conditions are: y(0) = 0.3 for the RefM and φ 1 (0) = j0.3 for the RedM, all others initial conditions are equal to zero.

Figure 8 :

 8 Figure 8: Comparison between theoretical bifurcation diagram, obtained from Eq. (46) (gray dots stable fixed points and black dots for unstable fixed points) and maximum steady-state amplitudes obtained from numerical simulations of the RefM (12) (blue squares) and from numerical simulations of the RedM model (25) (red circles). Bifurcation diagram and maximum steady-state amplitudes are plotted for the variables (a) N 1 , (b) N 2 and (c) N 3 as a function of σ. When the trivial solution is unstable, parameter range in which PR SMR or Explosion are observed are specified on the figure. The parameters used are: λ y = 0.075, ω y = 1, λ δ = 0.2, ω δ = 2, S d = 4, T d = 1, S c = 2, T c = 2, µ = 0.1, α 1 = 0, α 3 = 10 and = 0.01.

and b 3 = b * 1 , we therefore consider only the equations of motion associated to y, a 1 and b 1 . For convenience, variables a 1 and b 1 are now noted simply a and b.

This can be shown for example by computing the power spectrum of the steady part of the signal.

The complexification is not necessary for the variable A(t) and B(t) because they are already complex variables.

In the reference[START_REF] Starosvetsky | Strongly modulated response in forced 2dof oscillatory system with essential mass and potential asymmetry[END_REF] the critical manifolds is called Slow Invariant Manifold.

Acknowledgements

This research was done within the framework of the industrial chair "Dynamique des Systèmes Mécaniques Complexes (Dynamics of complex mechanical systems)" financed by the Airbus Group Foundation.


Ganesh Kumar

Hari Shankar

Lal Das

Bertrand Tondu

Florent Forget

Jérôme Manhes

Olivier Stasse

Philippe Souères

P Soueres

Controlling a multi-joint arm actuated by pneumatic muscles with quasi-DDP optimal control

I. INTRODUCTION

Pneumatic actuators are inherently compliant and have very high power to weight ratio. These factors motivate researchers to use pneumatic actuation for exoskeleton, prosthesis, rehabilitation or even in walking robot and human robot interaction. The goal of this paper is to use an arm actuated by McKibben muscles to perform explosive tasks like hammering or kicking balls. Recent works [START_REF] Kumar | Fast, strong and compliant pneumatic actuation for dexterous tendon-driven hands[END_REF] advocates other kind of actuation system which are easier to control but come with a high price tag. Our goal was to demonstrate here that with minimal modification of an existing platform [START_REF] Tondu | A seven-degrees-of-freedom robot-arm driven by pneumatic artificial muscles for humanoid robots[END_REF], it is possible to perform such kind of explosive motions, and not fine manipulation as in [START_REF] Kumar | Fast, strong and compliant pneumatic actuation for dexterous tendon-driven hands[END_REF]. The Mckibben artificial muscle is known for its nonlinearities and hence pose a great control challenge. These non-linearities are mainly due to hysteresis, saturation and internal friction between fabrics. So far these challenges were dealt with traditional controllers like high gain PID controller, sliding mode controller for position control [START_REF] Tondu | Second order sliding mode control for an anthropomorphic robot-arm driven with pneumatic artificial muscles[END_REF]. These methods usually lead to stiff system dynamics with higher impedances. Recent developments in pneumatic actuation control [START_REF] Tassa | Modeling and identification of pneumatic actuators[END_REF], [START_REF] Kumar | Real-time behaviour synthesis for dynamic hand-manipulation[END_REF] are use Model Predictive Control (MPC). For instance in [START_REF] Andrikopoulos | Pneumatic artificial muscles: A switching model predictive control approach[END_REF], the control scheme is based on a switching Piece Wise Affine (PWA) system model approximation. The method is able to capture the high non-linearities of the Pneumatic Artificial Muscle (PAM). In [START_REF] Tassa | Modeling and identification of pneumatic actuators[END_REF] a linear formulation of the actuation system is proposed to simplify the algorithm implementation and makes it real-time.

The contribution of this paper is two fold.

-First it proposes a new model which provides a good compromise between accuracy and simplicity. 1 are with CNRS, LAAS, 7 avenue du colonel Roche, F-31400 Toulouse, France 2 are also with Univ de Toulouse, INSA, F-31400 Toulouse, France gkharish,btondu,ostasse,psoueres@laas.fr -The second contribution is to show that, on this model an iLQR control scheme can be used to generate with good precision positioning task and explosive movements. This is achieved with the introduction of an empirical model of pressure generation from the Intensity-Pressure converter (I/P). Each of the seven Degrees of Freedom (DoFs) of the manipulator arm of LAAS-CNRS is actuated by a pair of agonist-antagonist Mckibben muscles. The capabilities of the non-linear iLQR control to execute some simple tasks of reaching a point, maximizing the velocity of joint and end effector are demonstrated on this platform.

The paper is organised as follows: The dynamical model of joints actuated by Mckibben muscles is presented in section II. Section III describes the model of the anthropomorphic robot arm. The optimal control formulation and brief introduction to iterative Linear Quadratic Regulator(iLQR) is presented in section IV. Section V presents simulations and experiments.

II. DYNAMIC MODEL OF A JOINT ACTUATED BY MCKIBBEN MUSCLES

There have been several attempts to model the McKibben artificial muscle [START_REF] Tassa | Modeling and identification of pneumatic actuators[END_REF], [START_REF] Tondu | Modelling of the mckibben artificial muscle: A review[END_REF], [START_REF] Kothera | Experimental characterization and static modeling of mckibben actuators[END_REF], [START_REF] Tondu | Mckibben artificial muscle can be in accordance with hill skeletal muscle model[END_REF], [START_REF] Daerden | The concept and design of pleated pneumatic artificial muscles[END_REF]. The difficulties for getting an accurate model are due to a combination of complex phenomena during static and dynamic contraction: shape changing at the muscle tip which lose their initial cylindrical shape for a conical one, mobility and flexibility of the braided sleeve, elasticity of the inner rubber tube, exotic friction in the textile braided sleeve without forgetting, for the pneumatic version of the McKibben muscle, dynamic fluidic phenomena resulting from the artificial muscle volume variation during contraction. When we want to include such a physical model into the closed-loop control of a McKibben muscle actuator, a compromise must be found between the complexity of an accurate model and the time required for its computation. The need for an efficient but not too complex dynamic muscle model is all the more important in the case of artificial muscles for robot arms. Indeed the preservation of the actuator compliance imposes a joint direct-drive mode and so the consideration of a dynamic robot model, linked to the actuator model. The originality of our approach, described in this section, consists in including an original model of real pressure variation inside a simple McKibben muscle dynamic model before proposing a multi-variable pressure control of the actuator made of two similar antagonist muscles.

A. Muscles side dynamics

The following model was considered for relating the dynamic contraction force F of the artificial muscle to its control pressure P and contraction ratio. Where,

F(ε, P) = (πr 2 0 )P[a(1 -kε) 2 -b] (1)
ε = l o -l l o a = 3 tan 2 α o b = 1 sin 2 (α o )
l being the current muscle length and l o , r o , α o are the muscle geometric parameters corresponding to its initial length, radius and braid angle respectively (see [START_REF] Tondu | Modelling of the mckibben artificial muscle: A review[END_REF] for more details). The parameter k, slightly greater than 1, is empirically chosen for taking into account the conic shape at muscle tips. The viscous coefficient is also experimentally estimated for expressing the naturally damped behavior of the artificial muscle contraction due to the complex kinetic friction inside the braided sleeve. Although hysteresis phenomenon due to dry friction is not included in this model, we think it captures most of the static and dynamic behaviour of the artificial muscle.

B. Pressure side dynamics

Control pressure in a muscle is provided by an Intensity-Pressure (I/P) converter which translates a current value into a desired pressure value that has to be fed to the muscle. In literature, see [START_REF] Kelasidi | A survey on pneumatic muscle actuators modeling[END_REF] for a survey, there exist several models of pressure generated in a muscle in terms of the mass of air injected, volume of the muscle and some parameters of servo valves or I/P converter (for an example see [START_REF] Vanderborght | Torque and compliance control of the pneumatic artificial muscles in the biped "lucy[END_REF]. These models are highly non-linear and eventually become too cumbersome for devising control strategies. We propose a somewhat simpler empirical model to cover the dynamics of the pressure generation from the I/P converter. It is derived from SAMSON I/P converters reported in [START_REF] Boitier | Design and control of a 2-dof scara-type robot actuated by pneumatic mckibben artificial muscles[END_REF]. Following this model the instantaneous pressure p inside a muscle is represented by a damped second order differential function with a control pressure as input as follows:

p + 2w n ṗ + w 2 n p = w 2 n p des (2)
Note that the input of the intensity converter is a current value which is then scaled to the corresponding pressure control input P des in the pressure unit. The above equation is also non-linear as the natural frequency w n has been empirically identified as a function of the volume of the muscle, V .

w n = 2π f v (1/V ) (3)
f v is the empirically found parameter. The volume of a McKibben artificial muscle can be approached by the following cylindrical volume: V = πr 2 l, where r and l are the current radius and length of the artificial muscle. By using the relationship Eq.( 4) as reported in [START_REF] Tondu | Theory of an artificial pneumatic muscle and application to the modelling of mckibben artificial muscle[END_REF], we deduce an expression for the volume of the muscle.

r/r 0 = (1/sin(α o )) (1 -ε) 2 ) (4) V = πr 2 o l sin 2 (α o ) (1 -cos 2 (α 0 )(1 -ε) 2 ) (5)
For a second order system, the rise time is inversely proportional to its natural frequency and in the case of our pressure dynamic model, the natural frequency is inversely proportional to the volume. It implies that the bigger the muscle, the larger the rise time. Fig. 2 Fig. 2. Natural frequency and volume of the muscle with the joint angle of operation for a joint is small i.e range of θ is small, then the corresponding changes in the volume of the muscles as θ varies, would be small at that joint. In such case, a constant value for w n could be chosen for the muscles corresponding to their mean volumes. Here, another interesting point to notice is that for different muscles, pressure dynamics of a larger volume muscle will be slower. The above two points help in trading off between the operating range of the joint and the corresponding easiness of control strategy.

C. Agonistic-antagonistic joint actuator

Following the human arm model, a pair of artificial muscles can be set up in antagonistic fashion to drive a chained wheel of radius R. According to Fig. 3, the resulting actuator torque T can be written as follows:

T = R[F 1 (ε 1 , P 1 ) -F 2 (ε 2 , P 2 )] (6)
where F 1 and F 2 are the forces of muscles 1 and 2 respectively, defining the antagonist muscle pair. As illustrated in Fig. 3, the relationship between θ , the actuator angle and the contraction ratios for each muscle can be expressed as follows:

ε 1 = ε 10 + Rθ l 0 and ε 2 = ε 20 -Rθ l 0
, where ε 10 and ε 20 are the initial contraction ratios for muscle 1 and muscle 2 respectively corresponding to the zero-angular position. Moreover, we propose to specify the pressures in muscle 1 and muscle 2 as follows: P 1 = P 10 + ∆P 1 and P 2 = P 20 + ∆P 2 where P 10 and P 20 are respectively the initial pressure in the agonist-antagonist muscle and related to zero positioning of the joint. ∆P 1 and ∆P 2 , are the control pressures.

Using Eq.( 1), and neglecting the terms ε 2 1 and ε 2 2 , the torque at the joint can be expressed as

T = K 1 ∆P 1 -K 2 ∆P 2 -K 3 (P 1 + P 2 )θ -K 4 θ + K 5 , (7) θ
K 1 = (πr 2 o )R[a(1 -2kε 10 ) -b], K 2 = (πr 2 o )R[a(1 -2kε 20 ) -b], K 3 = 2(πr 2 o )R 2 ka/l o , K 4 = (πr 2 o )R 2 f v /l o , K 5 = (πr 2 o )R[(a -b)(P 10 -P 20 ) -2ka(P 10 ε 10 -P 20 ε 20 )].
If P 10 = P 20 and ε 10 = ε 20 , K 5 is equal to zero. The details of our robot-arm including the dimensions of muscles are given in [START_REF] Tondu | A seven-degrees-of-freedom robot-arm driven by pneumatic artificial muscles for humanoid robots[END_REF]. The antagonist muscle actuator is now considered as a MIMO-system whose inputs are the control pressures ∆P 1 , ∆P 2 , and outputs are both the θ and the actuator stiffness which is defined as the instantaneous ratio between the current torque variation and the current angular position variation (see Eq.( 9)). When no gravity effect is considered, the static equilibrium position of the actuator can be directly derived from Eq.( 7) with zero angular velocity.

θ equ = (K 1 ∆P 1 -K 2 ∆P 2 + K 5 )/(K 3 (P 1 + P 2 )) (8)
With associated stiffness σ equ expressed as

σ equ = - ∂ T ∂ θ = K 3 (P 1 + P 2 ) (9)
From Eq.( 8) and Eq.( 9), it is possible to remark that the equilibrium position can be changed while keeping the same stiffness by modulating ∆P 1 and ∆P 2 with a constant ∆P 1 +∆P 2 .

In the case of a symmetrical pressure variation in both muscles: ∆P 1 = -∆P 2 = ∆P , the actuator becomes a SISO-system whose corresponding torque T SISO , is now given by the following relationship:

T = (K 1 + K 2 )∆P -K 3 (P 10 + P 20 )θ -K 4 θ + K 5 ( 10
)
where stiffness at equilibrium position is now constant and equal to K 3 (P 10 + P 20 ). When the joint angle θ varies, and assuming that the actuator chain is inextensible, the contraction ratio of each muscle is known and, consequently, the pressure dynamic model, proposed in II-B, can be applied to each muscle. It is important to note that, in the case of relatively small joint muscles, the effect of volume variation with joint angle on pressure dynamics could be negligible. If the range of operation is small, a constant w n could be chosen for that joint.

As previously noted, pressure dynamics of a larger volume muscle will be slower which is, especially, the case of the muscles at the shoulder joint compared to the muscles at the elbow joint.

III. ROBOT DYNAMICS

This section presents the robot model formally including the rigid body model of the robot with its pressure dynamics. Let us consider a n degrees of freedom robot with generalized joint angle coordinates q ∈ R n . Each joint is actuated by 2 pneumatic muscles, so there will be 2n muscles, each one with a pressure P i ∈ R 2n . The robot dynamics can be represented in standard form as below: M(q) q +C(q, q) + G(q) = T (q, P)

P + 2C p Ṗ + G p P = G p P des (11)
M(q) ∈ R n×n is the mass inertia matrix of the robot, C(q, q) ∈ R n×1 is the coriolis and centrifugal terms, G(q) ∈ R n is the gravity related terms, T (q, P) ∈ R n is the torque generated by pneumatic muscles as described in Eq. [START_REF] Andrikopoulos | Pneumatic artificial muscles: A switching model predictive control approach[END_REF].

C p = diag[w n 1 , w n 2 , ...w n 2n ] and G p = diag[w 2 n 1 , w 2 n 2 , ...w 2 n 2n
] are the collection of coefficients of the pressure dynamics of each muscle and w n i is the natural frequency of the pressure dynamics at the i th muscle. P = [P 1 , P 2 .....P 2n ] T ∈ R 2n is the vector of current muscle pressure and P des = [P des 1 , P des 2 .....P des 2n ] T ∈ R 2n , is the vector of desired pressure P des i , the control input in pressure unit given to the I/P converter of the i th muscle. Any reference trajectory in q given to the robot has to respect the pressure dynamics whose bandwidth depends on w n . As discussed in the previous section, if the joint is actuated by larger muscles (having smaller bandwidth) and higher rise time, rapidly varying trajectory would be difficult to track and pose a limiting factor in performing explosive tasks. Limits to the operating range of the joints are obviously decided by the maximum contraction of the spanning muscles. In addition to that, w n could be reasonably considered as constant in a limited operating range, removing in this way the non-linearity in pressure dynamics. Apart from this, due to hardware limitation of I/P converter, the control input P des is bounded as follows:

U = P des ∈ R 2n : P des ∈ [P min , P max ] (13)
Where, P min , P max are the lower and upper bound on control input P des .

IV. OPTIMAL CONTROL FORMULATION

This section presents the optimal computational framework used to find the control sequences to perform a desired task. The optimal control problem is the minimization or maximization of a performance criterion with respect to the control under a set of constraints that arise from the physical limitation of the control action and from the plant dynamics.

A. State space representation

Let us represent the dynamics stated in Eq.( 11) and Eq.( 12) in state space form considering the state vector as x = [q, q, P, Ṗ] T .

ẋ = f (x, u) =    q M -1 (-C(q, q) -G(q) + T (q, P)) Ṗ -C p Ṗ -G p P + G p u    (14)
where, f is the non-linear function given by Eq.( 14) in state x and control u that gathers Eq.( 11) and Eq.( 12). In the present work, we consider the constraints on the state in the optimal control formulation which is discussed in the following subsection.

B. Treatment of Constraints

There exists a mechanical limit to each degrees of freedom of the arm. To have safe operation, we have introduced these limits as constraints on the state space inside the cost function of the optimal control formulation. The chosen cost function for the state space constraints are expressed in the following equations:

max = 1 -λ (x max -x) min = 1 -λ (x -x min ) C s = e λ * max + e λ * min (15)
where x min , x max are lower and upper limits on the state and λ is a constant. The above consideration for the cost function will ensure that the cost near the limits will be very high as evident in the Fig: 4 and hence the optimal solution will keep the system within the operating limits. Also, the control action In order to handle the constraints on the control, the control limited Differential Dynamic Programming (DDP) approach reported in [START_REF] Tassa | Control-limited differential dynamic programming[END_REF] is used. The control problem formulation is then expressed as determining an open-loop control input u = u(t, x) ∈ U which can minimize or maximize a cost function along a given time interval t ∈ [0, T ] and with initial state x(0) = x 0 . The most generic expression for the cost function can be written as follows:

J(x 0 ) = C f +C r +C s (16)
where

C f = h(x(T )) is the final cost, C r = T 0 c(x(t), u(t, x(t))
)dt is the integral of the running cost c(x, u) which encapsulates the task objectives and C s is the cost value imposed by the constraints. For a non-linear dynamics Eq.( 14) and non-quadratic cost Eq.( 16), optimal control solutions can be obtained using full DDP. However, as DDP is computationally expensive, an iterative LQR (iLQR) approach is considered [START_REF] Li | Iterative linearisation methods for approximately optimal control and estimation of non-linear stochastic system[END_REF]. The iLQR method is relying on linearizing the dynamics and approximating the cost function to quadratic form along the x trajectory. This control approach is briefly summarized in the next section.

C. Iterative Linear Quadratic Regulator (iLQR)

iLQR is initialized with a nominal control sequence and the corresponding state trajectory (x 0 , u 0 ). the dynamical system is then linearized as in Eq.( 17) and the cost function is approximated by the quadratic form Eq.( 18) and a local LQR problem is then solved. Using this solution, the states and the control sequence are improved iteratively.

δ x = Aδ x + Bδ u (17)
∆J = h T x δ x(T ) + δ x T (T )h xx δ x(T ) + T 0 c T x δ x + c T u δ u + δ x T c xx δ x + δ x T c xu δ u + δ u T c uu δ u (18)
where

A = ∂ f ∂ x and B = ∂ f ∂ u .
In Eq.( 17) and Eq.( 18) subscript x and u indicate that the function is partially derivated with respect to x and u. At every iteration, Eq.( 17) and Eq.( 18) are solved and (δ x, δ u) are deduced from the resolution of a modified Ricatti-type system. Then the new improved sequence is generated by x ← x + δ x and u ← u + δ u. When ∆J ≈ 0, the iLQR converges and gives an optimal control sequence u * ∈ U and the corresponding optimal state trajectory x * .

V. SIMULATION AND EXPERIMENTAL RESULTS

The experimental set up considered here is the manipulator of LAAS CNRS which is an anthropomorphic arm of seven degrees of freedom (DoF), where each joint is pneumatically actuated by a pair of Mckibben muscles. In the experiment presented in this paper, only two joints are controlled. These joints are the flexion θ 1 at the shoulder and the flexion θ 2 at the elbow. (see Fig. 5). So, for the experiments, the robot can be viewed as a 2 DoF manipulator with state defined by x = [θ 1 , θ1 , θ 2 , θ2 ] T . Each muscle is pressurized by an I/P converter which converts a current command to a reference pressure value. The objective of the experiments is to evaluate the performance of the iLQR control to achieve the following tasks with our multi-link pneumatic robot. 1) Final position control: We aim to compare the quality of the positioning control in the presence of different loads with the iLQR approach and with the feed-forward proportional control. 2) Capability to execute explosive movements by maximizing the link speed at a given time.For this second task we compare simulation results with the results of experiment executed with the real robot. For each task, we analyze the stiffness modulation and we show that the optimal control approach enhances the explosive motion capabilities by simultaneous modulation of position and stiffness. In order to do that, simulations and experiments are done in two cases. 1) Case-I: When the sum of the pressure in agonistantagonist pair of muscles at each joint is kept constant, i.e. P 1 + P 2 = Constant. 1) Case-II: When P 1 and P 2 are left independent. In this case, there will be two inputs for each joint.

Shoulder Flexion Muscle Pair

Elbow Flexion Muscle Pair

A. Task 1: Position control

The manipulator is given a final position with a load mass of m l = 0.1kg at the end effector. iLQR is used here to find the optimal path to reach the final goal (θ 1 , θ 2 ) = (28.8, 57.6)degrees. The same task is repeated with different load masses of m l = 0.1kg and m l = 1.0kg which is shown in Fig. 6.

The following cost function is considered for this task.

C f = Q f (x re f (T ) -x(t)) 2 , C r = Q T 0 ((x re f (T ) -x(t)) 2 + u 2 (t))dt, (19)
where, x re f is the final position for the two joints. The iLQR control uses Eq.( 14) and Eq.( 19) to solve for optimal control sequence u * . This optimal control sequence in forward application to Eq.( 14) will yield the needed trajectory. To compare the effectiveness of the iLQR approach, the task is executed using a Proportional-Integral controller with a feedforward term. The feed-forward term gives a desired pressure for each muscle at the joint needed to maintain a desired joint angle. Thus, the control action of the feed-forward PI controller can be defined as follows:

u(t) = P f eed + K p (e(t)) + K i t 0 (e(t))dt (20)
where, e(t) = x re fx(t), K p and K i are proportional and integral gains respectively. The simulation results are shown in Fig. 7. Response of feed-forward PI controller is shown in black dashed line is compared with the responses of optimal position control in Case-I when P 1 +P 2 = Constant (blue lines) and in Case-II when P 1 and P 2 are independent (magenta lines)(See Fig. 7).

From the position plots in Fig. 7, it appears that the iLQR gives a good compromise between keeping the stiffness low and minimizing the oscillations which results into a smooth motion. However, the feed-forward PI controller The plots compare the response of the joints and the stiffness profile when the end effector is loaded with mass 0.1kg (solid line) and 1.0kg (dashed line). For the heavier load, optimal control gives a stiffer stiffness profile. makes joint 1 very stiff and joint 2 very flexible leading to an overshoot and oscillations.

As discussed in section II-C, the stiffness can be adjusted for the same equilibrium position by changing the sum of pressures in the agonist-antagonist pair of muscles. As the sum increases, stiffness at the equilibrium increases which is evident at each joint (blue and brown lines in stiffness plots).

In order to find the most compliant position trajectory, i.e to find the best sum of pressure in the agonist-antagonist pair, P 1 , P 2 are left independent and Eq.( 7) has two independent inputs for each joint. In this case, the optimal solution gives the similar optimal position trajectory but with better stiffness profile (magenta line in stiffness plot in Fig. 7).

B. Task 2: Maximizing the link speed

The objective is to execute some explosive motions with the aim to perform in the future tasks such as ball throwing, kicking or hammering a nail. Such motions would require either maximizing the joint/link speed or end-effector speed [START_REF] Garabini | Optimality principles in variable stiffness control: Vsa hammer[END_REF], [START_REF] Braun | Robots driven by compliant actuators: Optimal control under actuation constraints[END_REF]. Here, two subtasks are presented: 1) Maximizing the angular speed of the elbow joint. The task is first simulated and then executed by the real robot. A comparison between simulation and the experimental results are shown in Fig. 8. 2) Maximizing the end effector speed: The objective of this task is to achieve maximum end-effector speed at terminal time T . Simulations are done for both Case-I (P1 1 + P 2 = Constant) and Case-II (P 1 and P 2 are independent). 1) Maximizing the angular speed: For maximizing the elbow joint's angular speed at final time T = 1, the task requires only the terminal cost, C f . But to minimize the control effort, a running cost C r involving only control pressures is used.

C f = -Q f (x 4 (T )) 2 , C r = Q u T 0 u 2 (t)dt, (21)
where, x 4 = θ2 is the angular velocity of the elbow joint which constitutes the fourth state variable. Q f and Q u are weights for the terminal cost and the running cost. The simulation results are shown in Fig. 8 Plots on the top compare the robot response (solid lines) with simulated response (dashed lines) at joint 1 (blue line) and joint 2 (red line). Bottom plots compare the velocity response of the robot and the simulation.

In simulation, the maximum joint velocity is reached at terminal time which is around -15 rad/sec. The robot response and the simulation response in position and speed show a good match. However, some discrepancies can be observed. Apart from some modeling error, this difference is probably due to the fact that the constraints on the state vector are not considered in the optimal control problem formulation. This is an improvement that we plan to do in future work. Snapshots of the video of the experiment are presented in Fig. 12 2) Maximizing the end effector speed: The objective is to maximize the end-effector speed of the manipulator at terminal time. The cost function, thus, comprises the terminal cost C f involving the end-effector velocity and the running cost C r involving only control efforts.

C f = -Q f ( f vel (x(T )) 2 , C r = Q u T 0 u 2 (t)dt, (22)
where, f vel (x(t)) is the function which computes end effector velocity from the robot kinematics. Simulations are done for the following two cases. In Case-I, the sum of pressures in the agonist-antagonistic pair is kept constant at the joint 1 to 3 bar and at the joint 2 to 4 bar. In this second simulation (Case-II), Comparison of the stiffness profile of each joint (joint 1: blue lines, Joint 2: red lines) in Case-I when P 1 + P 2 = Const (dashed lines ) and Case-II when P 1 , P 2 are independent (solid lines). It is evident from the plots that the motion in Case-II is indeed slightly more compliant than the motion in Case-I.

we let P 1 , P 2 independent.

Comparing the velocity plots and the corresponding stiffness plots in Fig. 11, it can be seen that the maximum end-effector velocity has increased from 3.8m/s in Case-I to 4.1m/s in Case-II. This is due to the fact that in Case-II, the optimal control has the freedom to modulate the compliance with two inputs at each joint and hence it was able to achieve higher speed.

The hardware set up used for the experimentation is the robot mentioned in [START_REF] Tondu | A seven-degrees-of-freedom robot-arm driven by pneumatic artificial muscles for humanoid robots[END_REF]

VI. CONCLUSION AND FUTURE WORK

We considered the problem of modeling and controlling a robot manipulator actuated by the pneumatic Mckibben muscles in an optimal control frame-work. Our first contribution is to propose a model which encapsulates the pressure dynamics in an efficient way. It is done by using a second order differential function whose bandwidth is dependent on the instantaneous volume of the muscles. Even though, the model does not include the effect of dry friction and hysteresis, it covers most of the static and the dynamic behavior. Using this model, along with the robot manipulator model, our second contribution is to show that the implementation of iLQR allows to perform efficient position control but also preserves the inherent compliance of the Mckibben muscles with respect to the more conventional control approaches. The optimal control formulation with this new model is used to perform various tasks like positioning control under different loads and maximizing the link speed. The results were reported in simulation first and then some of these results were validated on the hardware platform. In future work, we plan to extend the application on the real robot in order to perform explosive movements such as throwing a ball or hammering a nail. This will be done by embedding the state constraints in the optimal control formulation. A detailed analysis of stiffness variation and possibilities of stiffness control will then be possible.

Fig. 1 .

 1 Fig. 1. Muscle geometric parameters: l 0 initial length, r 0 radius and α 0 braid angle

2 ε 10 +Fig. 3 .

 2103 Fig. 3. Pulley-chain driven antagonist muscle actuator made of two identical McKibben artificial muscles. It has the possibility to adapt initial torque at initial zero θ -position by means of P 10 , P 20 , ε 10 and ε 20

Fig. 4 .

 4 Fig. 4. Cost function of the state constraints has to be admissible, i.e u ∈ U = [P min , P max ].In order to handle the constraints on the control, the control limited Differential Dynamic Programming (DDP) approach reported in[START_REF] Tassa | Control-limited differential dynamic programming[END_REF] is used. The control problem formulation is then expressed as determining an open-loop control input u = u(t, x) ∈ U which can minimize or maximize a cost function along a given time interval t ∈ [0, T ] and with initial state x(0) = x 0 . The most generic expression for the cost function can be written as follows:J(x 0 ) = C f +C r +C s[START_REF] Li | Iterative linearisation methods for approximately optimal control and estimation of non-linear stochastic system[END_REF]

Fig. 5 .

 5 Fig. 5. Picture of the robot-arm showing the two pairs of muscles, in shoulder and elbow flexion.

Fig. 6 .

 6 Fig.6. Figures on the left side represents the position response and stiffness variation at joint 1 and the figures on the right side represents the same for joint 2 in order to reach the target (θ 1 , θ 2 ) = (28.8, 57.6)degrees. The plots compare the response of the joints and the stiffness profile when the end effector is loaded with mass 0.1kg (solid line) and 1.0kg (dashed line). For the heavier load, optimal control gives a stiffer stiffness profile.

Fig. 7 .

 7 Fig. 7. Figures on the left side represents the position response and stiffness variation at joint 1 and the figures on the right side represents the same for joint 2.

Fig. 8 .

 8 Fig.8.Plots on the top compare the robot response (solid lines) with simulated response (dashed lines) at joint 1 (blue line) and joint 2 (red line). Bottom plots compare the velocity response of the robot and the simulation.

Fig. 9 .Fig. 10 .Joint 1 (

 9101 Fig. 9. Case-I: P 1 + P 2 = Const. Plots on the top show the simulated position response of joint 1 (blue line) and Joint 2 (black line). Control effort is shown in the bottom plots. The control inputs, Input 1 and Input 2 are applied on joint 1 and joint 2 respectively.

Fig. 12 .

 12 Fig. 12. The pneumatic manipulator arm executing the task of maximizing the angular speed of the elbow.

  whose control modules have been upgraded recently. The important components of the current set up are I/P converters, encoders, NI data acquisition devices and the development computer running real time control software. 1) I/P converter : It is a Samson I/P 6111 manufactured by Samson Corporation, Frankfurt. This I/P is rapid and produces output pressure in range 0 to 5 bar. The bandwidth of the I/P is volume dependent. There are 7 joints and a gripper. Each joint and the gripper are actuated by a pair of Mckibben muscles and each muscle is controlled by one I/P so, there are 16 IP converters. 2) Encoders/Potentiometers: There are seven potentiometers, one at each Dof to get the joint angle data. 3) The data Acquisition device: CompactRIO from National Instruments is used for this purpose. This is a reconfigurable and embeddable chassis, integrating an intelligent real-time controller. It has NI9205 module which is an Analog-to-Digital-Converter (ADC) with up to 32 channels, selectable acquisition range (200mV, 1V, 5V and 0V ) and a max sampling frequency of 250kS/s. This device samples the seven potentiometers dedicated to the positioning of the arm. There is a NI9265 card which has three modules of Digitalto-Analog-Converters (DAC). They have the particularity to provide analog current outputs. Each module is able to manage up to 4 channels in the 0 to 20mA range at 100kS/s. 4) The development computer: The development computer runs the high level control algorithm for the robot. It communicates via the NI module using standard UDP protocol at a frequency of 1kHz in real time. The computer has linux 14.04 with Xenomai real time kernel running in parallel. For the present experimentation, control software are running at 200Hz frequency.

VII. ACKNOWLEDGEMENTS

This work was partially supported by the European FP7 Project EUROC under Grant Agreement 608849 and the ERC Advanced Grant Actanthrope. We thank Jerome Manhes for his support in hardware development, Florent Forget and Kevin Giraud Esclasse for their help in software development.