Selective Pyragas control of Hamiltonian systems

We consider a Newton system which has a branch (surface) of neutrally stable periodic orbits. We discuss sufficient conditions which allow arbitrarily small delayed Pyragas control to make one selected cycle asymptotically stable. In the case of small amplitude periodic solutions we give conditions in terms of the asymptotic expansion of the right hand side, while in the case of larger cycles we frame the conditions in terms of the Floquet modes of the target orbit as a solution of the uncontrolled system.


1.
Introduction. In non-linear sciences stabilization of periodic solutions is an important problem. In [13] Pyragas suggested a simple method to stabilize a periodic solution of a system by introducing feedback of the form K(x(t − T ) − x(t)) where T is the period of the solution that is to be stabilized. An important feature of this control is that it is non-invasive in the sense that it does not change the periodic solution itself but only its stability properties. This method has been implemented in a wide variety of applications, see e.g. [14,16,17,18].
In most applications Pyragas control is used to stabilize an unstable orbit which is hyperbolic. That is to say that at least one of its Floquet multipliers is outside the unit circle, and the gain matrix K must be sufficiently large to 'move' these unstable eigenvalues inside the unit circle. On the other hand, since introducing feedback transforms the system to a delay differential equation, for small K infinitely many eigenvalues appear near zero. However, with the increasing gain these additional eigenvalues may become difficult to control. Therefore a delicate balance must be maintained: too small gain and the original unstable eigenvalues cannot be controlled, too large gain and the the additional eigenvalues may become unstable. For a discussion of sufficient conditions under which small orbits born near a Hopf bifurcation point can be stabilized see [9]. For necessary conditions for stabilization of generic orbits see [8]. In [4,6,7,15] these results were extended to equivariant systems.
In this paper we extend the Pyragas control scheme to Hamiltonian systems. Due to the additional symmetry inherited from the Hamiltonian structure, Floquet multipliers come in pairs, µ, µ −1 . For this reason, periodic solutions are either unstable or neutrally stable. In this work we concentrate on the case of neutrally stable solutions. The multiplier µ = 1 will always have multiplicity at least 2, and it is generic for periodic orbits to foliate a surface in the phase space. However, generically different periodic solutions lying on the same surface have different periods and hence proper selection of the parameter T in the Pyragas control scheme can select between them. For convenience in this paper we concentrate on Newton equations of the typeẍ + ∇V (x) = 0. Our goal is to transform a neutrally stable periodic orbit of the uncontrolled system to the exponentially stable orbit by implementing feedback of the form K(x(t − T ) − x(t)) where T is the period of the targeted orbit to be stabilized and K is an arbitrarily small gain matrix. This is in contrast to the situation which is generic in non-Hamiltonian systems where control can never be arbitrarily small. We discuss sufficient conditions from two perspectives.
In Section 2 we treat the case of small amplitude solutions where sufficient conditions for exponential stability are framed in terms of the asymptotic expansion of ∇V (x) at the equilibrium x = 0 up to the third order. The importance of the third order expansion is motivated by the Lyapunov centre theorem and the normal form for the Hopf bifurcation. Section 3 deals with arbitrary periodic solutions, but instead of being framed in terms of the asymptotic expansion of the field ∇V (x), we use conditions on the Floquet modes of the targeted orbit as a solution of the uncontrolled system.
Each section is supplemented with examples. In a brief Section 4 we show how the conditions in Section 3 agree with those from Section 2 in the case of small amplitude cycles. Finally we present conclusions. An Appendix containing the derivation of the normal form for the Hopf bifurcation in a delayed system is included at the end of the paper for convenience of the reader.
2. Stabilization of small periodic orbits.
2.1. Main statement. Consider the system with a sufficiently smooth potential V , where x τ (t) = x(t − τ ). Let x * be an equilibrium of the uncontrolled system Without loss of generality, assume that x * = 0, i.e. ∇V (0) = 0. Denote by H = D 2 V the Hessian of V and by H 0 the value of the Hessian at zero, H 0 = D 2 V (0). We will use the third order expansion of ∇V at zero: with a bilinear map L : R N × R N → R N and a trilinear map Q : Here and henceforth, (·, ·) is the usual scalar product and | · | is the corresponding Eucledian norm in R N . Denote by Sp(H 0 ) the spectrum of the symmetric matrix H 0 . We make the following assumptions about the uncontrolled system (2).
• Stability of the zero equilibrium: The Hessian H 0 is positive definite.
Hence, the spectrum of H 0 is positive: • Non-resonance condition: For some k and all n = k, 1 ≤ n ≤ N , the ratio ω n /ω k is not an integer number. Further ω n = ω m for n = m. Denote by e n a unit eigenvector of the Hessian H 0 corresponding to the eigenvalue ω 2 n and set where I denotes the identity matrix.
• Non-zero Lyapunov quantity: By the Lyapunov Center Theorem, the non-resonance condition ensures that the uncontrolled system (2) has a continuous branch (one-parameter family) Γ k of small periodic orbits with the (minimal) periods satisfying T → 2π/ω k as x C → 0 (e.g. [11]). We note that if the non-resonance condition is satisfied for several ω k , then there are several branches of small periodic orbits near zero. In particular, if no frequency is a multiple of any other frequency, then there are n branches of small periodic orbits with pairwise different asymptotic periods.
The proof of the theorem uses the following simple statement. Lemma 2.2. Let A and B be two N × N matrices and 0 < ε 1. Assume that A is diagonal. Then, Proof. Consider the function g(ε) = det (A + εB). By Taylor's formula, where by the chain rule with C(ε) := A + εB.
Here where E ij is the N × N matrix with 1 at the i-th row and the j-th column, and zero everywhere else. Recall that A is diagonal, therefore Hence, Combining this with (9) and (10) gives (8).
Proof of Theorem 2.1. Denote by S an orthogonal matrix that diagonalizes H 0 : The characteristic equation of the linearization of the controlled system (1) with K = εD at zero can be written as whereD = S T DS. For ε = 0, the roots of this equation are ±iω n , n = 1, . . . , N . Let us set τ = 2π/ω k , take n = k, and consider a small perturbation λ = iω n + ρ n of the root iω n for 0 < ε 1. Substituting the expression λ n = iω n + ρ n in (11) and using Lemma 2.2, we see that the zero order terms vanish, and keeping the first correction terms, we obtain 2iω n ρ n + ε(1 − e −2πiωn/ω k )D nn + o(ε + |ρ n |) = 0, whereD nn = e n ·De n is the corresponding diagonal entry of the matrixD. This implies Re λ n = Re ρ n = −e n · De n sin(2πω n /ω k ) 2ω n ε + o(ε), n = k.
Now, we fix an ε > 0 and note that the characteristic equation (11) has the roots ±iω k for τ = 2π/ω k . Considering a small perturbation τ = 2π/ω k − ξ k δ of the delay, for the perturbed root λ k = iω k + ρ k of (11) we obtain where δ is considered as a small parameter and ε is fixed. Hence, Relations (6), (7) combined with (12), (13) ensure that the controlled system (1) with the delay τ considered as the bifurcation parameter and K = εD satisfies the conditions of the Hopf Bifurcation Theorem (e.g. [5]) for any sufficiently small fixed ε > 0. Indeed, relations (7) and (13) imply that the complex conjugate eigenvalues λ k (τ ),λ k (τ ) cross the imaginary axis transversally at the points ±iω k for τ * = 2π/ω k . At the same time, relations (6) and (12) imply that all the eigenvalues λ n (τ ),λ n (τ ) with n = k, 1 ≤ n ≤ N have negative real parts for any τ close to τ * , while the smallness of ε ensures the same property for all the other eigenvalues. Hence, equation (1) has a continuous branch (one-parameter family)Γ of small periodic orbits with (minimal) periods satisfying T → 2π/ω k and τ → 2π/ω k as x C → 0. Further, stability of these small periodic orbits is determined by the Lyapunov quantity [3]. A computation of this quantity is included in the Appendix for convenience of the reader. In particular, it is shown that condition (7) ensures that the Hopf bifurcation is supercritical, i.e. all the small orbits of the branchΓ are asymptotically stable. Finally, the Hopf Bifurcation Theorem implies that all the small periodic orbits of (1) that exist for τ close to τ * and have periods close to τ * belong toΓ (e.g. [5]), in particular the branch Γ k of small periodic orbits of the uncontrolled system (2) is included inΓ. 1 This completes the proof. Remark 1. An asymtotic analysis presented in Appendix shows that the period of a small periodic orbit of the uncontrolled system (2) scales with the amplitude of the orbit as Combining this with (7) and assuming that (6) holds, we see that if |T − τ * | scales quadratically with the amplitude of the orbit and T decreases with x C , then the stabilization of small periodic orbits is achieved when e k · Ke k > 0. On the other hand, if T increases with x C , then the stabilization is achieved when e k · Ke k < 0.

2.2.
Example: Duffing oscillator. Consider the scalar equation Without the control (κ = 0), the scaling of the period of the periodic orbit near the origin can be found, for example, by the Poincaré-Lindstedt method [10]: where prime denotes the derivative with respect to the rescaled timet = ωt and where r measures the amplitude of the small periodic solution. The first correction of order r 3 has no secular terms for A = 3/8. Hence, the frequency increases and the period T decreases with the amplitude r of the periodic orbit: According to Remark 1, the condition κ > 0 ensures the non-invasive stabilization of small periodic orbits provided that κ is sufficiently small and τ < 2π is close to 2π. Similarly, a sufficiently small control with κ < 0 and τ > 2π, τ ≈ 2π stabilizes a small periodic orbit of the equation

2.3.
Example: Two coupled Duffing oscillators. Now, let us consider two coupled Duffing oscillators: with the potential Here the eigenvalues and eigenvectors of the Hessian H 0 are Using the diagonalizing coordinate transformation Two branches of periodic orbits can be obtained by the Poincaré-Lindstedt method using the following expansions: with r 1. Elimination of the secular terms in the first correction is achieved for for the first branch with ω ≈ √ 2 and for the second branch with ω ≈ √ 0.7. Since γ 1,2 > 0, the frequency increases and the period of the periodic orbit decreases on each branch with the amplitude r of the orbit. According to (14), it means that the Lyapunov quantity (5) for each branch satisfies ξ 1,2 > 0.
Let us assume a diagonal control matrix Using (16) and ξ 1 > 0, the conditions (6), (7) for the branch of small periodic orbits with ω ≈ √ 2 have the form According to Theorem 2.1, these conditions ensure that the branch with ω ≈ √ 2 stabilizes provided that κ 1 , κ 2 are sufficiently small in absolute value.
Similarly, the branch of small periodic orbits with ω ≈ √ 0.7 stabilizes if the opposite inequalities hold: On the other hand, if then each of the branches will be stabilized if 3. Stabilization of large periodic orbits.

Main statement.
We consider now a branch of periodic orbits of the uncontrolled system (2) further away from zero. Consider a particular periodic orbit Since the orbit is included in the branch, it has a characteristic multiplier µ = 1 of algebraic multiplicity greater than 1. We make the following assumption.
• Genericity assumption: The characteristic multiplier µ = 1 of the periodic orbit x * has geometric multiplicity 1 and algebraic multiplicity 2. This is the generic situation for a periodic orbit embedded into a surface of such orbits. More precisely, we write the uncontrolled system in the phase space: consider the periodic solution (x * (t), p * (t)) = (x * (t), −ẋ * (t)), and consider the linearization of system (19) on the periodic solution: with y = (x, p) and where H(x * (t)) = H T (x * (t)). Hence, B(t) ≡ B(t+T * ). Floquet modes are solutions of (20) of the form where µ n is a (complex) characteristic multiplier. In particular, denote u 0 (t) = x * (t), then y 0 (t) = (u 0 (t), −u 0 (t)) is the periodic Floquet mode, which satisfies and has the characteristic multiplier µ 0 = 1. Also, according to the genericity assumption above, equation (20) has a solution (v 0 (t), −v 0 (t)) satisfying which can be called the generalized Floquet mode. Further, assume that • Stability assumption: All the characteristic multipliers µ n = 1 of (20) are simple and satisfy |µ n | = 1. In other words, the periodic solution x * (t) of the uncontrolled system is stable.
Let µ denote the complex conjugate of µ ∈ C and u · v = u T v for u, v ∈ C N .
Then, for every sufficiently small ε > 0, the periodic solution x * (t) of the controlled system (1) with τ = T * and K = εD is orbitally asymptotically stable.
Note that in (23), (24), be the linearization of the controlled system (1) with the delay τ = T * and the gain matrix εG near x * . Suppose that µ is a characteristic Floquet multiplier of system (25), hence this system has a solution (Floquet mode) with a periodic q(t) ≡ q(t + T * ). Since the delay in (25) equals the period of q, we see that y µ satisfies the ordinary differential systeṁ Denoting by Ψ ε,µ (t) the fundamental matrix of (27), we conclude that the monodromy matrix Ψ ε,µ (T * ) of (27) also has the characteristic multiplier µ. Hence, this multiplier µ is a root of the characteristic equation Denote by Φ(t) the fundamental matrix of system (20). Clearly, Φ(t) = Ψ 0,µ (t). Suppose that the monodromy matrix Φ(T * ) of system (20) has a simple eigenvalue µ * = 1 on the unit circle, hence µ * is a root of (28) for ε = 0. Let us consider the perturbation µ = µ * + ρ of this root for small ε. Note that the fundamental matrices Φ(t) and Ψ ε,µ (t) are related by the identity where we omit the terms of order o(ε). In the following, the approximate equality also means that o(ε) terms are omitted. Let us denote by S the transition matrix to a basis in which the matrix Φ(T * ) assumes the Jordan form Λ; for convenience, we also agree that S −1 maps the eigenvector e * corresponding to the simple eigenvalue µ * to the vector S −1 e * = e 1 := (1, 0, . . . , 0) T ∈ R 2N . In this new basis, using (29), we can rewrite (28) as with µ = µ * + ρ. Since Λ 11 = µ * , this implies Taking into account that y * (t) = Φ(t)e * = Φ(t)Se 1 is the Floquet mode of (20) corresponding to the characteristic multiplier µ * and ϕ † T e 1 is the adjoint Floquet mode corresponding to the multiplier 1/µ * = µ * and normalized by the condition ϕ † * (t) · y * (t) ≡ 1, we see that Therefore, the stabilization condition |µ * +ρ| < 1, which is equivalent to Re (ρµ * ) < 0 for small ρ, is Now we consider the perturbation µ = 1 + ρ of the eigenvalue 1. By assumption, the Jordan form Λ of the matrix Φ(T * ) has the Jordan block Further, y * (t) = Φ(t)e * = Φ(t)Se 1 is the Floquet mode of (20) corresponding to the multiplier 1 and ϕ † T e 2 is the adjoint Floquet mode. Therefore and the stabilization condition reads Notice that ϕ † * (t) · y * (t) ≡ 0 automatically, and relations (32), (33) imply the following normalization condition in terms of z * = Φ(t)g * : Thus, conditions (31) and (34) ensure the stabilization. Finally, let us rewrite this conditions using the specific structure of our system and the specific structure of the control matrix G. First, note that due to the block structure of the matrix (21), each eigenfunction has the form y * (t) = (u * (t), −u * (t)) where u * has values in R N . Further, we see that if an eigenfunction y * (t) = (u * (t), −u * (t)) corresponds to an eigenvalue µ, then ϕ † * (t) = α(u * (t), u * (t)) is an adjoint eigenfunction corresponding to the same eigenvalue for any α = 0. For µ = µ * = 1, these eigenfunctions should satisfy the normalizing condition .

3.2.
Example. Again, we consider two coupled Duffing oscillators with the potential (15). We have already established that there are two branches of periodic solutions emanating from the origin. Close to the origin the frequencies of these branches are ω ≈ √ 2 and ω ≈ √ 0.7. Moreover we have shown that small periodic orbits can be stabilized by the diagonal control matrix K = diag (κ 1 , κ 2 ) satisfying conditions (17) and (18) for ω = √ 2 and ω = √ 0.7, respectively. Using numerical continuation with the delay τ as a continuation parameter we can "travel" along the branch of interest and study the stability of the periodic solutions away from the origin. The period of the orbit equals the delay, T = τ . Figure 1 presents the branch of periodic orbits starting from the zero equilibrium at τ = 2π/ √ 0.7. The period decreases and the amplitude increases along the branch. Depending on the value of κ = (κ 1 , κ 2 ), different parts of the branch are stable/unstable. Panels (A) and (B) present the stability of the solutions along the branch forκ = (−0.02, 0.012) andκ = (−0.02, 0.02), respectively. Bothκ andκ satisfy conditions (18), therefore the small-amplitude solutions are asymptotically stable. For κ =κ the orbit destabilizes at the point B (see Figure 1(A)), while for κ =κ all the orbits remain stable for all the values of τ considered (see Figure  1(B)). Continuation of the branch was performed using DDE-BIFTOOL package for MATLAB ® [1, 2]. Conditions (23), (24) were evaluated numerically for three points A, B and C of the branch, see Figure 1. For each point, these conditions define a sector with the vertex at the origin in the plane of control parameters (κ 1 , κ 2 ). The three corresponding sectors, A 1 OA 2 , B 1 OB 2 , C 1 OC 2 , are shown in Figure 2. In particular, if κ belongs to the interior of the sector A 1 OA 2 and |κ| is sufficiently small, then Theorem 3.1 ensures that the periodic orbit corresponding to the point A is asymptotically stable for the controlled system when the delay is set to τ = τ A . Figure 2 suggests that the sector becomes more narrow with the decreasing τ . Correspondingly, sinceκ belongs to C 1 OC 2 , the stable part of the branch extends from its origin at least to the point C for the controlled system with K = diag(κ 1 ,κ 2 ). On the other hand, for K = diag(κ 1 ,κ 2 ), the stable part of the branch extends only to the point B becauseκ lies on the boundary of the sector B 1 OB 2 . This agrees with Figure 1. 4. How do stability conditions for small and large cycles agree? For a small cycle, the Hessian is close to its constant value at zero. Therefore the Floquet modes are close to the eigenfunctions at zero. In particular, for n = k, we have u * (t) ≈ e iωnt e n . Using this approximation, and we see that (24) is equivalent to (6). On the other hand, the solution is approximately x * (t) ≈ re k cos(ω k t) and for n = k, we have u * (t) =ẋ * (t) ≈ −rω k e k sin(ω k t). If we assume in addition to (5) that the period T is a strictly monotone function of the amplitude x C of the periodic orbit along the branch (for small cycles), and the derivative ∂T /∂ x C is well-defined and is non-zero, then the sign of this derivative is opposite to the sign of ξ k (cf. (14)). To be definite, assume that ξ k > 0, hence T decreases with x C along the branch. Then, v * (t) ≈re k cos(ω k t) whererr > 0. Therefore, ≈ ω 2 k r 2 sin 2 (ω k t) e k · De k , and we see that (23) implies e k · De k > 0, which agrees with (7). In the case ξ k < 0, we have v * (t) ≈re k cos(ω k t) withrr < 0, which leads to (7) again.

5.
Conclusions. We considered an extension of the Pyragas delayed control method from general differential equations to Newton systems that possess a branch of neutrally stable periodic orbits. We proposed sufficient conditions which allow arbitrarily small Pyragas control to transform the stability of one selected periodic solution from neutrally stable to exponentially stable. In the case of small amplitude periodic solutions our conditions were given in terms of the asymptotic expansion of the vector field, while in the case of general periodic orbits we expressed the conditions in terms of the Floquet modes of the target orbit as a solution of the uncontrolled system. The theorem were illustrated by analytical and numerical examples. The results can be naturally extended to Newton systems with controls K(ẋ τ −ẋ) and to Hamiltonian systemsẋ = J∇H(x) with controls K(x τ − x).
Appendix: Computation of the normal form. Let us consider the delay τ > 0 in equation (1) as a bifurcation parameter, which is varied near the point τ * = 2π/ω k . The characteristic equation of the linearization of (1) at zero reads Assume that conditions (5), (6) hold and K = εD with a sufficiently small ε > 0. Then, as shown in the proof of Theorem 2.1, a pair of complex conjugate eigenvalues λ k (τ ),λ k (τ ) crosses the imaginary axis transversally at the points ±iω k for τ * = 2π/ω k , while all the other eigenvalues have negative real parts. We now use the multiple time scale method to compute the normal form of system (1) near the zero equilibrium x = 0 for τ close to τ * . According to (3), we have the following third order expansion of equation (1): We seek a uniform second-order approximate periodic solution in the form where T 0 = t, T 2 = s 2 t, and s is a nondimensional book keeping parameter; h.o.t. stays for higher order terms in s. Note that it is not necessary to relate the approximate solution with the time scale T 1 = st because a second order approximation will require independence of T 1 in the solution (e.g. [12], page 166). Using the notation for the time derivatives, from (39) one obtains Further, introducing the notation τ = τ * + s 2 δ where δ is the detuning parameter, Now, substituting (39)-(41) into (38) and comparing the like powers of s up to the cubic terms yields the following equations, where for notational simplicity we omit the time scales of the variables except for the delayed ones: Since τ = τ * is a Hopf bifurcation point for (38), the linearization (42) has a periodic solution of the form where the real eigenvector e k corresponds to the eigenvalues ±iω k of (37) for τ = τ * = 2π/ω k . Substituting this expression in equation (43) and using e iω k τ * = 1, we obtain a particular solution with a k , b k defined by (4). Further, substituting the expressions for x 1 and x 2 into (44), we obtain −Ā e −iw k T0 e k + L (Ae iw k T0 +Āe −iw k T0 )e k , a k (A 2 e 2iω k T0 +Ā 2 e −2iw k T0 ) + b k AĀ + 1 6 Q(e k , e k , e k ) A 3 e 3iω k T0 +Ā 3 e −3iω k T0 + 3A 2Ā e iω k T0 + 3AĀ 2 e −iω k T0 = − τ * A e iw k T0 + τ * Ā e −iw k T0 + iω k δAe iw k T0 − iω k δĀe −iw k T0 Ke k .
Note that the nonhomogeneous equation (46) has terms proportional to e iw k T0 , which sum to χ(T 2 )e iw k T0 with χ = 2iω k A e k + L(e k , a k + b k ) + 1 2 Q(e k , e k , e k ) A 2Ā + τ * A Ke k + iω k δAKe k .
At the same time, iw k is an eigenvalue of the associated homogeneous equation. Since x 3 is a periodic solution of (46), it does not contain secular terms, which is only possible if the equation is solvable with respect to φ = φ(T 2 ), where φ(T 2 )e iw k T0 is the first harmonic in the Fourier expansion of x 3 . The solvability requires the orthogonality of χ with the eigenvector e k of H 0 corresponding to the eigenvalue ω 2 k : e k · Re χ = e k · Im χ = 0.
Substituting the above expression for χ in these relations leads to the following equation A + iω k δ e k · Ke k 2iω k + τ * e k · Ke k A + e k · L(e k , a k + b k ) + 1 2 Q(e k , e k , e k ) 2iω k + τ * e k · Ke k A 2Ā = 0, which is the normal form of the Hopf bifurcation for equation (1) with the parameter τ . Using the notation (5), this is equivalent to A + iω k δ e k · Ke k 2iω k + τ * e k · Ke k A + ξ k 2iω k + τ * e k · Ke k A 2Ā = 0.
The normal form (47) has a periodic orbit iff δξ k < 0. This orbit has the simple form A = A 0 e iω0t (a relative equilibrium). Computing A 0 , ω 0 from (47) and using the relations (39) and (45), one obtains the asymptotic relationship between the amplitude x C of a small periodic orbit of (1) and its period T , which is equal to the delay: Furthermore, if δe k · Ke k < 0, then the periodic orbit of the normal form (47) is asymptotically stable and so is the periodic orbit x(t) of (1). Due to δξ k < 0, this asymptotic stability condition is equivalent to (7) yielding the conclusion of Theorem 2.1. Finally, equation (47) with K = 0 provides asymptotics of small periodic solutions of the uncontrolled system (2). In particular, the relationship (48) between the period and the amplitude holds, which justifies Remark 1.