Stefano Lenci

Giuseppe Rega

Optimal control of the homoclinic bifurcation in buckled beams: Infinite dimensional vs reduced order modeling

Keywords: Buckled beams, infinite dimensional system, homoclinic bifurcation, Holmes & Marsden theorem, nonlinear normal modes, reduced order models, optimal control of chaos

A method for controlling nonlinear dynamics and chaos is applied to the infinite dimensional dynamics of a buckled beam subjected to a generic space varying time-periodic transversal excitation. The homoclinic bifurcation of the hilltop saddle is identified as the undesired dynamical event, because it triggers, e.g., cross-well scattered (possibly chaotic) dynamics. Its elimination is then pursued by a control strategy which consists in choosing the best spatial and temporal shape of the excitation permitting the maximum shift of the homoclinic bifurcation threshold in the excitation amplitude-frequency parameters space.

The homoclinic bifurcation is detected by the Holmes & Marsden's theorem [1981] constituting a generalization of the classical Melnikov's theory. Two classes of boundary conditions are identified: for the first, the Melnikov function is exactly the same as obtained with the reduced order models, while for the second, which is more general, this is no longer true, and the nonlinear normal modes theory is used. Based on this distinction, the control method is then separately applied to the two cases, and the optimal spatial and temporal shapes of the excitation are determined.

A detailed comparison of the infinite vs finite dimensional models is performed with respect to the control features, and it is shown that, depending on the boundary conditions, the control based on the reduced order model provides either exact or engineering acceptable results, although more systematic investigations are required to generalize the last conclusion.

Dipartimento di Ingegneria Strutturale e Geotecnica, Università di Roma "La Sapienza" via A. Gramsci 53, 00197, Roma, Italy, Giuseppe.Rega@uniroma1.it

Introduction

Nonlinear dynamics of many mechanical structures are governed by partial differential equations (PDEs), i.e., they are infinite dimensional. These equations are very difficult to be dealt with, or they are over-complicated, as the actual dynamics activates only few spatial modes, so that reduced order, finite dimensional, models governed by ordinary differential equations (ODEs) are usually introduced. Then, the question arises of the reliability of the involved approximations in providing accurate, qualitative and/or quantitative, descriptions of the true dynamics [START_REF] Rega | Multiple resonances in suspended cables: direct versus reduced-order models[END_REF] A c c e p t e d m a n u s c r i p t 2 [START_REF] Blekhman | Vibrational Mechanics[END_REF][START_REF] Nayfeh | Nonlinear Interactions[END_REF][START_REF] Steindl | Methods for dimension reduction and their application in nonlinear dynamics[END_REF]. Of course, the answer strongly depends on the kind of investigated phenomena, since, depending on various circumstances, reduced order models can be able to exactly, accurately or poorly describe various nonlinear dynamical features.

In this work this question is addressed with reference to the optimal control of the homoclinic bifurcation (herein identified as the first zero of the Melnikov function, see Remark 1 in the following) occurring in the finite planar dynamics of buckled beams, where the effects of the truncation due to the reduced dynamics are addressed against the actual full dynamics. In this respect, it is worth to point out that the theory and the development of this work are not related to the old and well-established theory of optimal control [START_REF] Kirk | Optimal Control Theory: An Introduction[END_REF], which has different objectives and uses different techniques. However, it is aimed at modifying, i.e. controlling, the system dynamics in an optimal way, and this is why we call it "optimal control."

The homoclinic bifurcation of the hilltop saddle (unstable rest position) is considered as the undesirable event to be eliminated, because it is well known that it triggers appearance of cross-well (scattered) dynamics, and that it is at the base of the principal manifestations of chaotic dynamics [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems and Bifurcation of Vector Fields[END_REF][START_REF] Holmes | Strange Attractors and Chaos in Nonlinear Mechanics[END_REF][START_REF] Thompson | Nonlinear Dynamics and Chaos[END_REF], although in general being not directly responsible for onset of chaotic attractors.

First, the homoclinic bifurcation is detected in a full infinite dimensional framework. Use is made of the Holmes & Marsden's theorem [1981] constituting an extension of the classical Melnikov's theory for analytical measuring of the distance between stable and unstable manifolds.

This permits detecting any kind of perturbations of the unperturbed homoclinic orbits, and therefore constitutes a straightforward generalization of the classical finite dimensional results [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems and Bifurcation of Vector Fields[END_REF][START_REF] Wiggins | Global Bifurcation and Chaos: Analytical Methods[END_REF].

Its elimination is then pursued by a control method previously developed by the authors [START_REF] Lenci | Controlling Nonlinear Dynamics in a Two-Well Impact System. Parts I & II[END_REF], 2004a], which consists of optimally modifying the shape of the excitation in order to shift as far as possible the critical threshold in the excitation amplitude-frequency parameter space.

It belongs to the family of control of chaos techniques [START_REF] Chen | From Chaos to Order: Methodologies Perspectives and Applications[END_REF][START_REF] Ott | Controlling Chaos[END_REF][START_REF] Fradkov | Chaos Control Bibliography[END_REF] sharing the common feature of exploiting dynamical chaotic properties to obtain specific results -usually, but not always, elimination of chaos -, irrespective of the nature of the required regularization and of the actual tools used to this aim.

The method has been formerly applied to various low dimensional mechanical systems and models [Lenci & Rega, 2003a, 2003b, 2004b[START_REF] Lenci | Control of pull-in dynamics in a nonlinear thermoelastic electrically actuated microbeam[END_REF][START_REF] Cao | Primary resonant optimal control for homoclinic bifurcations in single-degree-of-freedom nonlinear oscillators[END_REF], and its practical effectiveness has been confirmed by numerical simulations [Lenci & Rega, 2003a, 2003b]. The application of the method to an infinite dimensional system, which is one of the goals of this work and is herein pursued for the first time, appears therefore as the natural development of previous researches,

A c c e p t e d m a n u s c r i p t

3 which are herein embedded in a more general theoretical framework. Apart from its intrinsic interest, this extension also permits to confirm the great generality and versatility of the control method, already highlighted in [Lenci & Rega 2004a, 2004b].

The extension to infinite dimensional systems entails a major freedom in the choice of the optimal excitation. In fact, while in low dimensional systems only the temporal shape of the periodic excitation is modified, in the present case also the spatial shape can be changed, at least from a theoretical point of view. Thus, different optimization problems can be expected in principle, and this constitutes the main technical difference with respect to previous authors' works. However, it will be shown that the spatial shape can be preliminarily, and independently, determined in such a way that the temporal shape optimization proceeds just as in finite dimensional cases, so that the relevant results can be applied at once.

Another objective of the paper, as said, is to compare the results of the infinite dimensional analysis with those obtained with finite dimensional models. In this respect, it is known [START_REF] Troger | Nonlinear Stability and Bifurcation Theory[END_REF][START_REF] Steindl | Methods for dimension reduction and their application in nonlinear dynamics[END_REF]Rega & Troger, 2005] that obtaining reduced order models by the classical (linear) Galerkin method, which projects the dynamics on a planar subspace, may give incorrect results, even from a qualitative viewpoint (see, e.g., [START_REF] Arafat | Non-linear responses of suspended cables to primary resonance excitations[END_REF]), so that more refined analytical techniques are required to overcome this, and other, drawbacks.

A phenomenon of this type has been found in the present analysis. More precisely, it will be shown how for a class of system boundary conditions the classical Galerkin method correctly captures the homoclinic bifurcation, while for others it does not. Accordingly, the application of the control method to the first class can be correctly performed by means of the reduced order model, while in the other case different techniques, such as nonlinear normal modes, are preliminarily required in order to correctly detect the homoclinic bifurcation, and the reduced order model provides only approximate results.

In spite of this, however, we have found that, at least in a specific, yet technically meaningful, case the predictions made by the reduced order model are acceptable from the engineering viewpoint also for the second family of boundary conditions, although we have not proved this to be a robust property, so that it is not expected to hold for other parameter values or for other boundary conditions.

The paper is organized as follows. In Sect. 2 the mechanical model and the associated partial analyses of homoclinic orbit and homoclinic bifurcation, and on their consequences in terms of optimal control of the latter. Some mechanical and mathematical extensions of the analysis are then addressed in Sect. 7, and the paper ends with some conclusions (Sect. 8).

The mechanical model

The dimensionless partial differential equation governing the planar nonlinear dynamics of an initially straight buckled beam is

w & & +w′′′′+Γw′′-kw′′ ∫ 1 0 2 ) ' ( ζ d w =ε[F(z,t)-δ w & ], (1)
where w(z,t) is the time dependent transversal displacement of the point at abscissa z∈[0,1], dot and prime represent derivatives with respect to time and space, respectively, Γ>Γ cr is a dimensionless parameter and k is the dimensionless stiffness due to the membrane effects. The parameter Γ represents the first-order axial load (positive=compression, Γ cr =buckling critical parameter) associated with a finite prescribed end-displacement of the beam in the axial direction [START_REF] Villaggio | Mathematical Models for Elastic Structures[END_REF], δ is the coefficient of viscous damping, and

F(z,t)= ∑ = + N n n n t n z f 1 ) sin( ) ( ψ ω (2)
is the external T=2π/ω time periodic, spatially distributed, excitation. ε is a smallness parameter, introduced to stress that damping and excitations are small quantities; ε=0 is referred to as "unperturbed" or "conservative" problem.

Equation ( 1) is the simplest model which describes the finite amplitude dynamics of the beam. It is used because it permits feasible computations and analytical developments without losing the main mechanical phenomena, at least for moderate displacements, which are implicitly assumed without further remarks.

To make the problem (1) well-posed, initial and boundary conditions must be added. While initial conditions are not important, as we consider global system dynamics and not specific trajectories, boundary conditions (b.c. in the following) play a significant role in successive developments. Basically, we can have two different families of b.c., depending on whether they allow for a solution of the unperturbed problem in the form of separate variables or not.

Looking for an exact solution of (1) (with ε=0) in the form w(z,t)=a(t)b(z) requires b′′′′∝b′′∝b, namely, b′′=κb. If κ is positive, this equation has solution b(z)=c 1 sinh(z√κ)+ c 2 cosh(z√κ), while for κ negative the solution is b(z)=c 1 sin(z√(-κ))+c 2 cos(z√(-κ)). As we have only the two unknowns c 1 A c c e p t e d m a n u s c r i p t 5 and c 2 , in general there is no hope to satisfy all four b.c. associated to (1). This occurs only in some special cases, which constitute the first family. They are the hinged-hinged beam (Fig. 1a), which is the case considered by [START_REF] Holmes | A Partial Differential Equation with Infinitely Many Periodic Orbits: Chaotic Oscillations of a Forced Beam[END_REF] and [START_REF] Berti | Chaotic dynamics for perturbations of infinite dimensional Hamiltonian systems[END_REF], for which b(z)=√2 sin(πz), and the guided-hinged beam (Fig. 1b), for which b(z)=√2 cos(πz/2). The b.c. not fulfilling the previous requirement belong to the second family. Examples are the fixed-fixed (Fig. 2a), the fixed-hinged (Fig. 2b) and the fixed-free (Fig. 2c) beams.

The previous distinction has not only a mathematical character, but also an important dynamical meaning. In fact, for the first class of b.c. the invariant manifold Σ on which the homoclinic loop (which is the principal part of the dynamics herein considered) lies is flat (Sects. 3 and 4), while in the other case it is non planar (Sect. 6). This has strong consequences both on the amount of analytical computations and on the spatial modal shape associated with optimal control (Sect. 6).

In this paper, one representative case for each family is considered. The main part (Sects.

The homoclinic bifurcation by the Holmes & Marsden theorem

For buckled beams, the undeformed rest configuration is an unstable fixed point (saddle) in the appropriate, infinite dimensional, phase space. For conservative systems, the global unstable direction (which is unique for Γ in between the first and the second critical loads) of the hilltop saddle, which is a curve in the phase space and the swaying buckled spatial shape in the physical space, usually coincides with the corresponding global stable direction, and this special solution is the homoclinic orbit of the system. The remaining part of the phase space is a center manifold. A schematic picture of this situation (and of the invariant manifold Σ) is described in Fig. 3. When damping and excitations are added, the stable and unstable manifolds split, and, depending on the relative magnitude of damping and excitation, they keep disjoint or intersect. The intermediate critical case, which may be very complex and may involve various phenomena, corresponds to manifolds tangency, and represents the homoclinic bifurcation threshold, whose detection is central for the application of the control technique (Sect. 4), which is aimed at shifting this event as far as possible in parameter space.

The homoclinic bifurcation can be analytically detected by means of the [START_REF] Holmes | A Partial Differential Equation with Infinitely Many Periodic Orbits: Chaotic Oscillations of a Forced Beam[END_REF] abstract theorem (H.M. theorem in the following) -see the Remark 1 in the following -, which is a generalization to infinite dimensional systems of the classical Melnikov's method holding for finite dimensional systems [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems and Bifurcation of Vector Fields[END_REF]].

The theorem is quite technical, but the application to the case considered in this section (Fig. 1a)

A c c e p t e d m a n u s c r i p t 7 is relatively straightforward, and basically has been done in [START_REF] Holmes | A Partial Differential Equation with Infinitely Many Periodic Orbits: Chaotic Oscillations of a Forced Beam[END_REF]. Therefore, we refer to that work for the details and we simply summarize the hypotheses and the results, paying attention only to the aspects significant for the present work. The unique, important, modification is that the non-harmonic, but still periodic in time, excitation (2) is considered, while [START_REF] Holmes | A Partial Differential Equation with Infinitely Many Periodic Orbits: Chaotic Oscillations of a Forced Beam[END_REF] considered only harmonic excitation, i.e., N=1 in (2).

The first step consists in rewriting (1) as the evolution equation

dt dx =f 0 (x)+εf 1 (x,t)=[Ax+B(x)]+ε[A 1 x+f(t)], x∈X, (3)
in a Banach space X (see [START_REF] Brezis | Analyse fonctionnelle. Théorie et applications[END_REF] for the definitions of the Banach spaces used in the following). This is done by assuming x=(w,u) T , X= 2 0

H [0,1]×L 2 [0,1], i.e., x∈X ⇔ w∈ 2 0 H [0,1] and u∈L 2 [0,1], with the natural norm 2 X x =||x|| 2 =|w′′| 2 +|u| 2 = ∫ 1 0 2 ) ' ' ( ζ d w + ∫ 1 0 2 ζ d u
, and

Ax=         Γ - - ' ' ' ' ' ' w w u , B(x)=         ∫ 1 0 2 ) ' ( ' ' 0 ζ d w kw , A 1 x=         -u δ 0 , f(t)=         ) , ( 0 t z F , (4)
with f n (z)∈L 2 [0,1], which implies F(z,t)∈L 2 [0,1] for all t. Note that A:X⊃D(A)→X is an unbounded linear operator (i.e., there not exists a constant c such that ||Ax||<c||x||, ∀x∈D(A) ) with domain

D(A)={(w,u)∈H 4 ×H 2 |w=0
, w′′=0 and u=0 at z=0,1} dense in X (i.e., for every x∈X there exists an y∈D(A) which is arbitrarily close to x), B:X→X is a C ∞ nonlinear operator, A 1 :X→X is a bounded linear operator and f:IR→X is C ∞ and T-periodic.

There are four hypotheses to be satisfied.

(H1)

. A generates a C 0 one parameter group of transformations on X and the flow of ( 3) is defined for all t∈IR and for ε>0 sufficiently small.

This has been proved in [START_REF] Holmes | A Partial Differential Equation with Infinitely Many Periodic Orbits: Chaotic Oscillations of a Forced Beam[END_REF]] by theorems of [START_REF] Holmes | Bifurcation to divergence and flutter in flow induced oscillations: an infinite dimensional A c c e p t e d m a n u s c r i p t 30 analysis[END_REF] and implies that the dynamics is globally well defined in time.

(H2a) Assume that the unperturbed system is Hamiltonian with energy H:X→IR [Ω:X×X→IR is the weakly non-degenerate (i.e., Ω(u,v)=0 for all v implies u=0) associated symplectic form, such

that

Ω(f 0 (x),u)=dH(x)⋅u].
The symplectic form and the Hamiltonian are: (H2b) Assume there exists a symplectic 2-manifold Σ⊂X invariant under the flow of the unperturbed system and that on Σ the fixed point x=0 has an homoclinic orbit x h (t) (i.e., lim t→±∞ x h (t)=0, namely, lim t→±∞ w h (t)=0 and lim t→±∞ ) (t w h & =0, which actually are the "initial" conditions for the temporal boundary value problem).

Ω(x 1 ,x 2 )=Ω[(w 1 ,u 1 ),(w 2 ,u 2 )]= ∫ - 1 0 1 2 2 1 ) ( ζ d u w u w , (5a)
H(x)=H(w,u)= 2 1 |u| 2 - 2 Γ |w′| 2 + 2 1 |w′′| 2 + 4 k |w′| 4 . ( 5b
Hypothesis H2b is sketched in Fig. 3. For the b.c. of the first family the symplectic manifold is planar because the equation of motion admits solutions in the separate variable form (Sect. 2), while for those of the second family this is no longer true (see Fig. 3). In the considered case of Fig.

1a it is the plane (w,u)=(a 1 sin(πz), a 2 sin(πz)), (a 1 ,a 2 )∈IR 2 . The homoclinic orbit is

w h (z,t)=± π 2 k 2 π - Γ ) cosh( ) sin( 2 π π π - Γ t z , u h (z,t)= ) , ( t z w h & . ( 6
)
This of course requires Γ>Γ cr =π 2 . Indeed, we will assume in the following π 2 <Γ<4π 2 (see (H4)),

i.e., that the axial load belongs to the interval [Γ cr,1 , Γ cr,2 ] delimited by the first and second branching points. Note that, by symmetry, there are two homoclinic loops, one positive (called "right") and one negative (called "left").

(H3) Suppose that the linearized system dx

L /dt=[A+εA 1 ]x L +f(t) has a T-periodic solution x L (t,ε) such that x L (t,ε)=O(ε).
In this case x L (t,ε) can be computed exactly:

w L (z,t)= ∑ ∑ ∞ = =             Φ + 1 1 ) sin( sin( j N n n nj z j t n d π ω , u L (z,t)= ) , ( t z w L & , ( 7a
)
d nj = 2 2 2 2 2 2 2 2 ) ( ] ) ( [ ωεδ ω π π εγ n n j j nj + - Γ - , γ nj =2 ∫ 1 0 ) sin( ) ( ζ πζ ζ d j f n , Φ n =Φ n (ψ n ,ω,δ). (7b)
Eq. (7b) shows that (H3) is satisfied if and only if all the non-resonance conditions

n 2 ω 2 ≠j 2 π 2 (j 2 π 2 -Γ), for all n=1,2..N, j=2,3,…,∞ such that γ nj ≠0, (8)
hold. We remark that: (i) There is no non-resonance condition only for j=1, because of the next assumption (H4) π 2 <Γ<4π 2 , namely, the linear stiffness is negative on the first buckling mode, positive on the others. This explains why there is no non-resonance conditions for 1 d.o.f. models based on the first buckling mode [START_REF] Holmes | A Partial Differential Equation with Infinitely Many Periodic Orbits: Chaotic Oscillations of a Forced Beam[END_REF]]. (ii) There is a non-resonance condition to be satisfied only if γ nj ≠0, i.e., only if the j th buckling mode sin(jπz) is excited by the n th term f n (z) of the excitation.

C 2 ε≤dist( ε σ R ,|z|=1)≤C 1 ε, C 1 ,C 2 >0.
(H4c) The spectrum of A does not accumulate on 0.

A direct computation shows that the eigenvalues of e T(A+εA 1 ) are ε λ j T e , where

ε λ j = ] ) ( ) ( [ 2 1 2 2 2 2 2 π π εδ εδ j j - Γ + ± - . ( 9
)
The eigenvalues of e TA are 0 j T e λ , while the eigenvalues of A are 0 j λ . Thus, for the hinged-hinged beam, (H4) is satisfied if and only if π 2 <Γ<4π 2 . Generalizing to other b.c., this hypothesis is satisfied if and only if Γ cr,1 <Γ<Γ cr,2 , where Γ cr,j is the value of Γ for which 0 j λ vanishes. Note that the eigenvalues of A are the complex numbers λ satisfying λ 2 w+w′′′′+Γw′′=0 + b.c. for w(z)≠0 [START_REF] Holmes | A Partial Differential Equation with Infinitely Many Periodic Orbits: Chaotic Oscillations of a Forced Beam[END_REF]. Note also that a slightly different, but equivalent, definition of eigenvalues is used in [START_REF] Lenci | Dimension reduction of homoclinic orbits of buckled beams via nonlinear normal modes technique[END_REF].

This hypothesis implies that for the unperturbed system the saddle (rest position) has a one dimensional stable and a one dimensional unstable manifolds, which indeed coincide on the homoclinic loop, and an infinite dimensional center manifold (see Fig. 3). When the perturbations are added, the stable and unstable manifolds split, as said, and remain a one dimensional strongly stable and a one dimensional unstable manifolds. The center manifold, on the other hand, perturbs to a weakly stable manifold, so that the whole perturbed stable manifold becomes infinite dimensional.

The physical interpretation is that the system is dissipative outside u W ε (the perturbation of the unstable invariant manifold W u ), at least in the neighborhood of the saddle, and therefore it is expected that, if the excitation amplitude is not too high, the steady state dynamics contracts on s W ε (the union of the perturbation ss W ε of the stable invariant manifold W s and of the perturbation of the center manifold) due to damping. Thus, the steady state should be "similar" to that of the 1 d.o.f. model based on the first buckling mode, at least if the excitation amplitude is low. It will be proved in the following that this conjecture is mathematically exact for the b.c. of the first family, while it is engineering acceptable for the other family.

Once all the hypotheses have been verified, we are in a position of stating and applying the H.M. Before calculating the Melnikov function, we observe that, roughly speaking, the theorem says that simple zeros of M(m) imply the so-called "Melnikov chaos," i.e., fractal basins boundaries, transient chaos, sensitivity to initial conditions, onset of the chaotic saddle possibly representing the skeleton for chaotic attractors, etc., which are direct consequences of the stable and unstable manifolds intersection [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems and Bifurcation of Vector Fields[END_REF][START_REF] Holmes | Strange Attractors and Chaos in Nonlinear Mechanics[END_REF]. The existence of the chaotic attractor is however not guaranteed.

Remark 1. In this paper we identify the homoclinic bifurcation, i.e. the first tangency between the perturbed stable and unstable manifolds of the saddle ensuing from the rest position, with the threshold at which the Melnikov function M(m) has its first (double, and thus non-simple) zero.

Technically speaking, in some cases we can have Melnikov chaos even below this threshold, because when x=0 has a homoclinic tangle on one side (say, the intersection of W s,l and W u,l is nontrivial), the manifold W u,l winds itself around W u,r infinitely often, intersecting W s,r before W u,r does.

Thus, the "first" homoclinic tangency on the other side, which would actually correspond to the relevant homoclinic bifurcation, does no longer exist. However, it is difficult to estimate the size of this effect, which is typically small, and we consider acceptable our identification from a practical point of view.

With respect to the purposes of control, on the other hand, this fact is unessential. In fact, above the Melnikov threshold there certainly exists Melnikov chaos, and thus increasing this threshold is worthy in any case, even if some (spurious) chaos may occur for values of the excitation slightly lower than this critical threshold.

After some computations we obtain that the Melnikov functions for right and left homoclinic loops are given by

M r,l (m)=-δ ∫ ∫ ∞ ∞ - dt dz t z w h 1 0 2 ) , ( & + ∫ ∫ ∞ ∞ - + dt dz m t z F t z w h 1 0 ) / , ( ) , (
k 2 / 3 2 ) ( - Γ δ ∑ =         - Γ + N n n n n nm n k 1 2 1 2 cosh ) cos( π ω ψ ω π γ m . ( 10
)
It is worth to note that only the terms

γ n1 =2 ∫ 1 0 ) sin( ) ( ζ πζ ζ d f n
are involved in (10). This is due to the flatness of the invariant manifold Σ and no longer holds for the second family of b.c. (see eq. ( 41)).

In the following it is convenient to rewrite (10) in the form

M r,l (m)=const.[1± ) ( , 11 11 ω γ γ h cr h(m)], (11)
where

) ( , 11 ω γ h cr =δ k ω π 3 ) ( 4 2 / 3 2 - Γ         - Γ 2 2 cosh π ω , h(m)= ∑ = + N n n n nm h 1 ) cos( ψ , h n = 11 1 γ γ n         - Γ         - Γ 2 2 2 cosh 2 cosh π ω π ω n n . ( 12
)
We prefer to use the representation (11) because we assign γ 11 the role of overall excitation amplitude, while the relative amplitudes of the superharmonics γ n1 /γ 11 and the phases ψ n determine the (temporal) shape of the excitation [Lenci & Rega, 2003a, 2004a]. The effects of the superharmonics on the Melnikov functions are instead governed by the parameters h n , n=2,3,…,N (note that h 1 =1, and that the amplitude-free, oscillating part h(m) of M r,l (m) is 2π-periodic with zero mean value).

According to the H.M. theorem, we have homoclinic intersection of the right (left) manifolds if there exists m∈[0,2π] such that M r,l (m) have simple zeros, respectively. As in [Lenci & Rega 2003b], we see that the equation M l (m)=0 has a simple solution if and only if is called saved (i.e., controlled) region (Fig. 4), and represents the zone where the unharmonic excitation is theoretically effective. Its maximum enlargement constitutes the objective of the control method.

γ 11 > l h cr M ) ( , 11 ω γ = ) ( , 11 ω γ l cr , M l =max m∈[0,2π] {h(m)}, ( 13
To quantitatively measure the increment of the critical thresholds of the unharmonic excitations with respect to the harmonic reference case we introduce the gains [Lenci & Rega, 2004a, 2004b],

which are the ratios

G l = ) ( ) ( , 11 , 11 ω γ ω γ h cr l cr = l M 1 , G r = ) ( ) ( , 11 , 11 ω γ ω γ h cr r cr = r M 1 , (15)
A c c e p t e d m a n u s c r i p t

13 and which will be useful in the following.

Optimal control

In this section we apply to the buckled beam a method for controlling nonlinear dynamics and chaos developed by the authors [Lenci & Rega, 2004a, 2004b] and previously applied only to single d.o.f. mechanical systems [START_REF] Lenci | Controlling Nonlinear Dynamics in a Two-Well Impact System. Parts I & II[END_REF], 2003a, 2003b]. Those analyses are herein extended to the infinite dimensional system (1).

The idea of the method is to increase the ) ( , , 11 ω γ l r cr thresholds or, equivalently, to enlarge the saved region or increase the gain. This result is pursued by optimally varying both the spatial (Sect.

4.1) and temporal (Sect. 4.2) shapes of the excitation, and it is what we call "optimal control of homoclinic bifurcation." In fact, the two problems can be analyzed separately.

"Optimizing" the spatial shape of the excitation

Still referring to the first family of b.c., the optimal choice of spatial shapes f n (z) of the excitation is based on the observation that, because of the flatness of Σ, only the coefficients γ n1 are involved in the expression of the Melnikov function ( 10). Accordingly, only the spatial shape of the excitation parallel to the invariant manifold Σ plays a role in the Melnikov analysis, and then on control. This agrees with the physical considerations suggesting that, due to damping, the significant part of the dynamics occurs near Σ.

As a consequence, if we assume f n (z) orthogonal to Σ, i.e., f n (z)= ∑ ∞ =2

) sin( j nj z jπ γ , we have no homoclinic bifurcation at all, and the problem is meaningless. Thus, f n (z) must contain at least the term sin(πz). Indeed, we just assume f n (z)=γ n1 sin(πz), which implies that the useless coefficients

γ nj =2 ∫ 1 0 ) sin( ) ( ζ πζ ζ d j f n
, j=2,3,…∞, vanish identically. This entails two meaningful circumstances:

(i) being sin(πz) the only spatial component of f n (z) able to modify the homoclinic bifurcation, this choice permits to minimize the cost of control, because we only use what is effective;

(ii) having γ nj =0, j=2,3,…∞, the non-resonance conditions (8) become meaningless, so that the control, which would not hold in the present form for the frequencies not satisfying eqs. (8), is applicable everywhere. This is practically important, because if N is large enough (i.e., if we use a lot of superharmonics, which provide the best optimal gains (see Tabs. A1 and A2 in the appendix)), the resonant frequencies are almost dense in IR and the control would not be applicable at all.

The previous points show why we consider practically "optimal" the assumption f n (z)=γ n1 sin(πz).

A c c e p t e d m a n u s c r i p t 14

In fact, contrary to what happens with the temporal shape (see Sect. 4.2), here the choice does not ensue from a mathematical optimization problem, but just from mainly qualitative considerations.

From a practical point of view, the assumption f n (z)=γ n1 sin(πz) means that we do not excite the vibrating modes; from a theoretical point of view, on the other hand, it means that the perturbed manifold Σ ε remains planar and coincides with Σ.

Optimizing the temporal shape of the excitation

The issues of choosing the optimal temporal shape, i.e. the optimal superharmonics, is now addressed. It is quite clear from ( 13)-( 14) that if we want to enlarge as much as possible the saved region we have to increase the gains (15). However, as shown in [Lenci & Rega, 2003b, 2004a],

this requires some care, because the presence of two homoclinic orbits permits to choose among different control strategies. Indeed, we can control only the right (left) homoclinic bifurcation,

irrespective of what happens to the other, or we can control simultaneously the right and the left homoclinic bifurcations.

The two approaches are complementary rather than competing [Lenci & Rega, 2003b, 2004a]:

the first is aimed at obtaining a topologically "localized" control and provides "large" optimal gains (see Tab. A1), whereas the second is aimed at controlling, on average, the "whole" phase space, but provides smaller optimal gains (see Tab. A2). Maximize G r (G l ) by varying the coefficients h j and ψ j , j=2,3,… of h(m).

(16)

These problems have been solved in [Lenci & Rega, 2003b, 2004a], and the solution is ψ j =0, while the optimal gains and the optimal h j for increasing number N of added superharmonics are given in Tab. A1, which refers to the "right" case. The solution of the "left" case is that of Tab. A1

with even coefficients h 2j changed of sign.

Table A1 shows that, as expected, the optimal gain is an increasing function of the number N of controlling superharmonics. It is bounded from above by 2, which corresponds to the (mostly theoretical) case of infinite superharmonics and shows that it is possible to double the critical excitation threshold. This limit case is reported in Fig. 4. From a practical point of view, very satisfactory optimal gains can be obtained with even few controlling terms.

Using the coefficients of Tab. A1 in eq. ( 12) 2 gives the physical ω-dependent optimal excitations, which for N=2 and N=3 and in the "right" case are reported below:

            - Γ - Γ ) 2 /( cosh 2 ) 2 /( 2 cosh 2 2 π ω π ω sin(2ωt)], (17) F(z,t)=γ 11 sin(πz)[sin(ωt)+0.5528             - Γ - Γ ) 2 /( cosh 2 ) 2 /( 2 cosh 2 2 π ω π ω sin(2ωt)+0.1708             - Γ - Γ ) 2 /( cosh 3 ) 2 /( 3 cosh 2 2 π ω π ω sin(3ωt)].
The second of the previous excitations, divided by γ 11 sin(πz), is reported in Fig. 5a . a) optimal one-side "right" excitation, i.e., (17) 2 . Note that in this case there are two controlling superharmonics added to the reference basic harmonic. b) optimal global excitation, i.e., (19) 1 . In this case there is only one added superharmonic.

ii) "Global" control. To control simultaneously the right and the left homoclinic bifurcations, the right and the left gains, G r and G l , must be increased at the same time. Mathematically, this entails to increase their minimum value, namely, to solve the following optimization problem:

Maximize G=min{G r ,G l } by varying the coefficients h j and ψ j , j=2,3,… of h(m).

(18)

The solution of ( 18) is given by ψ j =0, h 2j =0 (which automatically guarantees that max{h(m)}= -min{h(m)}, i.e., that G r =G l ) and by the coefficients reported in Tab. A2 [Lenci & Rega, 2003b, 2004a].

As in the "one-side" control, the optimal gains are still increasing functions of N, but now they are strongly lesser than those of Tab. A1, showing how global control is theoretically less performant, even in the limit case N=∞ (see Fig. 4). The counterpart of this drawback is the possibility to control simultaneously both homoclinic bifurcations.

A c c e p t e d m a n u s c r i p t

16 By the coefficients of Tab. A2 and by eq. ( 12) 2 we can compute the first two ω-dependent optimal physical excitations (N=3 and N=5, respectively):

F(z,t)=γ 11 sin(πz)[sin(ωt)-0.1667             - Γ - Γ ) 2 /( cosh 3 ) 2 /( 3 cosh 2 2 π ω π ω sin(3ωt)], (19) F(z,t)=γ 11 sin(πz)[sin(ωt)-0.2323             - Γ - Γ ) 2 /( cosh 3 ) 2 /( 3 cosh 2 2 π ω π ω sin(3ωt)+0.0610             - Γ - Γ ) 2 /( cosh 5 ) 2 /( 5 cosh 2 2 π ω π ω sin(5ωt)].
The first of the previous excitations, divided by γ 11 sin(πz), is reported in Fig. 5b for various values of the parameter p=

) 2 /( 2 π ω - Γ
. By comparing the two pictures of Fig. 5 it is possible to note how the one-side optimal excitation has slightly more pronounced peaks, and that for increasing values of p the third order superharmonic tends to become predominant.

Comparison with the single d.o.f. reduced model

In [Lenci & Rega, 2003b] the control method has been applied to the single d.o.f. approximation of the dynamics of a buckled beam. In this section we compare the present much more general analysis with that reported therein.

The single d.o.f. approximation is obtained by the Galerkin method by projecting the dynamics onto the planar manifold spanned by the local stable and unstable directions, which for the first class of b.c. coincide with Σ. Therefore, we assume w (z,t)=a(t)sin(πz), (20) insert this in (1), multiply by sin(πz) and integrate from 0 to 1. This yields the following ordinary differential equation:

a& & -π 2 (Γ-π 2 )a+k 2 4 π a 3 =ε[ ∑ = + N n n n t n 1 1 ) sin( ψ ω γ -δ a & ], γ n1 =2 ∫ 1 0 ) sin( ) ( ζ πζ ζ d f n . ( 21
)
Equation ( 21) is the classical Duffing equation with negative linear stiffness (being Γ>π 2 ) which, in a different notation, has been investigated in [Lenci & Rega, 2003b], where it has been shown that the right and left homoclinic orbits are given by

a h (t)=± π 2 k 2 π - Γ ) cosh( 1 2 π π - Γ t . ( 22
)
Equation ( 22), together with (20), gives exactly (6). The classical, two dimensional, Melnikov

A c c e p t e d m a n u s c r i p t

17 function [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems and Bifurcation of Vector Fields[END_REF]] can then be computed, which is shown [Lenci & Rega, 2003b] to be exactly given by ( 10). This proves that, in the present case, even the single d.o.f. reduced order model is able to exactly capture the homoclinic bifurcation, and thus the theoretical developments and the practical implementations of the control method are identical for the finite and infinite dimensional models.

This is a consequence of the flatness of Σ, and no longer holds for the second family of b.c.

(Sect. 6). In any case, it must be emphasized that "the full power of the (Holmes & Marsden) theorem is necessary since in the infinite dimensional case, the perturbed manifolds … do not lie in Σ" [START_REF] Holmes | A Partial Differential Equation with Infinitely Many Periodic Orbits: Chaotic Oscillations of a Forced Beam[END_REF], and therefore they are only approximated by the single d.o.f. model.

These perturbations are generic and can be due to external excitation and/or nonlinear coupling between different modes, and their negligibility, now rigorously proved by the H.M. theorem, couldn't be easily guessed a priori.

Since in the specific case of Fig. 1a, and more generally for the first family of b.c., the application of the control method is identical for finite and infinite models, we can take advantage of the results in [Lenci & Rega, 2003b], where some numerical simulations are reported aimed at checking the theoretical predictions and the practical performance of control. In particular, the oneside and global controls are numerically compared, and it has been shown:

1) the numerical reliability of the Melnikov prediction of homoclinic bifurcations, which can be usefully applied in spite of its perturbative nature;

2) the control induced regularization of fractal basin boundaries, which confirms the first predicted effect of elimination of the homoclinic bifurcation;

3) the reduction of fractal erosion of the basins of attraction of confined attractors, which is very useful for reducing the sensitivity to initial conditions; 4) the effect of control also beyond theoretical expectations, permitting, for example, the delay of the in-well to cross-well chaos transition in the case of one-side control.

Different boundary conditions -non planar manifold Σ

In this section the more general case of b.c. belonging to the second family is considered to complete the application of the control method to infinite dimensional dynamical systems. The analysis herein developed is general, but to fix ideas, and to make explicit computations, we refer to the fixed-fixed case of Fig. 2a. Furthermore, we choose Γ=40, which belongs to the interval [Γ cr,1 , 39.478, 80.763] (hyp. (H4) of the H.M. theorem).

Γ cr,2 ]=[4π 2 , 8.183π 2 ]=[

Nonlinear normal mode technique for homoclinic solution

The main difficulty is that the invariant manifold Σ of hyp. (H2b) is not planar and it is difficult

A c c e p t e d m a n u s c r i p t

18 to be handled. To overcome this point, the nonlinear normal modes technique is used [START_REF] Shaw | Normal modes for nonlinear vibratory systems[END_REF][START_REF] Vakakis | Normal Modes and Localization in Nonlinear Systems[END_REF][START_REF] Nayfeh | Nonlinear Interactions[END_REF]. Together with other related methods, such as the approximate inertial manifold and the proper orthogonal modes, it is one powerful technique developed to obtain reliable reduced order models, which is a matter of intensive on-going research (see also a recently published Special Issue [Rega & Troger, 2005]). Although these techniques were originally developed to deal with nonlinear oscillations, in this paper their application is concerned with homoclinic orbits [START_REF] Lenci | Dimension reduction of homoclinic orbits of buckled beams via nonlinear normal modes technique[END_REF].

To apply the method, the first N eigenfunctions of the linear operator A are required. Because A is self-adjoint and its domain is dense, we get the exact manifold Σ in the limit N→∞. For w( 0)=w′( 0)=w( 1)=w′(1)=0 and for Γ=40 the first five eigenvalues of A are:

λ 1 =2.620687, λ 2 =44.085549i, λ 3 =103.232438i, λ 4 =181.871769i, λ 5 =280.309519i, ( 23
)
while the expressions of the associated eigenfunctions are given in the appendix (eq. ( A1)), and are depicted in Fig. 6. They have been normalized with 1 )] ( [

1 0 2 = ∫ ζ ζ d w j
, which simplifies the inertial terms in the following expression (25) of the Hamiltonian. Continuous lines are odd eigenfunctions, and dashed lines are even eigenfunctions.

Note that, according to Γ∈[Γ cr,1 , Γ cr,2 ], only the first eigenvalue is real, while the others are purely imaginary. This property, in particular, guarantees that the infinite non-resonance conditions n 2 λ 1 2 ≠λ n 2 , n=2,3,4…, which are required for the validity of the procedure, as shown in [START_REF] Lenci | Dimension reduction of homoclinic orbits of buckled beams via nonlinear normal modes technique[END_REF], are automatically satisfied.

Due to symmetric boundary conditions, eigenfunctions with odd index are symmetric (i.e., w(0.5+δ)=w(0.5-δ), see the continuous lines in Fig. 6), while those with even index are antisymmetric (i.e., w(0.5+δ)=-w(0.5-δ), see the dashed curves in Fig. 6). Consequently, the odd

A c c e p t e d m a n u s c r i p t

19 eigenfunctions span the symmetric part of the domain of A, while the even eigenfunctions span its complement, i.e., the anti-symmetric part of the configuration space. As the homoclinic solutions are expected to be spatially symmetric, because they approach the unstable symmetric eigenfunction w 1 (z) for t→∞, we neglect the even eigenfunctions. This has no consequences on the expected solution since even and odd terms are not coupled in the reduced order Hamiltonian, so that a solution initially symmetric remains symmetric for all times.

Accordingly, we directly assume a 2j =0, j=1,2,…, and look for the solution in the form

w(z,t)=a 1 (t)w 1 (z)+a 3 (t)w 3 (z)+a 5 (t)w 5 (z)+…, u(z,t)= w & (z,t). ( 24
)
In the following we consider terms only up to w 5 (z) to limit the computations, but of course the analysis can be trivially extended. By inserting (24) in the Hamilton function (5b) we get

H =0.5( 2 1 a & + 2 3 a & + 2 5 a & )+(-3.43400 2 1 a +5328.47013 2 3 a +39296.0690 2 5 a )+ +k(6.58790 2 1 a +49.53020 2 3 a +132.04221 2 5 a -12.58343a 1 a 3 -7.67163a 1 a 5 -26.57427a 3 a 5 ) 2 = =0.5( 2 1 a & + 2 3 a & + 2 5 a & )+V(a 1 ,a 3 ,a 5 ). ( 25
)
The equation of motion can then be obtained by the classical Hamiltonian formalism from (25):

j a& & =- j a V ∂ ∂ , ( 26
)
and it is trivial to show that they admit the saddle equilibrium point a j =0 with H=0.

Generalizing the treatment of [START_REF] Rand | Lecture Notes on Nonlinear Vibrations[END_REF] (see also the introductory paper in [Rega & Troger, 2005] for an overall framework), we identify a main, or master, variable (i.e.. modal amplitude) x and consider the others as secondary, or slave, variables (higher order modal amplitudes) y i .

Because in the buckled beam problem the unstable direction triggering the homoclinic orbit is just along the first eigenfunction, it is natural to assume x=a 1 and y 1 =a 3 , y 2 =a 5 and so on. Equations ( 26)

then become x& & =- x V ∂ ∂ , i y& & =- i y V ∂ ∂ . ( 27
)
The key idea of the method consists in assuming that the slave variables can be expressed as time-independent functions of the master one:

y i =y i (x). (28)
If we use only N terms in (24) the relation ( 28) provides an approximation of the searched manifold Σ (approximate manifold), while for N→∞ it exactly gives Σ (exact manifold). We will see that, in general, few terms are indeed sufficient to practically give an adequate approximation of Differentiating (28) with respect to time we get, by the chain rule,

i y & = dt dy i = dx dy i dt dx =y i ′ x & , i y& & = y i ′′ x & 2 +y i ′ x& & . ( 29
)
Combining (29) 2 with (27) yields

i y V ∂ ∂ =-y i ′′ x & 2 +y i ′ x V ∂ ∂ . ( 30
)
To apply the previous general treatment to the homoclinic orbit, we note that on these solutions we have H=0, namely,

0=( 2 x & + ∑ - = 1 1 2 N i i y & )+2V(x, y i (x))= 2 x & [1+ ∑ - = 1 1 2 ) ' ( N i i y ]+2V(x, y i (x)). ( 31
)
From ( 30) and ( 31) we can eliminate x & and finally get the strongly nonlinear ordinary differential equations permitting to determine the unknown functions (28):

2 Vy i ′′+[1+ ∑ - = 1 1 2 ) ' ( N i i y ]         ∂ ∂ - ∂ ∂ i i y V x V y ' =0. (32)
Equations ( 32) are the slave amplitudes equations which particularize the modal equations of [START_REF] Rand | Lecture Notes on Nonlinear Vibrations[END_REF]. Indeed, they are much more complex than the original system (27), and its introduction may appear useless. Notwithstanding, it contains only polynomial terms, and it is natural to look for a relevant solution in the Taylor series form [START_REF] Lenci | Dimension reduction of homoclinic orbits of buckled beams via nonlinear normal modes technique[END_REF]]

y i = k 1 [y i,1 (x√k)+y i,3 (x√k) 3 +y i,5 (x√k) 5 +y i,7 (x√k) 7 +…], (33)
where even terms have not been considered because of the symmetry of the system, which requires that if w(z,t) is a solution of (1), so is -w(z,t). Moreover, the constant term vanishes also because the manifold contains the trivial saddle, i.e., y i (0)=0. The parameter k has been explicitly singled out in (33) to have k independent numbers y i,l .

To determine the unknowns y i,l , we insert (33) in (32) and expand in Taylor series. Equating to zero the coefficients of x we get

c i [1+ ∑ - = 1 1 2 1 , ) ( N j j y ] y i,1 =0 (34)
(c i are constants different from zero). Since complex values of y i,1 cannot be accepted, the unique solution of ( 34) is y i,1 =0. This shows that, in a neighborhood of x=0, a unique manifold Σ is detected Equating to zero the coefficients of the successive powers of x we get a recursive series of linear algebraic equations permitting to determine the remaining constants y i,l . The first few of them are reported in Tab. 1, and they are used in Fig. 7 to draw the solutions y 1 (x) and y 2 (x).

l y 1,l → y 1 (x) y 2,l → y 2 (x) 1 0 0 3 0.154679×10 -1 0.128512×10 -2 5 -0.810240×10 -3 -0.108526×10 -4 7 -0.472535×10 -4 -0.261876×10 -6 Table 1. The Taylor coefficients y i,l of the slave amplitudes y i (x), i=1,2, for Γ=40. . The slave amplitudes y 1 (x) and y 2 (x) (up to the 7 th order) for Γ=40, and a schematic section of a hypothetically bending exact manifold Σ.

Several conclusions can be drawn from Tab. 1 and Fig. 7. The most important is that, at least for the considered value of Γ, the slave variables are orders of magnitude lesser than the master one.

Furthermore, the degree of smallness of y i (x) increases with its order i, so that high order modal amplitudes (say, greater then 2) are practically unessential.

When considering from an engineering point of view the previous comments and the fact that Σ vanishes cubically at x=0, we can conclude that the non-flatness of the manifold, while theoretically important, has a minor practical relevance. This does not appear completely surprising since the considered axial load is not much larger than the first buckling load (Γ/Γ cr,1 =1.013). In [START_REF] Lenci | Dimension reduction of homoclinic orbits of buckled beams via nonlinear normal modes technique[END_REF]] it has been shown how moving away from it, the non-flatness of the manifold becomes more important.

To check the reliability of the approximation of y 1 (x) and y 2 (x) we have computed the Remark 2. It is important to note that the previous treatment fails if the assumption (28) no longer holds, e.g., if the manifold Σ containing the homoclinic orbit bends (as schematically shown by the dashed curve of Fig. 7) and the orbit overcomes the turning point. As a general consideration, we can say that this does not occur for small values of Γ-Γ cr,1 , as assumed in this section, but this phenomenon might not happen even for not to small values of Γ-Γ cr,1 .

Minimal order homoclinic solution and refined order homoclinic bifurcation

By inserting the expression of y 1 (x) and y 2 (x) obtained in the previous section in the expression

(25) of the potential V(x, y i (x)) we get

V (x)= k 1 [-3.434000(x√k) 2 +43.400440(x√k) 4 -1.354661(x√k) 6 +0.194647(x√k) 8 -0.213589×10 -1 (x√k) 10 +0.629506×10 -3 (x√k) 12 +…] (35)
and the equation of motion of x is (27) 1 . Equation ( 35) contains the effect of up to the second higher (slave) mode, each relevant contribution being evaluated by accounting for up to the seventh order nonlinearity in the Taylor expansion (33). It describes an hardening potential with two symmetric wells, surrounded by two symmetric homoclinic loops x h (t), which are of primary interest for this work. In the phase space, they extend from x h (∞)=0 to x h (0)= x , where x is the unique non null solution of V(x)=0. We get from (35) x =0.281634/√k.

The 4 th order approximation of the potential ( 35)

V (x)= k 1 [-3.434000(x√k) 2 +43.400440(x√k) 4 ]=-3.434000x 2 +43.400440kx 4 (36)
corresponds to a classical Duffing oscillator -which can also be obtained via a linear Galerkin projection -, and only accounts for the effect of the master mode. To check the validity of this approximation we have computed x by (36) and we got x =0.281289/√k. As this is practically coincident with the "true" x , and because we are interested only in the dynamics in the range [0,

x ], we will use (36) in the following. In terms of the reduced order models of homoclinic orbits discussed in [START_REF] Lenci | Dimension reduction of homoclinic orbits of buckled beams via nonlinear normal modes technique[END_REF], this is equivalent to consider the minimal (or unrefined) one, which was shown to work very well close to the first buckling load, as we are doing herein.

The equation of motion associated to (36) and the embedded homoclinic orbits are , see [START_REF] Lenci | Dynamics of shallow elastic arches. Part I: chaotic response of harmonically shaped arches[END_REF])

a& & + 2 4 π (h 2 -4-3h 4 2 - h )a- 2 3 4 π (h- 4 2 - h )a 2 +π 4 a 3 =ε[ ∑ = + N n n n t n 1 1 ) sin( ψ ω γ -δ a & ], (47)
instead of the Duffing equation ( 21).

The control method has been applied to (47) in [Lenci & Rega, 2004b], which reports all the theoretical developments and shows how in this case the global control is no longer system independent, contrary to what happens for one-side controls and for global control of the symmetric (Duffing) system, where the solutions are system-independent.

The situation of global control is relatively more involved, and two different approaches have been followed, based on increasing of either gains (global control without symmetrization) or excitation amplitudes (global control with symmetrization). In the latter case the symmetrization of system dynamics can be actually "achieved" or simply "pursued" [Lenci & Rega, 2004b].

As a consequence of the flatness of Σ for the considered undeformed shape and boundary conditions, the results of [Lenci & Rega, 2004b] can be immediately applied to the infinite dimensional case. Of course, the extension to the second family, i.e. to different initial shapes and to different b.c., requires an analysis similar to that of Sect. 6, and will not be pursued here.

Berti and Carminati theorem

The homoclinic bifurcations of infinite dimensional systems can be analytically detected also by the [START_REF] Berti | Chaotic dynamics for perturbations of infinite dimensional Hamiltonian systems[END_REF] theorem (see also [START_REF] Alessio | Heteroclinic motions joining almost periodic solutions for a class of Lagrangian Systems[END_REF]), which indeed is an extension of the H.M. theorem.

The Berti and Carminati theorem permits to deal with the following cases, which do not fulfill the H.M. theorem hypotheses: 1) the excitation does not need to be periodic, but can also be quasi-periodic and even more generic;

2) parametric excitations can easily be considered, while the application of the H.M. theorem to this case requires some technicalities;

3) the zeros of the Melnikov functions do not need to be simple. This is mainly a technical extension with minor practical consequences, because the non simple zeros of M(m) are structurally unstable; 4) it can be applied to the case with more than one unstable direction, which is excluded by the H.M. theorem. This means that hypothesis (H4) can be relaxed, and we can consider an axial load Γ larger than the second branching point Γ cr,2 ; 5) it can be applied to a vectorial version of (1), this permitting a more accurate modeling of the In addition to the previous advantages, however, there is also one disadvantage. In fact, the proof of Berti and Carminati theorem is strongly related to the buckled beam with hinged-hinged b.c., and the extension to other b.c., even of the first family, does not appear trivial.

Conclusions and further developments

A method for controlling nonlinear dynamics and chaos through elimination, or better, shift in parameter space, of the principal homoclinic bifurcation has been applied to an infinite dimensional model of a buckled beam. The method takes advantage of the partially free temporal and spatial shapes of the excitation, which is however constrained to be periodic in time, and allows to determine the relevant optimal shapes through either the solution of appropriate mathematical optimization problems, or practical considerations.

The theoretical analysis is based on an abstract theorem of Holmes and Marsden, which is an infinite dimensional version of the classical Melnikov theory, permitting the analytical detection of the homoclinic bifurcation of the main saddle. This theorem has been applied to the present case, and the influence of the excitation shape on the Melnikov function has been determined.

Two mechanically and mathematically different cases have been identified and studied separately. From a dynamical point of view, the two cases differ because in the first the homoclinic loops lie on a flat manifold, while in the second they lie on a non-flat manifold. From a mechanical point of view, they correspond to different beam boundary conditions.

The simpler case of flat manifold has been investigated first. After detecting the homoclinic bifurcation threshold by the H.M. theorem, the control method has been applied. The "optimal" spatial shape of the excitation has been identified as the one parallel to the flat manifold, which solely plays a role in the Melnikov analysis, whereas the cases of one-side and global control have been discussed and compared with each other in addressing its optimal temporal shape. It has been shown that the application of the method entails the same optimization problem found in other discrete mechanical systems previously investigated by the authors. While confirming the great generality of the method, this permits to take advantage of the results of other works to get the optimal temporal excitations.

The results of the infinite dimensional analysis have then been compared with those obtained with the reduced (indeed, one d.o.f.) model. It has been shown that the latter is able to exactly detect the homoclinic bifurcations, although the full power of the H.M. theorem is required to take into account out of manifold perturbations, so that the application of the control method is exactly the same for both finite and infinite dimensional models. This permits to refer to other works for The second case of non-flat manifold has then been considered. The nonlinear normal modes theory has been used to detect the non planar manifold and the embedded homoclinic loops. The homoclinic bifurcation threshold has then be computed again by the H.M. theorem. It has been shown by a numerical example that for axial loads only slightly larger than the buckling load, the non planarity of the manifold, while theoretically important, is practically negligible.

The control method has then been applied also to this case, showing how it is possible to choose the spatial shape of the excitation in such a way to have the same optimization problem of the flat manifold case (but the non planarity of the manifold does not disappear at all). Thus, the application of the control method proceeds exactly as in the previous case.

The comparison with the reduced order model has also been performed. It has been shown that, because of the practical negligibility of the manifold non-flatness, the reduced order model practically gives the same results as the full model. This theoretically differs from the previous case, where the reduced order model gives exactly the correct results.

Finally, some hints about the possible extension of the infinite dimensional approach to control of homoclinic bifurcation of systems with initial curvature have been given.

As far as further developments are concerned, we can quote, among others, at least the following, which appear to be particularly worthy to improve the knowledge of the whole matter:

1) detailed application of Berti & Carminati theorem to deal with, e.g., non-periodic excitations;

2) physical experiments;

3) more involved mechanical models of the beam; 4) parametric excitations. Figure 7. The slave amplitudes y 1 (x) and y 2 (x) (up to the 7 th order) for Γ=40, and a schematic section of a hypothetically bending exact manifold Σ.

  , Costruzioni e Strutture, Università Politecnica delle Marche via Brecce Bianche, 60131, Ancona, Italy, Lenci@univpm.it GIUSEPPE REGA

  differential equation are set. The homoclinic bifurcation of the hilltop saddle is then analytically detected by the Holmes & Marsden theorem (Sect. 3) for the case of hinged-hinged beam. These results are then employed in Sect. 4, where the control method is applied and the optimal spatial and temporal shapes of the excitation are determined. The comparison with the single d.o.f. reduced in Sect. 5. Section 6 deals with the more intriguing case of fixed-fixed boundary conditions via the nonlinear normal mode technique, by dwelling upon the infinite vs reduced order

Figure 1 .Figure 2 .

 12 Figure 1. Examples of b.c. belonging to the first family. a) hinged-hinged: w(0)=w′′(0)=0, w(1)=w′′(1)=0, b) guided-hinged: w′(0)=w′′′(0)+Γw′(0)=0, w(1)=w′′(1)=0.

  developed with reference to the hinged-hinged b.c., because it permits easier computations and easier development of the control technique. Successively (Sect. 6), to complete the analysis, the more general case of a fixed-fixed beam is considered as a representative of the second family of b.c..

  Figure3. A schematic picture of the dynamical system structure in the unperturbed case. The center manifold is perpendicular to Σ at the saddle.

  symplectic form will play a key role in the Melnikov function derived in the following, which is central for the purposes of this work.

  For ε=0, e TA has a spectrum consisting in two simple real eigenvalues e ±λT , λ≠0, with the rest of the spectrum on the unit circle.(H4b) For ε>0, e T(A+εA 1 ) has a spectrum consisting in two simple real eigenvalues inside the unit circle |z|=1 and verifying the estimates

  that M(m) has a simple zero as a function of m. Then for ε sufficiently small, the stable s W ε and the unstable u W ε manifolds of the saddle x L (0,ε) of the stroboscopic Poincarè map associated with the flow intersect transversally.

,

  )while the equation M r (m)=0 has a simple solution if and only if M r =-min m∈[0,2π] {h(m)}. depicted in Fig.4in the frequency/amplitude parameter space (ω, γ 11 ), are identified as the homoclinic bifurcation thresholds for the right (left) homoclinic loop. Above these critical curves right (left=M l /M r =const.), and coincide only for the excitations satisfying M l =max{h(m)}= -min{h(m)}=M r , like, for example, the harmonic excitation (N=1), for which h(m)=cos(ωt+ψ 1 ) and M r =M l =1. These last equalities and (13)-(14) further show that the curve ) . 4, represents the coinciding right and left homoclinic bifurcation thresholds in the case of harmonic excitation, which is considered as a reference to measure the improvement obtained with control in the next sections.

  Figure 4. The curves ) ( , 11 ω γ h cr

  have to solve the following mathematical problem of optimization:

  t)=γ 11 sin(πz)[sin(ωt)+0.3535

Figure 5 .

 5 Figure 5. The functions F(z,t)/γ 11 sin(πz) for various values of the parameter p= ) 2 /( 2 π ω -Γ and for N=3. a) optimal one-side "right" excitation, i.e., (17) 2 . Note that in this case there are two controlling superharmonics added to the reference basic harmonic. b) optimal global excitation, i.e., (19) 1 . In this case there is only one added superharmonic.

Figure 6 .

 6 Figure6. The first five normalized eigenfunctions w j (z) of the fixed-fixed case for Γ=40.Continuous lines are odd eigenfunctions, and dashed lines are even eigenfunctions.

  , and that it cubically approaches its tangent manifold at x=0. Thus, the nonflatness of Σ is expected to be engineering important only "far enough" from the saddle.

  Figure7. The slave amplitudes y 1 (x) and y 2 (x) (up to the 7 th order) for Γ=40, and a schematic section of a hypothetically bending exact manifold Σ.

  y i,25 and we have verified that in the range of interest x∈[0, 0.3] (see Sect. 6.2) the error is lesser than 0.005%, so that the functions of Fig.7are practically coincident with the true slave amplitudes y 1 (x) and y 2 (x).

  at checking the practical performances of the method.

Figure 1 .

 1 Figure 1. Examples of b.c. belonging to the first family. a) hinged-hinged: w(0)=w′′(0)=0, w(1)=w′′(1)=0, b) guided-hinged: w′(0)=w′′′(0)+Γw′(0)=0, w(1)=w′′(1)=0.

Figure 2 .

 2 Figure 2. Examples of b.c. belonging to the second family. a) fixed-fixed: w(0)=w′(0)=0, w(1)= w′(1)=0, b) fixed-hinged: w(0)=w′(0)=0, w(1)=w′′(1)=0, c) fixed-free: w(0)=w′(0)=0, w′′(1)=w′′′(1)+Γw′(1)=0.

Figure 3 .

 3 Figure3. A schematic picture of the dynamical system structure in the unperturbed case. The center manifold is perpendicular to Σ at the saddle.

FigureFigure 5 .

 5 Figure 4. The curves ) ( , 11 ω γ h cr

Figure 6 .

 6 Figure 6. The first five normalized eigenfunctions w(z) of the fixed-fixed case for Γ=40. Continuous lines are odd eigenfunctions, and dashed lines are even eigenfunctions.

The previous expression permits to compute an approximation of the homoclinic loop of the full dynamics, which can now be given the form w h (z,t)=x h (t)w 1 (z)+y 1 (x h (t))w 3 (z)+ y 2 (x h (t))w 5 (z)+… (39) accounting for the dynamics on the higher modes. According to the terminology of [START_REF] Lenci | Dimension reduction of homoclinic orbits of buckled beams via nonlinear normal modes technique[END_REF], we are now performing a refined order analysis of the homoclinic bifurcation.

Using this expression in the Melnikov function of the H.M. theorem we get, after some computations

where

Note that the constant β 0 and the functions β j (ω) do not actually depend on k. They are even functions of ω, and in the interval [0,+∞[ they are monotonically decreasing and exponentially approach zero for ω→∞.

Comparisons: flat vs non-flat manifolds, infinite vs finite dimensional models

Expressions (40) and (41) make explicit the contribution of non-flatness of Σ to the homoclinic bifurcation, which is taken into account by the terms containing y i (x), which are missing in (10).

There are contributions both in the constant part β 0 and in the oscillating part of the Melnikov functions, due to β 3 (ω), β 5 (ω),…, and measured by the coefficients γ n3 , γ n5 ,….

A c c e p t e d m a n u s c r i p t 24

To quantitatively measure these contributions, expressions (41) must be determined explicitly. β 0 and β 1 (ω) can be computed exactly, and in the reference case Γ=40 they are given by β 0 ={0.603839×10 -1 +0.428994×10 -6 +0.299111×10 -8 +…}

The coefficients β j (ω), j>2, on the other hand, can be evaluated only numerically. Therefore, supported by the fact that β j (ω) are monotonically decreasing in [0,+∞[, we consider only the maximum values β j (0) of each function β j (ω), which are given below (β 1 (0) is reported for comparison):

Once more, we have found that the non-flatness of Σ, while having theoretical relevance, gives contributions to the Melnikov functions which are orders of magnitude smaller than its flat approximation (given by the dominant terms and obtained by the one d.o.f. linear Galerkin approximation), and it is negligible from an engineering point of view for the considered parameters.

The previous considerations not only permit to conclude that the non-flatness of Σ can be neglected, but also that the bifurcation analysis based on the 1D reduced order model (accounting for only the first term of β 0 and for β 1 ) gives the "same" results as the infinite dimensional analysis.

This conclusion was obtained also in Sect. 5 for the case of flat Σ, but here, like there, it is important to stress that the full powerfulness of the H.M. theorem is necessary to take into account out of manifold perturbations, even if it has been shown a posteriori that they are practically negligible, a fact that couldn't easily be guessed a priori.

The difference is that for flat Σ the one d.o.f. model gives mathematically correct results in terms of Melnikov functions (Sect. 5), while for non-flat Σ it gives practically correct results, so that the theoretical background is different. However, it can be expected that, moving away from the first buckling load -which would entail considering a refined homoclinic orbit [START_REF] Lenci | Dimension reduction of homoclinic orbits of buckled beams via nonlinear normal modes technique[END_REF]] -, could evidence more marked differences of homoclinic bifurcation herein considered with respect to the linear Galerkin approximation.

Optimal control

Following the guidelines of Sect. 4, the application of the optimal control method to the second family of b.c. is now discussed.

The influence of the spatial shape of the excitation is considered first. As in Sect. 4.1, we assume In addition to the previous properties, and as in Sect. 4.1, the choice γ nj =0, j=3,5,…, permits to get rid of the non-resonance conditions (8) of hypothesis (H3), and thus contributes to enlarge the range of applicability of the method.

It is worth stressing the different meaning of the same assumption f n (z)=γ n1 w 1 (z) on the spatial shape of the excitation used for flat and non-flat manifold cases. In the former, it has no consequence in terms of the Melnikov function and only permits to reduce the cost of control (Sect.

4.1). In contrast, in the latter, it is mandatory in order to strongly simplify the Melnikov function (40) (see eq. ( 44)), and the successive analysis. In both cases, however, this choice does not arise from mathematical optimization problems, as the choice of the temporal shape, but just from practical considerations.

The previous assumption permits to rewrite the Melnikov function as

which is formally identical to (10). Then, the representation ( 11) is still valid with

which have the same meaning of Sect. 3.

This has strong consequences. The first is that the analysis of Sect. 4.2, where the optimal temporal shape of the excitation has been obtained, immediately applies to the present case.

Consequently, we refer to that Section for the determination of the excitation permitting the optimal shift of the homoclinic bifurcation in parameter space.

The second is that this result further confirms the generality of the proposed control method, which applies to a variety of mechanical systems and models within the same framework, as discussed in detail in [Lenci & Rega, 2004a].

To conclude this section, we note that the previous points have mainly a theoretical interest. In fact, according to the final considerations of Sect. 6.3, the application of the control method can be done as in the 1D Galerkin approximation, which is the subject of [Lenci & Rega, 2003b]. This

A c c e p t e d m a n u s c r i p t

26 would provide just the first term in the expression (42) 1 of β 0 , thus neglecting the contributions of higher order modes.

Extensions and generalizations

Various extensions and generalizations of the analysis of this paper are possible in principle and worthy to be investigated. In this section only two of them, the first of mechanical nature and the second of theoretical nature, are discussed.

Shallow arches

In addition to the buckled, initially straight beam, the control method can be applied to other mechanical systems. The application to shallow elastic arches, whose nonlinear dynamics is governed by the PDE equation [START_REF] Fung | Buckling of Low Arches or Curved Beams of Small Curvature[END_REF][START_REF] Lenci | Dynamics of shallow elastic arches. Part I: chaotic response of harmonically shaped arches[END_REF]]

is straightforward. In ( 46) y(z) represents the dimensionless undeformed shape of the arch, which can also be interpreted as the shape perturbation to a theoretical straight beam, while the other symbols have the same meaning of the previous sections.

First we note that if the arch is hinged at both ends and the initial spatial shape is harmonic, i.e., y(z)=hsin(πz) then the invariant manifold Σ of the theorem is flat and we are in the first family (see Sect. 2), which now involves not only boundary conditions but also undeformed shapes. The nonlinear dynamics and the homoclinic bifurcation of this case with harmonic (in time) excitation have been investigated in [START_REF] Lenci | Dynamics of shallow elastic arches. Part I: chaotic response of harmonically shaped arches[END_REF] by the H.M. theorem. The single d.o.f. model obtained by the Galerkin method was also obtained, and it was shown to provide the same results of the infinite dimensional analysis, according to the fact that the manifold Σ is planar.

The extension of the analysis of [START_REF] Lenci | Dynamics of shallow elastic arches. Part I: chaotic response of harmonically shaped arches[END_REF] to the non-harmonic, but still periodic in time, excitations, which is necessary for the application of control, can be obtained by superposition, just as done in Sect. 3, and will not be detailed. We only note that the sag h of the arch plays the role of the axial load Γ of the buckled beam.

The main mechanical difference between the shallow arch and the buckled beam is that the non vanishing initial shape destroys the symmetry of the system, which accordingly has also quadratic nonlinearities (see eq. ( 46)). This is at the base of the main dynamical differences, which are the non coincidence between the right and left homoclinic loops [START_REF] Lenci | Dynamics of shallow elastic arches. Part I: chaotic response of harmonically shaped arches[END_REF], eq. ( 7)],

which in turn implies different values of the homoclinic bifurcation threhsolds [Lenci & Tarantino, 1996, eq. (10)].

The quadratic terms is present also in the reduced order model, and the resulting equation is the

Appendix

In this appendix some results used in the text are reported.

Tables A1 andA2 have been obtained in [Lenci & Rega, 2003b, 2004a]. The normalized eigenfunctions of Sect. 6.1, depicted in Fig. 6, are analytically given by: w 1 (z)=-0.011200sin(6.31091z)-0.80798cos(6.31091z)+0.17021sin(0.41526z)+0.80798cos(0.41526z);

w 2 (z)=0.72180sin(8.27104z)-1.10927cos(8.27104z)-1.12007sinh(5.33011z)+1.10927cosh(5.33011z);

w 3 (z)=0.88169sin(11.18713z)-1.06911cos(11.18713z)-1.06890sinh(9.22778z)+1.06911cosh(9.22778z); (A1) w 4 (z)=0.93579sin(14.24669z)-1.04433cos(14.24669z)-1.04434sinh(12.76590z)+1.04433cosh(12.76590z);

w 5 (z)=0.95954sin(17.34999z)-1.03044cos(17.34999z)-1.03044sinh(16.15618z)+1.03044cosh(16.15618z);