Low-thrust displaced orbits by weak Hamiltonian-Structure-Preserving control

Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

Contents lists available at ScienceDirect

Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns

Research paper

Low-thrust displaced orbits by weak Hamiltonian-Structure-Preserving control Xiao Pan, Ming Xu∗, Yunfeng Dong School of Astronautics, Beihang University, Beijing 100191, China

a r t i c l e

i n f o

Article history: Received 5 May 2017 Revised 21 October 2017 Accepted 26 February 2018 Available online 27 February 2018 Keywords: Low-thrust displaced orbits Weak hamiltonian-structure-preserving control Homoclinic orbits Transit orbits

a b s t r a c t For the requirements from some missions, the spacecraft is placed near the unstable region around the hyperbolic equilibrium. To stabilize the motion, a novel concept of the weak Hamiltonian-Structure-Preserving (HSP) control parameterized by the Coriolis acceleration is introduced. Based on the dynamical model of low-thrust displaced orbits, this paper searches for the mechanism in orbital stability change by the weak HSP control and provides specific strategies on how the control is implemented. Since the weak control has no effect on the topology of the original system, two invariant equilibria, a hyperbolic one and an elliptic one are solved. Dynamical system techniques are employed to investigate the controlled motions near the two equilibria, illustrating the Lyapunov orbit in different Coriolis acceleration cases and presenting two basic periodic modes measuring the bounded trajectories in the interior region. Controlled trajectories in the neck region are analyzed for their characteristics, classification and transition to make preparation for the stability change. One of the important contributions of this paper is to numerically demonstrate that the weak HSP control preserves the homoclinic orbits to the hyperbolic equilibrium as well as to the Lyapunov orbit, and another is to achieve the orbital transfer within, or even beyond KAM tori based on two methods in determining whether the controlled trajectory integrated from an initial point can exhibit transit. © 2018 Elsevier B.V. All rights reserved.

1. Introduction The concept of counter-acting gravity through a thrust vector was first proposed by Dusek in 1966 [1], who discovered that a spacecraft could be kept in a stable orbit around the collinear points through the use of continuous propulsive thrust. Large families of displaced orbits for a generic low-thrust propulsion were generated by considering the dynamics of the two-body problem in a rotating frame of reference [2,3]. The displaced orbits have been identified by solar sail or electronic thrusters, including three basic types of circular orbits by a thrust-induced acceleration [2–4], body-fixed hovering orbits by open-loop control [5], quasi-periodic displaced trajectories by a fixed thrust along the rotation axis of planet [6], a sequence of individual Keplerian arcs connected by slight impulse propulsion [7], elliptic displaced orbits with advanced thrust model [8], displaced geostationary orbit using hybrid sail propulsion [9], and so on. Such displaced non-Keplerian orbits (NKOs) have a diverse range of potential applications in both Sun-centered and planet-centered studies, such as the solar physics and an Earth synchronous orbit for continuous observations and space weather monitoring [10], displaced orbits above

∗

Corresponding author. E-mail addresses: [email protected] (X. Pan), [email protected] (M. Xu), [email protected] (Y. Dong).

https://doi.org/10.1016/j.cnsns.2018.02.036 1007-5704/© 2018 Elsevier B.V. All rights reserved.

X. Pan et al. / Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

283

Nomenclature

ω ρ z hz

κ

E U n I J

factor quantifying Coriolis acceleration displaced orbit radius displaced orbit height angular momentum thrust acceleration energy of the system potential function thrust direction identity matrix symplectic matrix

Earth-Moon L2 for lunar far-side communication [11], and a new family of NKOs displaced above or below the Earth’s equatorial plane for the increasing number of available slots for geostationary communications satellites [12]. However, all of the mission orbits near hyperbolic equilibria are unstable, so that linear and nonlinear strategies have been investigated by a number of authors for the station-keeping requirements. McInnes [13] presented a passive control scheme of solar sail by designing the configuration to render the unstable subset of orbits linearly stable. Bookless [14] made analysis for the dynamics and control of the motions near the equilibrium, but his many conclusions are linear and local. Scheeres et al. [15] constructed a powerful tool Hamiltonian-Structure-Preserving (HSP) control that can stabilize the system in the sense of Lyapunov even in nonlinear situations. Bookless and McInnes [16] first developed the linear quadratic regulator technique in solar sail applications, which was then applied by Qian [17] to a nonlinear dynamical system of displaced orbits with a relatively favorable control precision. Benefiting from the non-dissipative structure [18], the control of preserving Hamiltonian structure is preferable to be developed to stabilize the motions near the hyperbolic equilibrium. To extend the work of Scheeres, Xu and Xu [19] separated the HSP controller into two parts: the strong HSP characterized by the manifolds (stable, unstable, and center) and position feedback, and the weak HSP by the Coriolis acceleration and velocity feedback. They proved that a 2-degree-of-freedom Hamiltonian system can be stabilized by the strong HSP. Inspired by Xu and Xu [19], we focus on the investigation of the weak HSP control, where the first aim of our research is to investigate the influence of such station-keeping strategy on the dynamical behavior of motions near equilibria. Then, the second aim is to present how this control is implemented to stabilize the unstable motion. In this paper, the displaced non-Keplerian orbits of low-thrust spacecraft are controlled using different Coriolis acceleration and velocity feedback. The developed control scheme is named as weak HSP control, for it presents the weak Hamiltonian-Structure-Preserving characteristics. Similarly, dynamical system techniques are used to analyze the nonlinear dynamics of this displaced orbit in a controlled system. Firstly, the weak controller is constructed and applied to a twobody problem, and two of the hyperbolic and elliptic equilibria which are independent of the control are solved. Then, the influence of the weak HSP control on the motions near the two equilibria is investigated by means of dynamical system techniques. The interior region is presented to be filled with the bounded trajectories, which can be measured by two basic periodic modes, referred to as the first type of periodic orbit (denoted as P.O.I) and the second type of periodic orbit (denoted as P.O.II). Subsequently, numerical simulations demonstrate that the stable and unstable manifolds in the interior region still connect homoclinically with each other. Next, the trajectories in the neck region are analyzed for their characteristics, classification and transition to make preparations for the stabilization of unstable orbits near hyperbolic equilibrium. Finally, the orbital transfer between different bounded orbits as well as the stable and unstable orbits is exemplified to show the function and application of the weak HSP control implemented. 2. Nonlinear dynamics of controlled displaced orbit 2.1. Hamiltonian dynamics of displaced orbit by low thrust In general, displaced NKOs can be achieved by seeking equilibrium solutions to the two-body problem in a rotating frame of reference R, as shown in Fig. 1, where the frame of reference R rotates with constant angular velocity ω relative to an inertial frame I. It is assumed that the spacecraft of mass m at position r has the active propulsion generating thrust T in direction n, and the magnitude of the thrust-induced acceleration is constant. This paper considers the two-body dynamics of a spacecraft by ignoring the higher harmonics of the gravitational potential and then generates NKOs displaced above the Earth using the low-thrust propulsion. To simplify,  the system is nondimensionalized through a characteristic length Re, i.e., the Earth radius, and a characteristic time T = Re3 /μ, where μ is the gravitational parameter of the planet chosen to be unity, so that the nonlinear dynamical equations of motion for a spacecraft in the polar coordinates (ρ , z, θ ) are given by [14]:

ρ¨ = h2z /ρ 3 − ρ /r3

(1a)

284

X. Pan et al. / Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

Fig. 1. Displaced non-Keplerian orbit reference frames.

z¨ = −z/r 3 + κ

(1b)

θ¨ = −2ρ˙ θ˙ /ρ

(1c)

where  ρ (ρ ≥ 0) denotes the orbit radius, z is the displacement distance parallel to the Sun line, θ is the azimuthal angle, r = ρ 2 + z2 is the distance between the central body and spacecraft, hz represents the constant angular momentum hz = ρ 2 θ˙ , and κ is the scaled thrust acceleration, in this paper it is fixed as κ = 1.0545 × 10−4 [14].

According to Eq. (1), the system can be treated as a two and half degree of freedom (DOF) system because θ has no effect on ρ and z. Therefore, the system of the spacecraft can be further simplified as a plane problem in which only ρ and z are used as variables to rewrite the dynamical model as follows:

⎧ ⎨ρ¨ = − ∂ U ∂ρ ⎩z¨ = − ∂ U ∂z

(2)

where the potential function U is

U=

1 h2z − − κ · z. r 2ρ 2

(3)

The energy of this system can be expressed as

E=

1 2 1 2 ρ˙ + z˙ + U. 2 2

(4)

Let q = [ρ , z]T . The equilibria for the system (2) can be achieved by



q = q0 , q˙ = 0 . ∇U = 0

(5)

Its variation equation near the equilibria is

δ q¨ + Uqq δ q = 0

(6)

where Uqq is the second derivative of U with respect to q. Based on the eigenvalues of Uqq , the stability of the equilibria can be classified into two categories [6]: the elliptic equilibrium if Uqq has two negative eigenvalues, and the hyperbolic equilibrium if Uqq has one negative and one positive eigenvalue. By choosing the angular momentum hz = 6.6285, these two equilibria, a hyperbolic equilibrium (saddle) denoted by Lu and an elliptic equilibrium (center) denoted by Ls , are solved for the system. As specified by Xu and Xu [6], Eq. (2) describes an energy surface achieved by setting the system energy as a constant E0 , and the energy surface is defined as

(κ , hz , E0 ) = { (ρ , z, ρ˙ , z˙ )|E (ρ , z, ρ˙ , z˙ ) = E0 }.

(7)

X. Pan et al. / Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

285

Fig. 2. Four basic configurations of Hill’s region.

where  is the invariant three-dimensional manifold. Then, Hill’s region [20] is defined as the projection of the energy surface onto the position space

R(κ , hz , E0 ) = { (ρ , z )|U (ρ , z ) ≤ E0 },

(8)

where the boundary of R(κ , hz , E0 ) is termed the zero-velocity surface. For different system energies, there are four basic configurations of Hill’s region, as shown in Fig. 2. Certain terminologies referred to later are defined as follows: when E ≤ U(Lu ), the allowed bounded region surrounding Ls is termed Island, which can be formalized as I (κ , hz , E0 ) = {(ρ , z )|U (ρ , z ) ≤ U (Lu ), z ≤ zu }, where zu is the z coordinate component of hyperbolic equilibrium Lu ; and Mainland is defined as the allowed unbounded region extending to infinity, formalized as M (κ , hz , E0 ) = {(ρ , z )|U (ρ , z ) ≤ U (Lu ), z > zu }; when E > U(Lu ), a Neck Region emerges, and the allowed region is termed Byland, which can be formalized as B(κ , hz , E0 ) = {(ρ , z )|U (Lu ) < U (ρ , z ) ≤ E0 }. 2.2. Weak Hamiltonian-Structure-Preserving control for low-thrust displaced orbits Most astrodynamical problems (e.g., planar circular three-body problem, and planar solar sail three-body problem, our problem presented in this paper) can be classified as the Hamiltonian system. For the requirements from some missions, the spacecraft placed near the unstable region around the hyperbolic equilibrium should be applied with station-keeping techniques to prevent escape. A powerful tool is the Hamiltonian-Structure-Preserving control which is constructed as

Tc = T δ r + K δ r˙



(9a)





T = − σ 2 G1 u+ uT+ + G2 u− uT− − γ 2 G3 uuT + u˜u˜T K =ωJ



(9b) (9c)

where ± σ are the hyperbolic eigenvalues associated with the unstable/stable manifolds u ± , ± γ i are center eigenvalues associated with the center manifolds u and u˜, G1 , G2 and G3 are, respectively, referred to as the gains of unstable, stable and center manifolds, J is the symplectic matrix, ω is used to quantify the Coriolis acceleration, referred as the factor quantifying Coriolis acceleration, δ r is the position feedback, and δ r˙ is the velocity feedback. The fact that T is a symmetry matrix and K is a skew symmetry matrix guarantees that the linear feedback controller preserves the Hamiltonian structure [15]. According to Xu and Xu [19], a 2 degree-of-freedom Hamiltonian system can be

286

X. Pan et al. / Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

stabilized using only the invariant (stable, unstable, and center) manifolds and the position feedback, i.e., the first term of Tc in Eq. (9a). Such control strategy is implemented by transforming the hyperbolic equilibrium (saddle) to an elliptic one (center) with the poles assignment and then a new type of quasi-periodic orbit referred to as a stable Lissajous orbit is obtained. Compared with the classical HSP controller of Scheeres [15], the controller developed by Xu and Xu [19] makes further improvement on the basis of retaining the strong characteristics of HSP control. To sort out, the controller using only the position feedback is referred later as the strong HSP. In comparison, the second term of Tc with only the variable Coriolis acceleration and velocity feedback is referred later as the weak HSP. It is obvious that the strong HSP control implemented to the system changes the topology from the saddle to the center, and reduces the number of the equilibria to only one. Furthermore, the system’s energy and Hill’s region become different. Any unstable trajectory is stabilized at the cost of forcing all the unstable region to the stable one. However, the weak HSP control makes no difference to the topology of the original system, for it does not change the Uqq , and as a result, the Hill’s region and the two equilibria always remain the same. To stabilize an unstable trajectory, the weak HSP control only needs to vary the factor ω slightly, which can effectively change the direction of some (not all) trajectories and help to avoid the unstable condition (shown as the Fig. 22). Compared to the strong control changing all unstable trajectories, the weak one takes less effort and is a well-leveraged control which stabilizes motion at the cost of changing only a handful of trajectories. Even though on some degree, velocity is not as easy to measure as the position, the velocity feedback is also required in the classical HSP controller, i.e., Eq. (9a), and then it is clear the weak control has no more requirements than the classical one. Besides, the velocity information can be obtained from position information with the filtering methods, and especially when the position data is enough, the velocity data can be solved by differential method. Furthermore, for the navigation scenario using inertial instruments, it is easier to yield the velocity information from measuring acceleration than the position. Extensive researches have been made for the strong HSP [15,19,21], but less attention was paid on the weak HSP control. Building on the dynamics of displaced orbits over a planet, this paper firstly focuses on the investigation of the weak HSP control to probe into the mechanism of stability change and presents specific utilized strategies on how the control can be implemented. Considering applying the weak HSP control to the displaced orbits above a planet, the controlled dynamical equation of the spacecraft’s motion is written as



ρ¨ z¨

⎡ ∂U ⎤

− ρ˙ ∂ρ ⎣ ⎦ = + ωJ . z˙ ∂U − ∂z

(10)

And the controlled variation equation near the equilibria is

δ q¨ − ωJ δ q˙ + Uqq δ q = 0

(11)

Then the problem of the displaced orbit by the weak HSP control is parameterized by the factor ω, and the Eq. (10) is a general dynamic system where the original uncontrolled model can be regarded as a particular case, i.e., ω = 0. So that in the following sections, the uncontrolled scenario can be simply denoted by a zero ω and the controlled scenario is denoted by a nonzero ω. 3. Controlled motions near hyperbolic equilibrium 2.1. Lyapunov orbit According to the Lyapunov Center Theorem [22], there exists a one-parameter family of periodic orbits emanating from the equilibrium point Lu with one pair of imaginary eigenvalues, referred to as the Lyapunov orbit. Based on the differential correction used in searching for the halo orbits in the circular restricted three body problem (CR3BP), an iterative algorithm to generate Lyapunov orbit is proposed as follows [23]. A good estimation of the periodic orbit is required by the initial state of the iteration algorithm, which is denoted as X0 = [ρ0 z0 ρ˙ 0 z˙ 0 ]T with the period of T0 . The slight variations δρ , δ ρ˙ and δ z˙ are governed by the conservation of energy as follows. The energy for the initial state X0 is given by

E0 =

1 2 1 2 ρ˙ + z˙ + U0 2 0 2 0

(12)

U0 =

 1 hz − − κ z 0 , r 0 = ρ0 2 + z 0 2 . 2 r0 2 ρ0

(13)

2

The three variations based on the constant energy can be obtained as follows:

ρ˙ 0 δ ρ˙ 0 +z˙ 0 δ z˙ 0 +

∂ U0 δρ = 0 ∂ ρ0 0

(14)

which can be rewritten as

δ z˙ 0 = C1 δρ0 + C2 δ ρ˙ 0

(15)

X. Pan et al. / Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

where

C1 =

1 z˙ 0



hz

2

ρ0 3

−

ρ0 r0 3

 , C2 = −

ρ˙ 0 z˙ 0

.

287

(16)

The motion is integrated by the iteration algorithm from X0 to the state of z f = z0 denoted as X f = [ρ f z f ρ˙ f z˙ f ]T , with the integration time T near one orbital period. Therefore, the correction δ X0 to the initial state is required to narrow the deviation, as given by

ϕT0 +δT0 (X0 + δ X0 ) = X0 + δ X0

(17)

where ϕt (X) is the general solution of the Lyapunov orbits at any time t starting from the initial state X0 ; the correction T to the initial state is available for the components as δ X0 = [δρ0 0 δ ρ˙ 0 δ z˙ 0 ] ; and T0 + δ T is the updated period of the corrected state of Xf (δ T0 is the correction to the initial period T0 ). We expand Eq. (17) to achieve the analytical first-order expression as

ϕT0 +δT0 (X0 + δ X0 ) − ϕT0 (X0 ) = MT0 δ X0 +

∂ϕ | δT ∂ t X0 0

(18)

∂ ϕT

where MT0 = ∂ X 0 is the state-transition matrix of the controlled nonlinear dynamics at initial state X0 . 0 Rewriting Eq. (18) yields the following matrix equation:

⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ρ0 + δρ0 − ρ f δρ0 ρ˙ ⎢0 ⎥ ⎢0 ⎥ ⎢z˙ ⎥ ⎣ρ˙ + δ ρ˙ − ρ˙ ⎦ = MT0 ⎣δ ρ˙ ⎦ + ⎣ρ¨ ⎦δ T0 0 0 0 f z˙ 0 +C1 δρ0 + C2 δ ρ˙ 0 − z˙ f C1 δρ0 + C2 δ ρ˙ 0 z¨

(19)

where C1 and C2 are the constraints between δρ 0 , δ ρ˙ 0 and δ z˙ 0 to preserve the energy conservation of the Hamiltonian system; then, matrix Z is introduced as

 Z=

MT0 11 − 1 MT0 31 MT0 41

MT0 13 MT0 33 − 1 MT0 43

MT0 14 MT0 34 MT0 44 − 1

ρ˙ ρ¨



z¨

⎡

1 ⎢0 ·⎣ C1 0

0 1 C2 0

⎤

0 0⎥ 0⎦ 1

(20)

where MT0 i j is the element of MT0 . Therefore, the correction δ X0 can be formulated as

⎡ ⎤ ⎡ ⎤ δρ0 ρ0 − ρ f ⎣δ ρ˙ 0 ⎦ = Z −1 · ⎣ρ˙ 0 − ρ˙ f ⎦. δ T0 z˙ 0 − z˙ f

(21)

Since the correction δ X0 obtained from Eq. (21) is still the first-order approximation, several iterations are required to refine the initial state with the step-by-step corrections on a periodic orbit. The numerical execution indicates that 5–9 iterations are appropriate. The convergence of the iterative algorithm depends on the initial value for the first iteration; thus, two steps for a satisfactory approximation of the true value are presented as follows. a) The linear oscillations generated by Eq. (5) are suitable to solve the initial value; thus, δ X˙ = Aδ X is converted into δY˙ = YδY through a similar transformation of matrix X = V Y , where Y is the diagonal matrix composed of A’s eigenvalues, V is a matrix composed of A’s eigenvectors, and A is constructed through using the linearized equations of motion from Eq. (11) as



A=

0

I

−Uqq

ωJ



.

(22)

4×4

b) With the constraint δ z = 0, the undetermined coefficient in the solution of the linear equation δY˙ = YδY can be determined; then, the initial value X0 can be obtained as well. With different factor quantifying Coriolis acceleration ω, the eigenvalues’ number and signs (positive or negative) of A remain the same, but the magnitude changes, which make the periodic solution present a closed trajectory in different size and shape. Under the same energy, the Lyapunov orbits in the case of ω = 0.005 (blue circle) and ω = 0 (black line) are shown in Fig. 3. It is indicated that in the uncontrolled system, the periodic orbit can reach the zero-velocity surface of the Hill’s region due to the zero velocity of two ends in the orbit, so that the Hill’s region is separated into two parts by the periodic orbit. However, there are no zero-velocity points in the periodic orbit of the controlled system, which results in the interspace

288

X. Pan et al. / Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

Fig. 3. Lyapunov orbit emanating from the hyperbolic equilibrium with different ω: E = −0.011397. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

between the periodic orbit and the zero-velocity surface. In order to distinguish the transit orbits and non-transit ones discussed in Section 5, two concepts: interior region and exterior region are defined here based on the Lyapunov orbit. When E ≤ U(Lu ), the region inside the island is referred to as the interior region, and the region inside the mainland is referred to as the exterior region. When E ≥ U(Lu ), the Hill’s region is separated into two parts by the uncontrolled Lyapunov orbit (ω = 0), where the region below the uncontrolled Lyapunov orbit is referred to as the interior region and the region above is the exterior region. Worth of noting that the interior and exterior regions are determined by the system energy rather than the factor ω. In other words, for different factor ω but same energy, the interior and exterior regions remain the same which are separated by the uncontrolled Lyapunov orbit, as shown in Fig. 3. Remark 1. The controlled Lyapunov orbit generated from the equilibrium presents in a circle shape in position space while the uncontrolled one presents in a line shape. The Lyapunov orbit emanating from the equilibrium configuration is periodic. Without control, i.e., when ω = 0, the integration time transformation t → −t, i.e. positive integration time → negative integration time, does not change the dynamic configuration in Eq. (10), which means the trajectories integrated backwards and forwards from the same initial point with null velocity being identical (the existence of null velocity points can be found in Ref. [6]). So integrating Lyapunov orbit for more than one trajectory period, the trajectory presents in a line shape. However, for the controlled system, the dynamical Eq. (10) is no longer symmetrical about t due to ω = 0, so that the Lyapunov orbit integrated for one period presents in the shape of circle rather than line. 2.2. Invariant manifolds According to the Invariant Manifold Theorem of Hyperbolic Trajectory [24], the stable and unstable manifolds are asymptotically close to or far from the hyperbolic periodic orbit. The stable or unstable manifolds can be achieved numerically from an appropriate initial state by the backward or forward integration in time, listed as follows. For a point XL on the Lyapunov orbit, the unstable manifolds Ws + that asymptotically approach XL in first order are determined by the unstable eigenvector Qs + produced from the state-transition matrix on this point. Through numerically integrating backwards in time, the manifolds can be identified. The initial state for integration is calculated as follows:

Xs+ = XL + d · Q s+

(23)

where d is the estimating distance between the invariant manifolds and the Lyapunov orbit in the direction of Qs + . The magnitude of d should be sufficiently small to avoid violating the linear estimate, but not so small that magnitude damps the integration due to the asymptotic nature. Based on the fact that the invariant manifolds of the equilibrium Lu in the uncontrolled system are homoclinic orbits [6], the relationship between the contact ratio of invariant manifolds and the parameter d in such case is demonstrated to determine how the manifolds contact with the equilibrium Lu . The contact ratio is characterized by the position difference of stable and unstable manifolds on the Poincaré section, i.e., ρ = ρ s −ρ u , where ρ s and ρ u are the ρ components of stable and unstable manifolds on the Poincaré section, respectively. As shown in

X. Pan et al. / Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

289

Fig. 4. Invariant manifolds by the weak HSP control: (a) relationship between the contact ratio of the invariant manifolds of Lu and the parameter d in the ω = 0 case; (b) the parts of the invariant manifolds close to the Lyapunov orbit. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 4(a), the ρ decreases with decreasing d, but the change is negligible when d < 10−4 . Considering the numerical errors, the infinitesimal ρ of the order of 10−9 is equivalent to the fact that the invariant manifolds of equilibrium Lu coincide completely. Therefore, a suitable value of d is obtained as 10−4 . A similar procedure can be employed to approximate the unstable manifolds. Then, the parts of invariant manifolds by the weak HSP control projected onto position space are shown in Fig. 4(b), where the yellow set denotes unstable manifolds, and the green set denotes stable manifolds. Remark 2. The weak HSP control results in an angle between the tangential lines of stable and unstable manifolds originating from the hyperbolic equilibrium, which can be computed by the angle between the stable eigenvector Ys− and the unstable eigenvector Ys+ of the hyperbolic equilibrium. According to Koon et al. [25] and Conley [26], the nonlinear solution near the equilibrium point can be expressed by periodic and exponential terms in which the exponential expansion of manifolds can be written as follows: +∞ 

x(t ) =

Yms+ emλt +

m=0

+∞ 

Yns− e−nλt

(24)

n=0

where the first term denotes the unstable manifold, the second term denotes the stable manifold, and the first-order coefficients Y0 s+ and Y0 s− are equal to the linear eigenvectors Ys+ and Ys− . Thus, the angle ϕ between the tangential lines of the hyperbolic equilibrium’s stable and unstable manifolds can be determined as follows:





∞ s− −nλt ( +∞ Y s+ emλt ) · ( +n=0 Y e )  +∞ ns−  ) = arccos(Y s+ · Y s− ) ϕ = arccos(lim +m∞=0 m s + mλt  ·  −nλt  t→0  Y e Y e m=0 m n=0 n

(25)

Based on the limitation equation above, ϕ can be considered as the angle between the first order or linear eigenvectors Ys+ and Ys− . Referring to Eq. (22), the matrix A can be expanded as

⎡

0 ⎢0 A =⎣ a b

0 0 b c

⎤

1 0 0 −ω

0 1⎥

ω⎦

(26)

0

where a is the first row and first column element of the matrix −Uqq , b is the first row and second column element of the matrix -Uqq , and c is the second row and second column element of the matrix −Uqq . Considering the structural feature of matrix A, there are a pair of conjugate complex eigenvalues and a pair of real eigenvalues. Since Ys+ and Ys − are the eigenvectors of the real eigenvalues, we compute only the real eigenvalues ±λ as a function of ω.



λ=





c + a − ω2 +



2   ω2 − c − a − 4 ac − b2 2

The eigenvectors Ys+ and Ys − of the linearized equations at point Lu can also be expressed with ω as

(27)

290

X. Pan et al. / Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

Fig. 5. Angle between the tangential lines of stable and unstable manifolds originating from the hyperbolic equilibrium: (a) the tangential lines of invariant manifolds in position space with different ω; (b) the variation of angle between invariant manifolds with different ω. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Y s+ =

1 −λ

  1   −λ

1 −λ



1 −λ

−ωλ + b λ2 − c





−ωλ + b λ2 − c

1

 1

T ωλ + b λ2 − c

 ,Y s− = ωλ + b  λ2 − c 



1

λ   1   λ

 ωλ + b 1 λ λ2 − c   1 ωλ + b 1 λ λ2 − c 1



−ωλ + b λ2 − c

T

−ωλ + b λ2 − c

   

(28)

Thus, the angle between the eigenvectors Ys+ and Ys − can be computed as a function of ω by ϕ = arccos(Y s+ · Y s− ). The relationship between the factor ω and the angle ϕ = arccos(Y s+ · Y s− ) is investigated for ω from −0.1 to 0.1, as shown in Fig. 5(b). The angle increases with increasing ω, and the angle only depends on the magnitude of ω rather than its sign. In addition, the angle between the stable eigenvector Ys− and the unstable eigenvector Ys+ will generally flatten out when ω > 0.01. The invariant manifolds of the Lyapunov orbit are coupled to the equilibrium point; thus, the stable and unstable manifolds of the Lyapunov orbit go in separate directions initially under the control of weak HSP as well. And the angle between the invariant manifolds of the Lyapunov orbit can be approximately measured by the angle between the center line or outer boundary of the manifolds. It is expected that the angle between the invariant manifolds of the Lyapunov orbit shares the similar relationship with the factor ω, where the larger ω is, the larger the angle is. 4. Controlled motions near elliptic equilibrium 4.1. Bounded trajectories near elliptic equilibrium When the energy of system is lower than U(Lu ), that is, the island and the mainland are disconnected; all trajectories near Ls will be bounded inside the island. To analyze the influence of the weak HSP control on the motions near elliptic equilibrium Ls , we employ the Poincaré section to demonstrate the family of bounded orbits. Due to the conservation of energy E, the motions can be constrained to three dimensions; then they can be further simplified into a two-dimensional Poincaré section by introducing the transversal plane z = zs , where zs is the z coordinate of Ls . In other words, the Poincaré mapping is plotted by all points (ρ , ρ˙ ) when the flow crosses the transversal plane z = zs from bottom to top. It is found from the Poincaré mappings in Fig. 6 that when E < U(Lu ), there are only two fixed points and many closed curves surrounding the fixed points regardless of the factor ω. The Lyapunov orbits of Ls in two directions are referred to as P.O.I (green line) and P.O.II (black line) in position space. Because of the periodic characteristics, the orbit crossing the transversal plane z = zs from bottom to top only plots one stable fixed point in the Poincaré section, while the other bounded orbits without the periodic characteristics can form a closed curve. In the uncontrolled case described in Ref. [6], the two orbits P.O.I and P.O.II are regarded as the basic periodic modes for the motions of bounded trajectories, and any trajectory will be coupled to both of the modes. The same conclusion can also be achieved for the weak controlled case, except for the difference that both the basic periodic modes in position space are closed trajectory shapes rather than curves as shown in Fig. 6(b), (d) and (f). In Fig. 6, the upper fixed point in (ρ , ρ˙ ) space

X. Pan et al. / Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

291

Fig. 6. Bounded trajectories near Ls when E = 0.9U(Lu ): (a) Poincaré section ω = 0, (b) basic periodic modes ω = 0, (c) Poincaré section ω = 0.001, (d) basic periodic modes ω = 0.001, (e) Poincaré section ω = 0.005, (f) basic periodic modes ω = 0.005. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

292

X. Pan et al. / Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

Fig. 7. Basic modes of trajectories near Ls when E = U(Lu ) and the invariant manifolds of Lu : (a) ω = 0.005, (b) ω = 0.

corresponds to the green Lyapunov orbit P.O.I in (ρ , z) space; the lower fixed point in (ρ , ρ˙ ) space corresponds to the black Lyapunov orbit P.O.II in (ρ , z) space. 4.2. Homoclinic orbits to hyperbolic equilibrium and Lyapunov orbit Increase the energy E gradually until E = U(Lu ); then, the hyperbolic equilibrium becomes the unique point connecting the island and the mainland. Fig. 7(a) demonstrates the basic modes with E = U(Lu ) and the invariant manifolds of Lu when ω = 0.005; it is found that the basic mode P.O.I totally overlaps with the invariant manifolds of the hyperbolic equilibrium. In other words, the stable and unstable manifolds are demonstrated to coincide in position space, and the basic mode P.O.I can evolve into the invariant manifolds staying in the interior region when the energy E increases to U(Lu ). In particular, the case of ω = 0 demonstrated in Fig. 7(b) shares the same conclusion that was proven in Ref. [6], which is concluded in the following remark. Remark 3. Regardless of ω, the stable and unstable manifolds of hyperbolic equilibrium Lu by the weak HSP control coincide in position space but with opposite velocities. Due to the non-integrable feature of the controlled dynamical equation, the series solution of differential equations developed in Ref. [27] is employed in this paper to search the homoclinic orbits of hyperbolic equilibrium Lu in the following system:

X˙ = f (X ) = AX + u

(29)

where X is the variable including the position and its velocity, i.e., X = [ρ zρ˙ z˙ ]T , f is the nonlinear dynamics with the first and second components equal to ρ˙ and z˙ as well as the third and fourth components equal to ρ¨ and z¨, u is the nonlinear part extracted from f, and A is the Jacobian matrix defined in Eq. (22). We transform the initial point of the Eq. (10) to (ρ u , zu , 0, 0); then, the elements of X are assumed to have the corresponding series solution as follows:

ρ (t ) =

⎧ +∞  ⎪ ⎨ ak ekτ1 t t ≥ 0 k=1

+ ∞  ⎪ ⎩ a k ekτ2 t t < 0 k=1

, z(t ) =

⎧ +∞  ⎪ ⎨ bk ekτ1 t t ≥ 0 k=1

+ ∞  ⎪ ⎩ b k ekτ2 t t < 0 k=1

, ρ˙ (t ) =

⎧ +∞  ⎪ ⎨ ck ekτ1 t t ≥ 0 k=1

+ ∞  ⎪ ⎩ c k ekτ2 t t < 0 k=1

, z˙ (t ) =

⎧ +∞  ⎪ ⎨ dk ekτ1 t t ≥ 0 k=1

+ ∞  ⎪ ⎩ d k ekτ2 t t < 0

(30)

k=1

where ak , bk , ck , dk , a k , b k , c k , and d k are undetermined coefficients, and

τ 1 and τ 2 are the negative and positive eigenvalues of the Jacobian matrix A, respectively. Firstly, considering the positive-time orbit passing (ρ u , zu , 0, 0), the substitution of Eq. (30) into Eq. (29) yields

⎡ ⎤ a1

⎢b ⎥ (τ1 I − A )⎣ 1 ⎦= 0. c1 d1

(31)

X. Pan et al. / Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

293

Fig. 8. The relationship between the order and its coefficients: (a) the coefficients of ak , bk , ck , and dk ; (b) the coefficients of ak ’ , bk ’ , ck ’ , and dk ’ .

Because det(τ1 I − A ) = 0, Eq. (31) has the non-zero solutions [a1 b1 c1 d1 ]T , which can be expressed by the parameter ξ . The same approach can be employed to yield [ak bk ck dk ]T , (k > 1), as

⎡k−1 ⎤ ⎡ ⎤ j=1 a j ak ρ˙ k−1 −1 ⎢kj=1 bj⎥ ⎢bk ⎥ ⎢z˙ k−1 ⎥ ( k τ1 I − A ) ⎣ ⎦ = ⎣ − A⎢k−1 ⎥. ⎦ ⎣ j=1 c j ⎦ ck ρ¨ k−1 k−1 dk z¨k−1 d ⎡ ⎤

j=1

(32)

j

It is clear that det(kτ1 I − A ) = 0, (k > 1 ), i.e., [ak bk ck dk ]T , (k > 1), can be determined solely by Eq. (32). Thus, the positivetime orbit is established by Eqs. (30)–(32). Similarly, the negative-time orbit passing (ρ u , zu , 0, 0) can also be solved with the method above. [a k b k c k d k ]T , (k > 1) can be solely determined by [a 1 b 1 c 1 d 1 ]T , which is expressed by the parameter η. Additionally, ξ and η are solved by + ∞ + ∞   ak = a k . k=1

k=1

The coefficients in the series solution are illustrated in Fig. 8 with respect to their orders, which indicates that all eight coefficients will converge to zero; the homoclinic orbit of the hyperbolic equilibrium Lu is demonstrated in Fig. 9. The discussion above only focuses on the invariant ω in the spacecraft’s dynamics; however, the factor qualifying Coriolis acceleration ω can also be time-varying. A simple example is made here to show the behavior of the trajectory under timevarying ω: the stable and unstable manifold of hyperbolic equilibrium Lu are integrated with the ω varying as

ω (t ) = ω0 cos ( · t ).

(33)

where ω0 and ϖare the magnitude and frequency of ω respectively, which are set as 0.01 and 0.0066 in this section. The manifolds of Lu under such cosine-varying ω are demonstrated in Fig. 10, and it is indicated that even though the factor ω is varying with time, the stable and unstable manifolds are still the homoclinic orbits to Lu . Then, homoclinic orbits to the Lyapunov orbit are researched with the numerical method developed by Llibre et al. [28] and Koon et al. [29]. For the hyperbolic characteristics of Lyapunov orbit are inheriting from the hyperbolic equilibrium Lu , the existence of the homoclinic orbits to Lu and homoclinic orbits to Lyapunov orbit can serve as references for each other. Here, the Poincaré section is redefined as the transversal section zp = 14, which is relatively far from the Lyapunov orbit, with the trajectory transiting from top to bottom. The invariant manifolds of the Lyapunov orbit intersect the section zp = 14 in a curve diffeomorphic to a circle. Fig. 11 shows the sections of the stable manifolds (solid line) and the unstable manifolds (dash line) of the Lyapunov orbit at zp = 14 with ω = 0.007 and different energies. According to the numerical results in Fig. 11, even though the energy is variable, the Poincaré sections of the stable and unstable manifolds coincide with each other, which demonstrates the invariant manifolds of the Lyapunov orbit are homoclinic orbits in such case. It is also noted that the size of the closed trajectory on the Poincaré section decreases with decreasing energy and particularly the trajectory degenerates into one point plotted in the Poincaré section when the energy decreases to U(Lu ). This means that the equilibrium case can be viewed as an extension of the Lyapunov orbit with a reduced energy. Therefore, Remark 4 is proposed as follows.

X. Pan et al. / Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

z

294

Fig. 9. Homoclinic orbit to the hyperbolic equilibrium Lu by the series solution: ω = 0.005.

Fig. 10. Homoclinic orbits to hyperbolic equilibrium under time-varying ω: ω 0 = 0.01 and ϖ = 0.0066.

Remark 4. The set formed by stable and unstable manifolds of Lyapunov orbit in the interior region are essentially the homoclinic orbits to Lyapunov orbit. To illustrate the stable manifolds (indicated by green set) and the unstable manifolds (yellow set) of the Lyapunov orbit (black line) connecting homoclinically with each other, four cases of different ω are shown in Fig. 12.

5. Transit and non-transit trajectories by weak HSP control Sections above analyze the influence of the weak HSP control on the motions near two equilibria, which demonstrates most characteristics of the original system can be preserved by such control. The main focus of this section is to investigate the effect of the weak control on the orbital stability.

X. Pan et al. / Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

295

Fig. 11. Poincaré section to search homoclinic orbits to the Lyapunov orbit under different energies: E1 = −0.0119308, E2 = −0.01193716, E3 = −0.0119383, U(Lu ) = −0.01193867.

5.1. Orbit classification and characteristics Firstly, the motions governed by the controlled displaced orbital dynamics can still be classified into the following four categories based on the different characteristics of the orbits, as shown in Fig. 13:

(a) Lyapunov orbit (black line); (b) stable/unstable manifolds (dash lines), which asymptotically tend to the Lyapunov orbit; (c) non-transit orbits (green lines), which stay in the interior or exterior region with infinite time; (d) transit orbits (red lines), which have only finite time to stay in the interior region and will then spread to infinity in the exterior region. When ω = 0, the periodic orbit separating the non-transit orbits inside the interior region from the non-transit orbits inside the exterior region can be considered as the critical boundary, beyond which the trajectories inside the interior region will transit to the exterior region and vice versa. Thus, in this case, the transit orbit must go through the Lyapunov orbit. However, different from ω = 0 which makes the periodic orbit not only be projected to a curve but also reach the zerovelocity surface due to the zero velocity of two ends in the curve, the nonzero ω does not separate the Hill’s region into two parts. So the transit trajectories in this case can go through the Lyapunov orbit (shown as the left red line in Fig. 13) or the interspace between Lyapunov orbit and zero-velocity surface (shown as the right red line in Fig. 13). Stable and unstable manifolds tend to the Lyapunov orbit asymptotically. Without control, the invariant manifolds of the periodic orbit act as the boundary of all the transit trajectories, which can be viewed as the critical orbits between transit and non-transit orbits. In contrast, in the controlled system, transit trajectories are always beyond the invariant manifolds in the neck region due to the interspace between the Lyapunov orbit and the zero-velocity surface. To classify the transit and non-transit orbits, the criterion based on the accessible positions by these orbits of ω = 0 developed in Ref. [6] is extended to the controlled case of ω = 0. Except for the invariant manifolds, the other regions in  are separated into three areas termed I, II and III, and they are defined based on the accessible regions by transit and non-transit orbits, where Area I denotes the region reached by all the transit orbits, which means the points in Area I must stay on transit orbits, and all the transit orbits will pass Area I; Area II and Area III denote the regions reached by all the non-transit orbits in the interior and exterior region, respectively, which means the points in Area II or Area III must stay on non-transit orbits, and all the non-transit orbits in the interior or exterior region will pass Area II or Area III, respectively. Since system presented by Eq. (10) is considered on an energy surface , Areas I, II and III are all three-dimensional. Thus, to completely investigate these three-dimensional areas, two scenarios, (ρ , z) space and (ρ , ρ˙ ) space, are separated to discuss the existence of Areas I, II and III as well as the transition and non-transition behavior.

296

X. Pan et al. / Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

Fig. 12. Homoclinic orbits to the Lyapunov orbit in different cases: (a) ω = 0, E = −0.011937; (b) ω = 0.005, E = −0.011934; (c) ω = 0.01, E = −0.011924; (d) ω = 0.05, E = −0.011908. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

5.2. Classification of the transit and non-transit orbits by accessible positions in the hill’s region of (ρ , z) space In (ρ , z) space, the three areas are projected onto the Hill’s region to distinguish the transit orbits from the non-transit orbits based on the (ρ , z) position. In the uncontrolled system [6], the invariant manifolds can reach the zero-velocity surface. In other words, the accessible Area I by the transit orbits and Areas II and III by the non-transit orbits are separated from each other by the intersections of the manifolds and the zero-velocity surface, which can be found by searching the zero-velocity points in the invariant manifolds. However, no intersection occurs in most controlled cases. Therefore, not all the transit orbits are restricted inside the invariant manifolds, which makes the three accessible areas in the ω = 0 case different from those in the ω = 0 case. The variations of Areas I, II and III with increasing ω are addressed as follows, where the Areas I, II and III labeled in Fig. 14 denote their projection on the Hill’s region. Firstly, different from the Area I of ω = 0 appearing only at the bottom of the byland, the Area I of ω = 0 can also appear in the interspace around the neck region, which is bounded by the zero-velocity surface and invariant manifolds. The bottom and neck-region parts of Area I (denoted by lower Area I and upper Area I, respectively) are separated until ω increases to the critical factor ωc , which is defined as the lowest ω at which lower Area I and upper Area I are connected and there is no intersection between the invariant manifolds and the zero-velocity surface in the interior region. The parameter ω c is determined by the energy E as shown in Fig. 15. It is indicated that when E-U(Lu ) > 5 × 10−6 , i.e., E > −0.0119337, ωc appears and that a larger E implies a smaller ωc . Additionally, ωc remains largely the same until it decreases to 0.002. Therefore, the two parts are connected with each other when ω is larger than its critical factor, as shown in Fig. 14. Secondly, different from the Area II of ω = 0 appearing only on both sides of the manifolds (denoted by side Area II), the Area II corresponding to ω = 0 can appear in the middle of donut-like manifolds (denoted by middle Area II) in the interior region as well. The side Area II is bounded by the zero-velocity surface and invariant manifolds, and it is separated from Area I by the intersections of the zero-velocity surface and the manifolds, as shown in Fig. 14(a) and (b). As ω increases,

X. Pan et al. / Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

297

Fig. 13. Four categories of orbits in the neck region of the position space: white region represents the Hill’s region. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

the side Area II becomes successively smaller until it disappears, while the middle Area II grows larger with the increasing angle between the manifolds. Similar conclusions can be achieved for the Area III in the exterior region. Thirdly, due to the fact listed in Fig. 5(b) that the variation range of angle between the stable and unstable manifolds is less than 1° when ω ≥ 0.01, the geometry of the invariant manifolds in all the cases of ω larger than 0.01 is similar to that of ω = 0.01. Thus, the topologies of Areas I, II and III in the ω > 0.01 case are the same as those at 0.01 as shown in Fig. 14(d), where the lower and upper Area I are connected, and the side Area II disappears with only the middle Area II remaining, and Area III merely appears in the middle of the Y-like manifolds in the exterior region. Based on the accessible positions discussed above, the points in Areas I, II, III are selected to check whether all of the transit and non-transit orbits follow the candidate rules. With the constraint of the given energy, the points (ρ , z) can derive out a series of four-dimensional state (ρ , z, ρ˙ , z˙ ) to integrate the orbit forwards and backwards. And the transitions of the trajectories are checked by examining the z component of trajectories larger than that of equilibrium point Lu or not. Based on the ergodic simulations of the factor ω, we present the following remark concerning the classification of the transit and non-transit orbits in (ρ , z) space as follows. Remark 5. Under the weak HSP control, Areas I, II and III presented in the (ρ , z) space are separated by the intersections of the manifolds and the zero-velocity surface. The interspaces around the neck region and at the bottom of the byland are Area I; the other interspaces in the interior (exterior) region, which are on the sides of the manifolds or in the middle of manifolds, are Area II (Area III). In (ρ , z) space, the existence of points whose locations are in the projection of Area I (Area II or Area III) on the Hill’s region is the necessary and sufficient condition for the transition (non-transition). When spreading to the exterior region, the transit trajectory will be finally restricted within the invariant manifolds in the exterior region. 5.3. Classification of the transit and non-transit orbits by transition region in the Poincaré region of (ρ , ρ˙ ) space In the previous section, we proposed a method in Remark 5 to determine whether the trajectory can transit by examining only whether there are points of the trajectory located in the projection of Areas I or II onto the Hill’s region. Since not all the transit orbits are restricted inside the invariant manifolds in the controlled case, the colored Poincaré section (i.e., Fig. 16) is demonstrated to investigate the relationship between the manifolds and the transition. Due to the conservative energy E, a Poincaré mapping is employed to reduce the four components (ρ , z, ρ˙ , z˙ ) of any state to the two components of (ρ , ρ˙ ). Any trajectory passing through the transversal section can be mapped into the pair (ρ , ρ˙ ) on this section because the other component z˙ can be solved from the conservation of energy E as 2

E=

1 1 2 1 2 h ρ˙ + z˙ + z 2 − − κ z p = constant 2 2 r 2ρ

(34)

298

X. Pan et al. / Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

Fig. 14. Accessible positions by transit and non-transit orbits with different values of ω: (a) ω = 0, E = −0.011937; (b) ω = 0.002, E = −0.0119366; (c) ω = 0.006, E = −0.011932; (d) ω = 0.01, E = −0.011924.

Since the term z˙ 2 is always larger than or equal to zero, all the points (ρ , ρ˙ ) plotted on the Poincaré section are governed by the inequality constraint as

z˙ 2 = 2E − ρ˙ 2 −

hz

2

ρ2

+

2 + 2κ z p ≥ 0 r

(35)

Thus, the Poincaré region is defined as the feasible region in (ρ , ρ˙ ) space satisfying this inequality constraint Eq. (35), and the critical circle is defined as the boundary of the Poincaré region, which can be formulized by the equality constraint 2

2E − ρ˙ 2 − hρz2 +

2 r

+ 2κ z p =0, shown as the black circle in Fig. 16.

X. Pan et al. / Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

299

Fig. 15. Critical factor ωc under different energy E: U(Lu ) = −0.01193867.

Fig. 16. Transition region of different ω: (a) ω = 0, E = −0.011937; (b) ω = 0.0065, E = 0.01193366; (c) ω = 0.008, E = −0.01193366; (d) ω = 0.01, E = −0.011924; the transversal section zp = zs . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

300

X. Pan et al. / Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

Fig. 17. Two types of transition regions. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Any point bounded by the critical circle in the transversal section z = z p represents two states associated with positive and negative z˙ p . We integrate from both of the two states to yield the whole trajectories and then determine whether these trajectories can transit to the exterior region by examining the z component of trajectories larger than that of the equilibrium point Lu . Subsequently, in the (ρ , ρ˙ ) space, we mark the points generating transit orbits with only positive values of z˙ p by red dots, mark the ones with only negative values of z˙ p by green dots, and mark the points that can generate transit orbits with both positive and negative values of z˙ p by yellow dots. Thus, all the painted regions in Fig. 16 are termed Transition Regions, where the red transition region is referred to as the positive transition, the green transition region is referred to as the negative transition, and the yellow transition region is referred to as the positive and negative transition. In addition, the left (or the lower in ω = 0 case) blue circle in Fig. 16 denotes the intersection of the stable and unstable manifolds crossing section z = z p from top to bottom; the right (or the upper in ω = 0 case) blue circle denotes the intersection of the stable and unstable manifolds crossing section z = z p from bottom to top. Both of the intersections are referred to as Mapped Manifolds. The location of the manifolds mapped on the Poincaré section is investigated under the traversal of factor ω. It is indicated that there exist only the two scenarios for the transition regions of all ω: one is the mapped manifolds locating inside and without any intersection with critical circle (Fig. 17(c) ω = 0), and the other is the mapped manifolds internally tangent to the critical circle (Fig. 17(d) ω = 0.01). The geometrical relationship between the accessible positions in (ρ , z) space and the transition region in (ρ , ρ˙ ) space is demonstrated in Fig. 17 along with the manifolds. In Fig. 17, the transversal section is located at zp = zs ; in (ρ , z) space, the black straight line is projected by the Poincaré region on the transversal section z = zp , the blue solid (dashed) lines denote the transit (non-transit) orbits integrated from the point in the negative transition region, the red solid (dashed) lines denote the transit (non-transit) orbits integrated from the point in positive transition region; and in (ρ , ρ˙ ) space, the red and green circles represent the mapped manifolds on the Poincaré section, and the black circle represents the critical circle. In the first scenario shown in Fig. 17(a) and (c), the mapped manifolds on the Poincaré section are located inside and have no intersection with the critical circle in (ρ , ρ˙ ) space. Clearly, the points in segments ab and e f with any admissible velocity ρ˙ are located outside the mapped manifolds in (ρ , ρ˙ ) space, and these points generate non-transit orbits regardless of whether z˙ p is negative or positive. For the segments bd and ce in the projection of the invariant manifolds, it is found that under only the condition of (ρ , ρ˙ ) located inside the lower branch of the mapped manifolds and a negative z˙ p , the point in segment bd can generate a transit orbit; similarly, only under the condition of (ρ , ρ˙ ) located inside the upper mapped manifolds and a positive z˙ p , the point in segment ce can generate a transit orbit. If (ρ , ρ˙ ) is located outside the mapped manifolds, the points in segments bd and ce generate the non-transit orbit whether z˙ p is negative or positive. In conclusion, according to the definition, Area II is the set of points that can generate non-transit orbits with both negative and positive z˙ p . In (ρ , z) space, Area II is partially represented by segments ab and e f because some points in Area II are in the position

X. Pan et al. / Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

301

Fig. 18. Transition duration of transit orbits remaining in the interior region: (a) ω = 0; (b) ω = 0.001. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 19. Bounded trajectory transfer from the Poincaré section: transversal plane zp = zs , z˙ p > 0. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

of the projection of invariant manifolds, while in (ρ , ρ˙ ) space, Area II of this transversal section is completely represented by the mapping region located outside the mapped manifolds. Essentially, the three-dimensional Area II is interconnected. Thus, the points in segments ab, and e f of position space and the points located outside the mapped manifolds of (ρ , ρ˙ ) space generate the non-transit orbits with both negative and positive z˙ p . In addition, an interesting conclusion can be obtained from Fig. 17(a) that any point located inside the lower branch with a negative z˙ p will generate a transit orbit arriving in Area I, shown as a blue solid line in position space, while the same point with a positive z˙ p will generate a non-transit orbit arriving in Area II, shown as a blue dash line in position space. The opposite conclusion can be achieved for the point inside the upper branch, shown as red lines in position space. In the second scenario shown in Fig. 17(b) and (d), the mapped manifolds are internally tangent with the critical circle, and the two yellow-colored regions are partitioned by the tangency. Similarly, with (ρ , ρ˙ ) located inside the left branch of mapped manifolds and a negative z˙ p , the point in segment bc generates the transit orbits, and with (ρ , ρ˙ ) located inside the right branch and a positive z˙ p , the point in de generates the transit orbits. The line segment cd in (ρ , z) space and the mapping region between the mapped manifolds in (ρ , ρ˙ ) space represent all the non-transit orbits arriving at the threedimensional Area II regardless of whether z˙ p is negative or positive. As the lower and upper Area I are connected as the

302

X. Pan et al. / Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

Fig. 20. Trajectory transfer from torus#1 to torus#3 from the position space. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 21. An example of transition control: E = 0.0119371. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

“gaps” between the manifolds and the zero-velocity surface shown in Fig. 17(b), the segments ab and e f in (ρ , z) space and the mapping region on the sides of the mapped manifolds in (ρ , ρ˙ ) space represent all the transit orbits arriving at the three-dimensional Area I regardless of whether z˙ p is negative or positive. In addition, since the scenario of transition regions is dependent on the transversal section, a conclusion regarding ω c can be drawn that the value of ω passing its critical values is the sufficient but unnecessary condition to create the positive and negative transition region. According to the description mentioned above, a new remark from the aspect of (ρ , ρ˙ ) space can be summarized as follows:

X. Pan et al. / Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

303

Remark 6. In (ρ , ρ˙ ) space of the controlled system, the initial point located inside the left (or the lower in the uncontrolled case) mapped manifolds with negative z˙ p or the right (or the upper in the uncontrolled case) mapped manifolds with positive z˙ p can generate a transit orbit; if there are no intersections between the critical circle and mapped manifolds, the parts outside the mapped manifolds represent the Area II in (ρ , ρ˙ ) space so that the initial point located internally can generate a non-transit orbit for either positive or negative z˙ p ; if the mapped manifolds are internally tangent to the critical circle, the new regions partitioned by tangency represent the Area I in (ρ , ρ˙ ) space so that an initial point located within it can generate a transit orbit for either positive or negative z˙ p , but the other regions represent the Area II in (ρ , ρ˙ ) space, within which an initial point can generate a non-transit orbit for either positive or negative z˙ p . Similar to Remarks 5 and 6 is useful to determine whether an orbit can transit or not merely based on the location of its initial point in the Poincaré region without extensive numerical computations. In summary, to determine whether the trajectory integrated from an initial point can transit between the interior and exterior regions, we present the first (ρ , z)space method based on the accessible positions by transit and non-transit orbits demonstrated in Remark 5 and Fig. 14 and the second (ρ , ρ˙ )-space method based on the locations of the manifolds mapped onto the Poincaré region as demonstrated in Remark 6 and Fig. 16. Then, an algorithm to determine whether a single point (ρ 0 , z0 ) can generate a transit orbit or not is listed as follows: if the point is located in one of the projections of Areas I, II and III in (ρ , z) space, we can employ the first method to determine whether the orbit integrated from the point transits directly; if it is located in other areas, we can define a Poincaré section at z = z0 and then use the second method to determine the transition. Furthermore, we consider examining whether the two methods are coherent by introducing Fig. 17 to explain how the transit or non-transit orbit evolves from the initial points in the Poincaré region and revealing the correspondence relationship between a threedimensional area projection on (ρ , z) space and (ρ , ρ˙ ) space. It is concluded that under the action of weak HSP control, the unstable trajectories can be stabilized through changing the direction and shape to avoid the transition region. 5.4. Transition duration of transit trajectories To investigate the properties of the transition region further, the transition duration is defined to measure how long the transit orbit stays in the interior region before it leaves. Taking the positive transition region in Fig. 16 as an example, the transition duration is computed for the total points in such region. As shown in Fig. 18, the area in dark blue has the shortest transition duration, while the transition duration of area in red is much longer. It is indicated that the shortest transition duration always occurs in the interior area of the transition region, which acts as a suitable reference to define the center of the transition region. More importantly, the point with the shortest transition time can be regarded as the optimal choice to design a transit orbit. 6. Application of the weak HSP into orbital transfers 6.1. Transfer within the KAM tori According to Remark 6 proposed above, the transition region changes with factor ω and so does the transition of one orbit integrated from the initial point with four dimensions. Although the transition and non-transition regions of different factor ω in Fig. 16 are presented under the different energies, the conclusion above still holds under the same energy. Thus, with the given energy constraint, we present the use of weak HSP control as the valuable tool to change the motion between different bounded orbits in Section. 6.1, and to stabilize the motions in the neck region in Section. 6.2. Since the control factor ω makes no change to the system energy, it will be demonstrated that the transfer between two KAM tori spanned by quasi-periodic trajectories characterized by ω can be made by Poincaré mapping method, and then large families of new artificial orbits are generated. The specific algorithm to the transfer of the quasi-periodic trajectories is addressed as follows: for any two tori, if there exists intersection point, orbital transfer between the two tori can be achieved at such point by changing ω; if there exists no intersection point, the intermediate patches are needed to connect these two separated tori so that this scenario can be transformed into the multi-segments intersecting patching. Under the same energy, any tori remained in the case of ω = 0 will never have intersection point. The typical example of the transfer between the separated torus#1 and torus#3 is shown in Fig. 19 from the Poincaré section, where torus#1 and torus#3 are the tori spanned by quasi-periodic trajectories of ω = 0, denoted by red and green curves respectively, and the intermediate torus#2 is spanned by quasi-periodic trajectories of ω = 0.009, denoted by the yellow curve. It is clear from Fig. 19, the transfer from the torus#1 to the final #3 can be achieved by the first transfer from #1 to #2 through altering ω from 0 to 0.009 at the intersection point P1 (or P4 ), and then the second transfer from #2 to #3 by recovering ω from 0.009 to 0 at the intersection point P2 (or P3 ). Worth of noting that the four intersection points P1 − 4 form a circulation loop, which can not only realize the orbital transfer from the torus#1 to #3, but also return back to torus#1 itself. The intermediate patches at P1 , and P2 give a candidate path transferring from the initial torus #1 to the final #3, which means the patching algorithm developed here is a powerful tool to obtain the transfer between the two separated tori without any intersection points. The specific three tori mentioned above are presented in the position space, shown as Fig. 20. It indicates that after twice transfers at P1 , P2 points, the quasi-periodic trajectory near the P.O.I basic periodic trajectory can be transformed to the one

304

X. Pan et al. / Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

Fig. 22. Transition control of point A1. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

near the P.O.II basic periodic trajectory, and the patched trajectory can reach more accessible positions occupied by the tori #1 and #3. Thus the application of using the weak HSP control to extend several natural segments to an artificial trajectory is practical. 6.2. Transfer beyond the KAM tori In addition to the transfer within KAM Tori, the factor ω can also be used to achieve the transfer between transit orbits and non-transit orbits, because the transition of one orbit integrated from the initial point may be changed by varying ω as well. Next, an example is provided to explain this point in detail. In Fig. 21, the Poincaré region at zp = 15 is bounded by the critical circle (black circle), and the transition regions in different cases of ω are demonstrated under the same energy, where the blue solid (dash) circles represent the positive (negative) transition region boundary in the ω = 0 case and the green solid (dash) circles represent positive (negative) transition region boundary in the ω = 0.005 case. When ω = 0, the point A1 is in the positive transition region; i.e., the orbit integrated from A1 will transit. However, when ω is changed to 0.005, point A1 is out of its corresponding transition region, which generates a non-transit orbit. The switch from the transition to non-transition of the orbit integrated from point A1 is clearly demonstrated in the position space of Fig. 22, where the blue line represents the KAM torus when ω = 0.005; the black line represents the transit orbit when ω = 0; the red line represents the zero-velocity surface when E = 0.0119371. Similarly, the orbit integrated from A2 can transit when ω = 0.005 but remain as a non-transit orbit when ω = 0. Since the location of the transition region varies with the factor ω, ω can be used as a control parameter to govern the transition or non-transition of the orbit, which is of a certain significance for trajectory design. 7. Conclusion This paper introduces a novel concept of the weak Hamiltonian-Structure-Preserving (HSP) control parameterized by the factor quantifying Coriolis acceleration, and investigates the mechanism of changing the orbital stability. The dynamics of controlled system is building on the displaced orbits above a planet generated by a low-thrust spacecraft, and two invariant equilibria, a hyperbolic one and an elliptic one, are solved from it. Then, the influence of the weak control on the motions near the two equilibria is investigated by means of dynamical system techniques. The interior region is filled with the bounded trajectories, which are coupled by two basic periodic modes. The existence of the homoclinic orbits to the Lyapunov orbit as well as the hyperbolic equilibrium is retained and numerically demonstrated. Controlled trajectories in the neck region are analyzed for their characteristics, classification and transition to prepare for the stabilization of unstable orbits near hyperbolic equilibrium. Finally, the orbital transfer between different bounded orbits as well as transfer between stable and unstable orbits are exemplified to show the function and application of the weak HSP control. The important contributions of this paper are as follows: angle between the tangential lines of stable and unstable manifolds originating from the hyperbolic equilibrium in the controlled system has been expressed as a function of the Coriolis

X. Pan et al. / Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

305

acceleration, so that the direction and shape of trajectory can be changed to keep away from the unstable region; it numerically demonstrates that the set formed by the stable and unstable manifolds of Lyapunov orbit in the interior region are still the homoclinic orbits to Lyapunov orbit, which work as the important boundary and criterion to distinguish the transition and non-transition region; due to the transition region varying with the control factor, the weak HSP control succeeds in changing the transition or non-transition of one orbit based on the two methods in determining whether the controlled trajectory can exhibit transit, which is of significance for trajectory design. Acknowledgments The authors would like to acknowledge financial support from the National Natural Science Foundation of China (11772024 and 11432001) and the Fundamental Research Funds for the Central Universities (YWF-16-BJ-Y-10). Appendix

Table 1 Coefficients of the series solution of the homoclinic orbit to the hyperbolic equilibrium. k

ak

bk

ck

dk

a k

b k

c k

d k

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

−7.7642 1.1055 −0.1072 −0.0999 −0.0743 −0.0523 −0.0356 −0.0232 −0.0141 −0.0076 −0.0029 0.0 0 05 0.0029 0.0046 0.0057 0.0065 0.0069 0.0072 0.0073 0.0073 0.0072 0.0071 0.0069 0.0067 0.0065 0.0063 0.0061 0.0058 0.0056 0.0054 0.0052 0.0050 0.0048 0.0046 0.0044 0.0043 0.0041 0.0040 0.0038 0.0037 0.0035 0.0034 0.0033 0.0032 0.0031 0.0029 0.0028 0.0028 0.0027 0.0026

−19.5120 8.5571 0.8268 0.4372 0.2921 0.2166 0.1699 0.1378 0.1143 0.0963 0.0821 0.0707 0.0613 0.0536 0.0472 0.0417 0.0371 0.0332 0.0298 0.0269 0.0243 0.0221 0.0202 0.0184 0.0169 0.0156 0.0144 0.0133 0.0124 0.0115 0.0107 0.0100 0.0094 0.0088 0.0082 0.0078 0.0073 0.0069 0.0065 0.0062 0.0058 0.0055 0.0053 0.0050 0.0048 0.0046 0.0043 0.0042 0.0040 0.0038

0.0034 −0.0010 0.0 0 01 0.0 0 02 0.0 0 02 0.0 0 01 0.0 0 01 0.0 0 01 0.0 0 01 0.0 0 0 0 0.0 0 0 0 0.0 0 0 0 0.0 0 0 0 0.0 0 0 0 0.0 0 0 0 0.0 0 0 0 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01

0.0084 −0.0074 −0.0011 −0.0 0 08 −0.0 0 06 −0.0 0 06 −0.0 0 05 −0.0 0 05 −0.0 0 04 −0.0 0 04 −0.0 0 04 −0.0 0 04 −0.0 0 03 −0.0 0 03 −0.0 0 03 −0.0 0 03 −0.0 0 03 −0.0 0 03 −0.0 0 02 −0.0 0 02 −0.0 0 02 −0.0 0 02 −0.0 0 02 −0.0 0 02 −0.0 0 02 −0.0 0 02 −0.0 0 02 −0.0 0 02 −0.0 0 02 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01 −0.0 0 01

−16.1696 8.1706 1.2119 0.6153 0.3865 0.2677 0.1961 0.1490 0.1162 0.0924 0.0746 0.0610 0.0504 0.0420 0.0354 0.0300 0.0256 0.0220 0.0190 0.0165 0.0144 0.0126 0.0111 0.0098 0.0087 0.0078 0.0070 0.0062 0.0056 0.0051 0.0046 0.0042 0.0038 0.0034 0.0031 0.0029 0.0026 0.0024 0.0022 0.0021 0.0019 0.0018 0.0016 0.0015 0.0014 0.0013 0.0012 0.0011 0.0011 0.0010

−15.6459 4.1010 0.7614 0.3959 0.2620 0.1944 0.1538 0.1266 0.1069 0.0919 0.0801 0.0705 0.0625 0.0559 0.0502 0.0453 0.0410 0.0373 0.0341 0.0313 0.0287 0.0265 0.0245 0.0227 0.0211 0.0196 0.0183 0.0171 0.0160 0.0151 0.0141 0.0133 0.0126 0.0119 0.0112 0.0106 0.0101 0.0096 0.0091 0.0087 0.0083 0.0079 0.0075 0.0072 0.0069 0.0066 0.0063 0.0061 0.0058 0.0056

−0.0070 0.0071 0.0016 0.0011 0.0 0 08 0.0 0 07 0.0 0 06 0.0 0 05 0.0 0 05 0.0 0 04 0.0 0 04 0.0 0 03 0.0 0 03 0.0 0 03 0.0 0 02 0.0 0 02 0.0 0 02 0.0 0 02 0.0 0 02 0.0 0 01 0.0 0 01 0.0 0 01 0.0 0 01 0.0 0 01 0.0 0 01 0.0 0 01 0.0 0 01 0.0 0 01 0.0 0 01 0.0 0 01 0.0 0 01 0.0 0 01 0.0 0 01 0.0 0 01 0.0 0 0 0 0.0 0 0 0 0.0 0 0 0 0.0 0 0 0 0.0 0 0 0 0.0 0 0 0 0.0 0 0 0 0.0 0 0 0 0.0 0 0 0 0.0 0 0 0 0.0 0 0 0 0.0 0 0 0 0.0 0 0 0 0.0 0 0 0 0.0 0 0 0 0.0 0 0 0

−0.0068 0.0035 0.0010 0.0 0 07 0.0 0 06 0.0 0 05 0.0 0 05 0.0 0 04 0.0 0 04 0.0 0 04 0.0 0 04 0.0 0 04 0.0 0 04 0.0 0 03 0.0 0 03 0.0 0 03 0.0 0 03 0.0 0 03 0.0 0 03 0.0 0 03 0.0 0 03 0.0 0 03 0.0 0 02 0.0 0 02 0.0 0 02 0.0 0 02 0.0 0 02 0.0 0 02 0.0 0 02 0.0 0 02 0.0 0 02 0.0 0 02 0.0 0 02 0.0 0 02 0.0 0 02 0.0 0 02 0.0 0 02 0.0 0 02 0.0 0 02 0.0 0 02 0.0 0 01 0.0 0 01 0.0 0 01 0.0 0 01 0.0 0 01 0.0 0 01 0.0 0 01 0.0 0 01 0.0 0 01 0.0 0 01

τ 1 = −0.0 0 043257 τ 2 = 0.0 0 043257.

306

X. Pan et al. / Commun Nonlinear Sci Numer Simulat 62 (2018) 282–306

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29]

Dusek HM. Optimal station keeping at collinear points. Progr Astronaut Aeronaut 1966;17:37–44. doi:10.2514/6.1969-906. McInnes CR. The existence and stability of families of displaced two-body orbits. Celestial Mech Dyn Astron 1997;67(2):167–80. McInnes CR. Dynamics, stability, and control of displaced non-Keplerian orbits. J Guidance Control Dyn 1998;21(5):799–805. doi:10.2514/2.4309. Dankowicz H. Some special orbits in the two-body problem with radiation pressure. Celestial Mech Dyn Astron 1994;58(4):353–70. doi:10.1007/ BF00692010. Scheeres DJ. Stability of hovering orbits around small bodies. Adv Astron Sci 1999;102:855–75. Xu M, Xu SJ. Nonlinear dynamical analysis for displaced orbits above a planet. Celestial Mech Dyn Astron 2008;102(4):327–53. doi:10.1007/ s10569- 008- 9171- 4. McInnes CR. Displaced non-Keplerian orbits using impulsive thrust. Celestial Mech Dyn Astron 2011;110(3):199–215. doi:10.1007/s10569- 011- 9351- 5. Niccolai L, Quarta AA, Mengali G. Electric sail elliptic displaced orbits with advanced thrust model. Acta Astronaut 2016. doi:10.1016/j.actaastro.2016. 10.036. Heiligers J, Ceriotti M, McInnes CR, et al. Displaced geostationary orbit design using hybrid sail propulsion. J Guidance Control Dyn 2011;34(6):1852–66. doi:10.2514/1.53807. McInnes CR, Simmons JFL. Solar sail Halo orbits I: heliocentric case. J Spacecraft Rockets 1992;29(4):466–71. doi:10.2514/3.25487. Simo J, Mcinnes CR. Solar sail orbits at the Earth–Moon libration points. Commun Nonlinear Sci Numer Simul 2009;14(12):4191–6. Baig S, Mcinnes CR. Light-levitated geostationary cylindrical orbits are feasible. J Guidance Control Dyn 2010;33(3):782–93. doi:10.2514/1.46681. McInnes CR. Passive control of displaced solar sail orbits. J Guidance Control Dyn 1998;21(6):975–82. doi:10.2514/2.4334. Bookless J, McInnes C. Dynamics and control of displaced periodic orbits using solar-sail propulsion. J Guidance Control Dyn 2006;29(3):527–37. doi:10.2514/1.15655. Scheeres DJ, Hsiao F-Y, Vinh NX. Stabilizing motion relative to an unstable orbit: applications to spacecraft formation flight. J Guidance Control Dyn 2003;26(1):62–73. doi:10.2514/2.5015. Bookless J, McInnes C. Control of Lagrange point orbits using solar sail propulsion. Acta Astronaut 2008;62(62):159–76. doi:10.1016/j.actaastro.2006. 12.051. Qian H, Zheng JH, Yu XZ, Gao D. Dynamics and control of displaced orbits for solar sail spacecraft. Chin J Space Sci 2013;33(4):458–64. McInnes AI. Strategies for solar sail mission design in the circular restricted three-body problem [M.S. Thesis]. West Lafayette: Purdue University; 20 0 0. Xu M, Xu SJ. Structure-preserving stabilization for Hamiltonian system and its applications in solar sail. J Guidance Control Dyn 20 09;32(3):997–10 04. doi:10.2514/1.34757. Wang SK, Lo MW, Marsden JE, Ross SD. Dynamical systems, the three-body problem and space mission design. In: International conference on differential equations; 20 0 0. p. 1167–81. Soldini S, Colombo C, Walker SJI. Solar radiation pressure Hamiltonian feedback control for unstable libration-point orbits. J Guidance Control Dyn 2017;40(6):1374–89. doi:10.2514/1.G002090. Meyer KR, Hall GR. Introduction to Hamiltonian mechanics and the N-body problem. In: Applied mathematical sciences, vol. 90. Springer-Verlag; 1992. p. 1–32. Xu M, Xu SJ. J2 invariant relative orbits via differential correction algorithm. Acta Mech Sinica 2007;23(5):585–95. Wiggins S. Introduction to applied nonlinear dynamical systems and chaos, 4. Springer-Verlag; 1990. p. 563. Koon WS, Lo MW, Marsden JE, Ross SD. Dynamical systems, the three-body problem and space mission design. In: International conference on differential equations; 20 0 0. p. 1167–81. doi:10.1142/9789812792617_0222. Conley CC. Low energy transit orbits in the restricted three-body problem. Siam J Appl Math 1968;97(4):732–46. doi:10.1137/0116060. Zhang QC, Tian RL, Wang W. Chaotic properties of mechanically and electrically coupled nonlinear dynamical system. Phys J 2008;57(5):2799–804. Llibre J, Martínez R, Simó C. Tranversality of the invariant manifolds associated to the lyapunov family of periodic orbits near L2, in the restricted three-body problem. J Differ Equ 1985;58(1):104–56. doi:10.1016/0 022-0396(85)90 024-5. Koon WS, Lo MW, Marsden JE, Ross SD. Heteroclinic connections between periodic orbits and resonance transitions in celestial mechanics. Chaos 20 0 0;10(2):427–69. doi:10.1063/1.166509.