The Hamiltonian Structure-Preserving Control and Some Applications to Nonlinear Astrodynamics

Abstract

A systematic research on the structure-preserving controller is investigated in this paper, including its applications to the second-order, first-order, time-periodic, or degenerated astrodynamics, respectively. The general form of the controller is deduced for the typical Hamiltonian system in full feedback and position-only feedback modes, which is successful in changing the hyperbolic equilibrium to an elliptic one. With the poles assigned at any different positions on imaginary axis, the controlled Hamiltonian system is Lyapunov stable. The Floquet multiplier is employed to measure the stability of time-dependent Hamiltonian system, because the equilibrium of periodic system may be unstable even though the equilibrium is always elliptic. One type of periodic orbits is achieved by the resonant conditions of control gains, and another type is making judicious choice in the foundational motions with different frequencies. The control gains are selected from the viewpoint of both the local and global optimizations on fuel cost. This controller is applied to some astrodynamics to achieve some interesting conclusions, including stable lissajous orbits in solar sail’s three-body problem and degenerated two-body problem, quasiperiodic formation flying on a 𝐽₂-perturbed mean circular orbit, and controlled frozen orbits for a spacecraft with a high area-to-mass ratio.

1. Introduction

It is well known that most of the astrodynamical problems could be classified as hyperbolic Hamiltonian systems, for example, circular restricted 3-body problem (CR3BP). There are the hyperbolic equilibria which have stable, unstable, and centre manifolds for hyperbolic Hamiltonian system, just like the collinear libration points 𝐿₁,𝐿₂, and 𝐿₃ in CR3BP [1], of which the unstable manifolds may cause instability in the Lyapunov sense.

Some astronautical missions, requiring increasing the coverage rate to the object (ground station or target spacecraft) or avoiding communication signals being lost in the sun, are preferring bounded motions near hyperbolic equilibria, such as lissajous (or halo) orbit generated in CR3BP [1, 2]. There presented two different approaches to generate bounded relative trajectories, of which one is obtaining the bounded orbits by making judicious choice in initial conditions. Certain zero-steady controllers could be applied to stabilize the nominal trajectories, owing to the hyperbolicity inheriting from the topological type of the equilibrium.

However, for the other case without a natural trajectory to track, some station keeping maneuvers are necessary to implement modification of the topology of the equilibrium. This paper proposes a structure-preserving controller to generate bounded trajectories around the equilibrium, which has many potential applications in astrodynamics. Firstly, a simple structure-preserving controller was constructed by Scheeres et al. [3], in which involved only the instantaneous stable and unstable manifolds, and the manifold’s gains are constrained to be equal. By choosing the proper manifolds’ gains, M. Xu and S. Xu [4] have demonstrated that the poles of the typical second-order Hamiltonian system can be allocated to any position along the imaginary axis and then obtained a new quasiperiodic orbit type referred to as a stable lissajous orbit. Therefore, the controller was powerful enough to work out some famous and difficult problems in sail’s restricted three-body problem that there exist the bounded trajectories near the off-axis solar sail equilibrium [4], and in maintaining the eccentricity in the neighborhood of the unstable fixed point. Besides, the controller implemented the stabilization on a saddle-node bifurcation point for the degenerated case, which makes the modified elliptic equilibrium unique for the controlled system to change the instability of the equilibrium [5]. Further work was implemented on the stabilization of the equilibrium for time-periodic system, which has time varying topological types and no fixed-dimensional unstable/stable/center manifolds [6]. The 2-dimensions temporally independent systems possess fixed-dimensions of unstable/stable/center manifolds and a permanent pair of hyperbolic eigenvalues. Tabuada and Pappas [7] focused on local accessibility preserving abstractions, and provided conditions under which local accessibility properties of the abstracted Hamiltonian system are equivalent to the local accessibility properties of the original Hamiltonian control system.

A systematic investigation on the structure-preserving controller is presented in this paper, including some applications in second-order, first-order, time-periodic, and degenerated Hamiltonian astrodynamics, respectively. The general form of the controller is designed for the typical second-order Hamiltonian system, which can work as full feedback and position-only modes. The controller is successful in changing the hyperbolic equilibrium (saddle) to an elliptic one (center) with the poles on the imaginary axis (marginal stability), and then the controlled Hamiltonian system can achieve the Lyapunov stability by means of the Morse lemma. In contrast to the time-independent Hamiltonian system, the equilibrium of periodic system may be unstable even though the equilibrium is always elliptic during its period. Therefore, it is necessary to check whether the Floquet multipliers lie on the unit circle in the complex plane, or whether the moduli of all the multipliers are equal to 1. Two different approaches to generate the periodic orbits are discussed, one of which is achieved by the resonant conditions of control gains for the resonant periodic orbits, and another is making judicious choice in the foundational motions with different frequencies. Furthermore, some selection techniques on the gains are considered from the point view of both the local and global optimizations.

This controller is applied to some astrodynamics to achieve some interesting conclusions, including stable lissajous orbits in solar sail’s three-body problem and degenerated two-body problem, quasiperiodic formation flying on a 𝐽₂-perturbed mean circular orbit, controlled frozen orbits for a spacecraft with high area-to-mass ratio.

2. Hamiltonian System and Hyperbolic Equilibrium

2.1. Typical Hamiltonian System in Astrodynamics

For the typical Hamiltonian system in astrodynamics, the Hamiltonian function has the following form as [4]: 1𝐻=2𝐩^𝑇𝐩+𝜔𝐩^𝑇𝐉𝐪−𝑉(𝐪,𝑡), (1) where 𝑉 is the general form of pseudopotential function relying only on the coordinate 𝐪 and time; 𝜔 is the characteristic parameter introduced by the Coriolis acceleration, such as the mean angular motion of the synodic frame used in CR3BP, and Local Vertical-Local Horizontal (LVLH) frame used in formation flying.

In general, the variables (𝐪,𝐩) can be expressed by the physical variables (𝐫,𝐫̇) from the Legendre transformation: 𝐪𝐩=𝐫𝐫̇𝐈0−𝜔𝐉𝐈, (2) where 𝐈 is the identical matrix, and 𝐉 is the symplectic matrix. The astrodynamics in physical variables can be expressed as 𝑥̈−2𝜔𝑦̇=𝜕𝑉,̈̇𝜕𝑥𝑦+2𝜔𝑥=𝜕𝑉.𝜕𝑦 (3) For a hyperbolic Hamiltonian system, there exist hyperbolic equilibria that have stable, unstable, and center manifolds [1], which are similar to the collinear libration points 𝐿₁, 𝐿₂, and 𝐿₃ in CR3BP. Hence, the equilibrium of the typical 2-dimensional Hamiltonian system in astrodynamics can be obtained from the following equations: 𝐫=𝐫_𝟎,𝐫̇=𝟎𝜕𝑉𝜕𝐫_𝐫𝟎=𝟎. (4) Denote 𝑉_𝐫𝐫 as the second derivative matrix of the pseudopotential function 𝑉 about the position vector 𝐫=[𝑥,𝑦], and the elements of the matrix have the general expression as 𝑉_𝑚𝑛=𝜕²𝑉/𝜕𝑚𝜕𝑛,(𝑚,𝑛)=(𝑥,𝑦). And then a hyperbolic (saddle) equilibrium results from 𝑉_𝑥𝑥⋅𝑉_𝑦𝑦−𝑉²_𝑥𝑦<0.

Generally, unstable manifolds of the hyperbolic equilibrium may cause instability. The spacecraft needs to be located in the unstable region around the hyperbolic equilibrium in order to meet the requirements of some missions. Thus, it is necessary to establish a stabilizing strategy for the hyperbolic equilibrium.

2.2. Typical Hyperbolic Equilibrium in Astrodynamics

For the typical second-order Hamiltonian system, different 𝑉 and 𝜔 can be used to define different astrodynamical problems, illustrated in the following section.

𝜔=1 and 𝑉=(1/2)(𝑥²+𝑦²)+𝜇/𝑟₁+(1−𝜇)/𝑟₂ define the circular restricted three-body problem (CR3BP) [8], where 𝜇 is the mass ratio between the primaries, 𝑟₁ and 𝑟₂, are, respectively, the distances between the spacecraft and the two primaries, that is, 𝑚₁ and 𝑚₂. There are five libration points, three collinear ones of which are hyperbolic, labeled as 𝐿₁,𝐿₂, and 𝐿₃ shown in Figure 1.

Figure 1 

Hyperbolic equilibria in circular restricted three-body problem: The three collinear equilibria are hyperbolic, labeled as 𝐿1,𝐿2, and 𝐿3. The equilibria in CR3BP are renamed as libration points.

𝜔=1 and 𝑉=(1/2)(𝑥²+𝑦²)+𝜇/𝑟₁+((1−𝜇)/𝑟₂)+𝑎⋅(𝑛_𝑥⋅𝑥+𝑛_𝑦⋅𝑦) define the planar solar sail three body problem [4], where 𝑎 is the solar pressure acceleration and 𝐧=[𝑛_𝑥𝑛_𝑦]^𝑇 is the sail surface normal vector. There are infinite equilibria existing in this system, which can be parameterized by the sail lightness number 𝛽: 𝑟𝛽=⁴¹⋅‖‖1−𝜇−∇𝑉³𝐫₁⋅−∇𝑉², (5) where 𝛽 is the sail lightness number which is the ratio of the radiation pressure acceleration to gravitational acceleration. This parameterization for Sun-Earth/Moon system generates level surfaces as demonstrated in Figure 2.

Figure 2 

Section of the level surfaces in Sun-Earth/Moon system. 𝛽0=0.5059.

𝜔=0 and 𝑈=ℎ²_𝑧/2𝜌²−1/𝑟−𝜅⋅𝑧 define the solar sail’s two-body problem [1], where 𝜅 is the constant determined by solar pressure; the polar coordinates (𝑧,𝜌) are parallel to the sun line and the minimal radius from the 𝑧 axis, respectively; ℎ_𝑧=ℎ_𝑧^max is the constant angular momentum directed along the sun line; 𝑟=𝜌²+𝑧². 𝜔=0 indicates that this Hamiltonian system is degenerated for the absence of the Coriolis acceleration and a unique double equilibrium due to ℎ_𝑧=ℎ_𝑧^max.

Considering the formation flying on a 𝐽₂-perturbed mean circular orbit, the linearized relative dynamics can be deduced as [6, 9] Δ𝐫̈+2𝐅⁻¹𝐅̇Δ𝐫̇+𝐅⁻¹𝐅̈Δ𝐫=−∇_Δ𝐫⋅∇𝑈_{𝐹Δ𝐫=𝟎}⋅Δ𝐫, (6) where Δ𝐫=[𝑥𝑦𝑧]^𝑇 is the relative position in the chief’s LVLH rotating frame, 𝑈 is the gravitational potential function of the deputy including 𝐽₂ perturbation, 𝐅 is the coordinate transformation matrix from the LVLH to inertial frames, and ∇ is the gradient vector described in orthogonal or spherical coordinates. Due to the osculating orbital elements, the relative dynamics on a 𝐽₂-perturbed mean circular orbit is time-periodic Hamiltonian system. This time-periodic system has the fixed equilibrium Δ𝐫=[0,0,0]^𝑇, that is, the location of chief; however, the equilibrium of the system has time-varying topological types and no fixed dimensional unstable/stable/center manifolds, which are quite different from the two-dimensional time-independent system with a permanent pair of hyperbolic eigenvalues and fixed dimensions of unstable/stable/center manifolds, as shown in Figure 3.

(a)

(b)

Figure 3 

Topological type of the equilibrium in the modified uncontrolled system: (a) the time history of real part of the equilibrium during the orbital period; (b) the time history of imaginary part of the equilibrium during the orbital period.

Some averaging techniques can be used to reduce the typical astrodynamics into a first-order Hamiltonian system. For a spacecraft with high area-to-mass ratio orbiting the Earth, its dynamics is strongly perturbed by the term of the gravitational field due to the Earth’s oblateness and by the effect of solar radiation pressure (SRP). The secular rate of the orbital elements due to SRP and 𝐽₂ term can be written as [10] 𝑒̇=−𝐶sin𝜙√1−𝑒²̇√𝜙=−𝐶cos𝜙1−𝑒²𝑒+𝑊1−𝑒²²−1, (7) where 𝐶 and 𝑊 are the constant parameters measuring the solar radiation pressure and the oblateness. The reduced system exists three equilibria, that is, (𝑒₁,0),(𝑒₂,0) and (𝑒₃,𝜋)(𝑒₁<𝑒₂<𝑒₃), where (𝑒₂,0) is hyperbolic and the others are elliptic, as shown in Figure 4.

Figure 4 

Hyperbolic regime in the eccentricity-𝜙 phase space evolution under the effects of 𝐽2+SRP: the bold lines represent the separatrices in correspondence of the saddle point and zero eccentricity Hamiltonian function.

3. Structure-Preserving Stabilization for Hamiltonian System

3.1. Controller Design Preserving Hamiltonian Structure

For the typical second-order Hamiltonian system, the variation equation around the equilibrium is: 𝛿𝐫̈−2𝜔𝐉𝛿𝐫̇−𝑈_𝐫𝐫𝛿𝐫=𝟎. (8) For this unstable system, there exist hyperbolic eigenvalues ±𝜎 associated with the stable and unstable manifolds 𝐮_±, and elliptic eigenvalues ±𝛾𝑖 associated with the centre manifolds 𝐮 and 𝐮.

The Hamiltonian structure-preserving controller then is constructed as 𝐓=−𝜎²𝐺₁𝐮₊𝐮^𝑇₊+𝐺₂𝐮₋𝐮^𝑇₋−𝛾²𝐺₃𝐮𝐮^𝑇+𝐮𝐮^𝑇,𝐓𝐊=2Δ𝐉,_𝐶=𝐓𝛿𝐫+𝐊𝛿𝐫̇, (9)(10) where Δ is used to change the Coriolis acceleration, 𝐺₁,𝐺₂ and 𝐺₃ are, respectively, the gains of unstable, stable, and centre manifolds. The symmetry of matrix 𝐓 and the skew symmetry of matrix 𝐊 guarantee the linear feedback controller preserves the Hamiltonian structure [3].

3.2. Stability Analysis

For the 2-dimensional and time-independent Hamiltonian system, some theorem and proposition about the stability of the controller can be deduced as follows.

Theorem 1. The poles can be assigned at any different positions in imaginary axis, and 𝐺₁,𝐺₂, and 𝐺₃ required are not unique.

Proof. Define the matrix as𝜎𝐀=²1+𝑢²₊𝜎²1+𝑢²₋2𝛾²1+𝑢𝑢𝜎²𝑢²₊1+𝑢²₊𝜎²𝑢²₋1+𝑢²₋2𝛾²𝑢𝑢1+𝑢𝑢𝜎²𝑢₊1+𝑢²₊𝜎²𝑢₋1+𝑢²₋𝛾²(𝑢+𝑢)1+𝑢𝑢, (11) where 𝑢₊ and 𝑢₋ are composing the stable and unstable manifolds of the hyperbolic eigenvalues ±𝜎 as 𝐮_±=11+𝑢²_±1𝑢_± (12) and 𝑢 and 𝑢 are composing the center manifolds of the elliptic eigenvalues ±𝛾𝑖 as 1𝐮=√1+𝑢𝑢1𝑢. (13) Thus 𝑉_𝑥𝑥,𝑉_𝑦𝑦, and 𝑉_𝑥𝑦 will be transformed by the controller: 𝑉^∼_∼𝑥𝑥𝑉_∼𝑦𝑦𝑉_𝑥𝑦=𝑉_𝑥𝑥𝑉_𝑦𝑦𝑉_𝑥𝑦𝐺−𝐀𝐺¹₂𝐺₃, (14) where 𝑉^∼_𝑥𝑥,𝑉^∼_𝑦𝑦, and 𝑉^∼_𝑥𝑦 can be assigned arbitrarily to satisfy the linearly stable conditions, so long as matrix 𝐀 is not regular. Obviously, 𝐀 cannot be regular due to the fact that 𝑢₊ is different from 𝑢₋ for the different eigenvalues +𝜎 and −𝜎.
Denote 𝜆²₁ and 𝜆²₂ as the solutions to the characteristic equation of the second-order dynamics, and then 𝐵 and 𝐶 can be expressed as 𝐵=−𝜆²₁+𝜆²₂,𝐶=𝜆²₁⋅𝜆²₂. (15) Consequently, 𝑉^∼_𝑥𝑥,𝑉^∼_𝑦𝑦, and 𝑉^∼_𝑥𝑦 can be solved from 𝐵=4𝜔²−𝑉_𝑥𝑥−𝑉_𝑦𝑦 and 𝐶=𝑉_𝑥𝑥⋅𝑉_𝑦𝑦−𝑉²_𝑥𝑦.
Furthermore, 𝜆²₁ and 𝜆²₂ are restricted by 0>𝜆²₁1≥−2𝐵≥𝜆²₂, (16) but if one sets 𝜆²₁=𝜆²₂1=−2𝐵 (17) then (𝜆+𝑖𝐵/2)²(𝜆−𝑖𝐵/2)²=0. Since this system has second-order elementary factors, the Jordan form of 𝑉𝟎𝐈_𝐫𝐫2𝜔𝐉 will have long-term dispersions because the matrix cannot be diagonalized.
Because of the nonuniqueness of 𝑉^∼_𝑥𝑥,𝑉^∼_𝑦𝑦, and 𝑉^∼_𝑥𝑦 solved from (9), 𝐺₁,𝐺₂, and 𝐺₃ are nonunique. Suppose that sets of 𝑉^∼_𝑥𝑥,𝑉^∼_𝑦𝑦,𝑉^∼_𝑥𝑦 and 𝑉_𝑥𝑥,𝑉_𝑦𝑦,𝑉_𝑥𝑦 generate the same values of 𝐵 and 𝐶 in (9), and the two sets have the following relationships as 𝑉_𝑥𝑥+𝑉_𝑦𝑦=𝑉^∼_𝑥𝑥+𝑉^∼_𝑦𝑦,𝑉_𝑥𝑥⋅𝑉_𝑦𝑦−𝑉²_𝑥𝑦=𝑉^∼_𝑥𝑥⋅𝑉^∼_𝑦𝑦−𝑉^∼2_𝑥𝑦. (18) It is obtained from (7) that 𝑉²_𝑥𝑥−𝑉^∼_𝑥𝑥+𝑉^∼_𝑦𝑦𝑉_𝑥𝑥+𝑉²_𝑥𝑦+𝑉^∼∼_𝑥𝑥𝑉_𝑦𝑦−𝑉^∼2_𝑥𝑦=0 (19) then the real solution is: 𝑉^∼_𝑥𝑥−𝑉^∼_𝑦𝑦²+4𝑉^∼2_𝑥𝑦−𝑉²_𝑥𝑦≥0 (20) we can fix 𝑉^∼2_𝑥𝑦<𝑉²_𝑥𝑦 to get two different sets generated from 𝐵 and 𝐶. Besides, the fact that there exist different sets generating the same values of 𝐵 and 𝐶 indicate 𝐺₁,𝐺₂, and 𝐺₃ are nonunique.

Hence, the controlled Hamiltonian can be expressed as 1𝐻=2𝐩^𝑇𝐩+𝜔^∼𝐩^𝑇1𝐉𝐪+2^∼𝜔²𝐪^𝑇(𝐪)−1𝐪−𝑉2𝛿𝐪^𝑇𝐓𝛿𝐪. (21) Expand 𝐻 around the equilibrium region to get 𝐻=𝐻₀+12𝛿𝐪^𝑇𝛿𝐩^𝑇𝜔²𝐈−𝑉_𝐫𝐫−𝜔^∼𝐉_∼−1𝜔𝐉𝐈𝛿𝐪𝛿𝐩2𝛿𝐪^𝑇𝐓𝛿𝐪+𝑂(3), (22) where 𝐻₀ represents the Hamiltonian value at the equilibrium and the first-order polynomial disappears because of the equilibrium.

Proposition 2. The Coriolis acceleration can be modified by 𝜔^∼=𝜔+Δ, which cannot stabilize the system independently for 2Δ⋅𝐉𝛿𝐫̇ does not change 𝑈_𝑥𝑥,𝑈_𝑦𝑦, and 𝑈_𝑥𝑦.

The controller developed here can transform the hyperbolic equilibrium (saddle) to an elliptic one (center) according to the theorems and propositions. The elliptic equilibrium has a linear symplectic transformation, which transforms 𝐻 to the following form: 𝐻𝐪^∼∼,𝐩^∼=𝐻₀+12𝜆𝑞₁^∼2₁+𝑝^∼2₁+12𝜆𝑞₂^∼2₂+𝑝^∼2₂+𝑂(3). (23) According to Morse lemma, there exists an analytic diffeomorphism around the equilibrium, which transforms 𝐻^∼ to the following form: 𝐾(𝐪,𝐩)=𝐾𝐪₀,𝐩₀+𝐪^𝑇𝐪+𝐩^𝑇𝐩; (24) thus the energy surfaces are locally diffeomorphic to a family of spheres which shrink down to the equilibria: 𝐾→𝐾(𝐪₀,𝐩₀). The Lyapunov stability follows because the trajectories are tangent to the energy surfaces [11]. Therefore, we have obtained that the nonlinear full dynamics is also Lyapunov stable.

Theorem 3. The gains can be chosen large enough to guarantee the modified elliptic equilibrium is the only equilibrium for the controlled system.

Proof. Since it is easy to demonstrate that 𝐓 is negative definite, the geometric structure of the Hill’s region for the controlled system will be changed by adding −(1/2)(𝐪−𝐪₀)^𝑇⋅𝐓⋅(𝐪−𝐪₀) to the Hill’s region which depends on 𝑉.
Another equilibrium 𝐿_𝑐 of the controlled system (if exist) will stay in Island or Mainland [1] because it must stay in the region where 𝜕_𝑥𝑉<0 and 𝜕_𝑦𝑉<0 for the negative definition of 𝐓.
Here it is just illustrated that the equilibrium 𝐿_𝑐 cannot stay in Mainland, and the negative inclusion for Island can be obtained by the similar way. We can verify whether 𝜕_𝐫𝑉=0 is true by increasing 𝑥 and 𝑦, since 𝐿_𝑐 is the bottom point of Mainland; moreover, we need to test which one descends faster between (1/2)(𝐪−𝐪₀)^𝑇⋅𝐓⋅(𝐪−𝐪₀) and 𝑉. The comparison is very tedious but simply, thus it will be ignored here.
For the controlled system, there exists only the modified elliptic equilibrium, which means the equilibrium has the global minimum pseudopotential. So the controller’s stabilizing region is the whole phase space.

Furthermore, the equilibrium of time-periodic system may be unstable even if the equilibrium is always elliptic during its period, compared with the independent Hamiltonian system. Consequently, it is necessary to verify if the modules of all the multipliers are 1, or if the Floquet multipliers lie on the unit circle in the complex plane. Considering the relative dynamics on a 𝐽₂-perturbed mean circular orbit, the Floquet stability of the periodic Hamiltonian system is obtained by the larger controller gain rather than the critical gain (𝐺₁=𝐺₂=𝐺,𝐺₃=Δ=0), as shown in Figure 5.

(a)

(b)

Figure 5 

Floquet multipliers of the Periodic Hamiltonian System: (a) the relationship between the moduli of Floquet multipliers and the controller gain 𝐺; (b) the Floquet multipliers and the unit circle in the complex plane.

3.3. Resonant Condition

A resonance happens when the eigenvalues of a controlled system have relation 𝜆₁∶𝜆₂=𝑚∶𝑛, where 𝑚<𝑛 and 𝑚 and 𝑛 are reduced with each other. Thus the trajectories will be periodic for any initial condition of spacecraft. We will give a particular solution for the resonance in the following section.

Suppose 𝜆₁∶𝜆₂=𝑚∶𝑛, and then: 𝑚²𝐵+𝐵²−4𝐶=𝑛²𝐵−𝐵²−4𝐶 (25) hence, 𝑚𝑛4𝜔^∼2−𝑉^∼_𝑥𝑥−𝑉^∼_𝑦𝑦=𝑚²+𝑛𝑉^2∼∼_𝑥𝑥𝑉_𝑦𝑦−𝑉^∼2_𝑥𝑦. (26) A particular solution can be constructed as 𝑉^∼_𝑥𝑥=𝑉^∼_𝑦𝑦=𝜔^∼2,_∼𝑉_𝑥𝑦=±𝜔^∼2⋅𝑛²−𝑚²𝑚²+𝑛². (27) If 𝜔^∼=0, a particular solution can be constructed as 𝑉^∼_𝑥𝑥=𝑉^∼_𝑦𝑦=−²_∼𝜔𝑉_𝑥𝑦=±𝜔²⋅𝑛²−𝑚²𝑚²+𝑛², (28) where 𝜔 is an arbitrary parameter.

Compared with the osculating frequencies of the controlled elliptic equilibrium for the global investigation during its period, the characteristic frequencies defined by the Floquet multipliers have more applications for the time-periodic Hamiltonian system. Floquet multipliers 𝑒^{±𝜃𝑘𝑖}(𝑘=1,2,3) are the eigenvalues of the linear Poincaré mapping [Δ𝐫,Δ𝐫̇]^𝑇(𝑡)→[Δ𝐫,Δ𝐫̇]^𝑇(𝑡+𝑇_𝑑) with the unit eigenvectors 𝝂₁𝝂₁,𝝂₂𝝂₂, and 𝝂₃𝝂₃. Characteristic frequencies 𝜃₂ and 𝜃₃ are approximately equal to zeros inheriting from two of osculating frequencies around 𝜛, while 𝜃₁ depends on the controller gain. The linear combinations of eigenvectors, (𝝂₁+𝝂₁),(𝝂₁−𝝂₁)𝑖,(𝝂₂+𝝂₂),(𝝂₂−𝝂₂)𝑖,(𝝂₃+𝝂₃), and (𝝂₃−𝝂₃)𝑖, are utilized to generate six different types of near-periodic relative orbits with their periods of (2𝜋/𝜃₁)𝑇_𝑑 and 𝑇_𝑑, respectively. The six different types of bounded relative orbits generated by the linear combinations can be then regarded as the foundational motions around the equilibrium.

Particularly, the linear combinations (𝝂₁+𝝂₁) and (𝝂₁−𝝂₁)𝑖 can be used as the initial values to generate the periodic relative orbits with their periods of 𝑁⋅𝑇_𝑑 when 𝜃₁ satisfies the resonant condition 𝑁⋅𝜃₁=2𝜋. To be stricter, all the relative periodic orbits, which produced by the judicious choice combined of the eigenvectors, are near periodic only for the nonlinear term that ignored by the linearized differential equations describing the relative dynamics.

It is worth mentioning that the six eigenvectors of the Floquet multipliers span the entire space of the relative position and velocity [Δ𝐫^𝑇,Δ𝐫̇^𝑇]^𝑇. Therefore, any initial relative positions and velocity vectors can be decomposed as [Δ𝐫^𝑇,Δ𝐫̇^𝑇]^𝑇=∑⁶_𝑖=1𝛼_𝑖𝝂_𝑖, where 𝛼_𝑖 is constant coefficients, and 𝝂_𝑖 is chosen among (𝝂₁+𝝂₁),(𝝂₁−𝝂₁)𝑖,(𝝂₂+𝝂₂),(𝝂₂−𝝂₂)𝑖,(𝝂₃+𝝂₃), and (𝝂₃−𝝂₃)𝑖. Hence, all the initial relative position and velocity can generate a bounded trajectory around the controlled elliptic equilibrium, which depends on the topological type of the equilibrium. Moreover, all the general trajectories are quasiperiodic for the different frequencies between the foundational motions and all the bounded trajectories are involved in the foundational motions caused by the eigenvectors.

3.4. Cost and Optimization

Because 𝐺₁,𝐺₂,𝐺₃, and Δ are nonunique, it is necessary to investigate the rule of choosing 𝐺₁,𝐺₂,𝐺₃, and Δ to stabilize the system. Denote 𝑈_𝜆1,𝜆2 as the collection of 𝐺₁,𝐺₂,𝐺₃, and Δ which will allocate the controlled system with the expected poles 𝜆₁ and 𝜆₂.

The controller’s output can be then determined as 𝐓_𝐶̇𝛿𝐫̇=𝐓⋅𝛿𝐫+𝐊⋅𝛿𝐫=𝐓𝟎𝟎𝐊𝛿𝐫. (29) Define 𝐳=[𝛿𝐫^𝑇𝛿𝐫̇^𝑇]^𝑇 to specify the sensitivity of the controller as 𝜅=‖𝐓_𝐶‖₂/‖𝐳‖₂.

The Frobenius norm which is consistent to the Euclidean norm of the vector, can be used to measure the controller’s sensitivity as 𝐓𝜅=^𝐶2‖𝐳‖₂≤𝐓𝟎𝟎𝐊_𝐹=tr𝐓^𝑇𝐓+8Δ². (30)

Define 𝐆_𝑚=𝐺₁𝐺₂𝐺₃𝐺₁𝐺₂𝐺₃𝐺₁𝐺₂𝐺₃𝐺₁𝐺₂𝐺₃^𝑇,𝐀_𝑚=𝐀₁₁𝐀₁₂𝐀₁₃𝐀₃₁𝐀₃₂𝐀₃₃𝐀₃₁𝐀₃₂𝐀₃₃𝐀₂₁𝐀₂₂𝐀₂₃ (31) and then deduce the matrix 𝐓 from 𝐆_𝑚 and 𝐀_𝑚 as𝐓=^{𝐺1⋅𝜎21+𝑢2++𝐺2⋅𝜎21+𝑢2−+𝐺3⋅2𝛾21+𝑢𝑢𝐺1⋅𝜎2𝑢+1+𝑢2++𝐺2⋅𝜎2𝑢−1+𝑢2−+𝐺3⋅𝛾2(𝑢+𝑢)1+𝑢𝑢}_{𝐺1⋅𝜎2𝑢+1+𝑢2++𝐺2⋅𝜎2𝑢−1+𝑢2−+𝐺3⋅𝛾2(𝑢+𝑢)1+𝑢𝑢𝐺1⋅𝜎2𝑢2+1+𝑢2++𝐺2⋅𝜎2𝑢2−1+𝑢2−+𝐺3⋅2𝛾2𝑢𝑢1+𝑢𝑢}=𝐀_𝑚⋅𝐆_𝑚 (32)

which has the following relationship as: ‖𝐓‖²_𝐹=tr𝐓^𝑇𝐓²≤𝐀_𝑚²_𝐹⋅𝐆_𝑚²_𝐹=𝐀_𝑚²_𝐹⋅4𝐺²₁+𝐺²₂+𝐺²₃. (33) Thus, put it into (33) and (30) to obtain: 𝜅≤𝐀_𝑚²_𝐹⋅4𝐺²₁+𝐺²₂+𝐺²₃+8Δ², (34) where ‖𝐀_𝑚‖_𝐹 can be calculated directly from the constant matrix 𝐀_𝑚 according to the invariant manifolds of the equilibrium.

For full feedback, 𝐿=‖𝐀_𝑚‖²_𝐹⋅4(𝐺²₁+𝐺²₂+𝐺²₃)+8Δ² can be treated as the optimization index to select more suitable values for the manifolds gains; however, for the position feedback, the optimization index can be set as 𝐿=𝐺²₁+𝐺²₂+𝐺²₃. With the constraint equation (15), the selection for controller’s gains can be considered as a nonlinear constrained optimization, as: min_{𝐺1,𝐺2,𝐺3∈,Δ𝑈𝜆12,𝜆}𝜆𝐿,s.t.𝐵+²₁+𝜆²₂=0,𝐶−𝜆²₁⋅𝜆²₂=0. (35) The programming functions “fmincon” and “confuneq” in MATLAB can solve this optimization problem.

As confirmed by zero real part of the eigenvalues in Figure 6, the equilibrium is elliptic during the period for time-periodic Hamiltonian system. Therefore, there exist three osculating imaginary eigenvalues ±𝜔₁(𝑡)𝑖,±𝜔₂(𝑡)𝑖, and ±𝜔₃(𝑡)𝑖 for the controlled equilibrium. The quasiperiodic relative trajectories around the equilibrium comprise three different oscillating motions with the osculating frequencies of 𝜔₁,𝜔₂, and 𝜔₃, which can be verified by the second-order term in the new Hamiltonian function yielded by a linear symplectic transformation: 𝐻𝐪^∼∼,𝐩^∼,𝑡=𝐻^∼₀()+1𝟎,𝟎,𝑡2𝜔₁(𝑡)𝑞^∼2₁+𝑝^∼2₁+12𝜔₂(𝑡)𝑞^∼2₂+𝑝^∼2₂+12𝜔₃(𝑡)𝑞^∼2₃+𝑝^∼2₃+𝑂(3). (36)

(a)

(b)

Figure 6 

Topological type of the uncontrolled equilibrium: (a) the time history of real part of the equilibrium during the orbital period; (b) the time history of imaginary part of the equilibrium during the orbital period.

According to the local optimization, less cost is obtained by small controller gain illustrated by (34); however, the smaller gain cannot maintain the relative position Δ𝐫 in a domain closer to the equilibrium, which may consume more fuels.

According to the global optimization, an averaging quadratic cost function is defined to measure the fuel consumption during [0,𝑡_𝑓] for the continuous thrust engine [12]: 1𝐽=2𝑡_𝑓^𝑡𝑓₀𝐓^𝑇_𝑐𝐓_𝑐1𝑑𝑡=2𝑡_𝑓^𝑡𝑓₀Δ𝐫^𝑇𝐆^𝑇𝐂^𝑇(𝑡)𝐂(𝑡)=1𝐆Δ𝐫𝑑𝑡2𝑡_𝑓⁶⁶_𝑖=1𝑗=1𝛼_𝑖𝛼_𝑗^𝑡𝑓₀𝜻_𝑖(𝑡)^𝑇𝐆^𝑇𝐂^𝑇(𝑡)𝐂(𝑡)𝐆𝜻_𝑖(𝑡)𝑑𝑡, (37) where 𝜻_𝑖(𝑡) is the position component of the periodic orbit which is developed by the eigenvector 𝝂_𝑖. The component in the 𝑖th row and 𝑗th column of the measuring matrix 𝐌(𝐺) is 𝐌_𝑖𝑗=(1/𝑡_𝑓)∫^𝑡𝑓₀[𝜻_𝑖(𝑡)^𝑇𝐆^𝑇𝐂^𝑇(𝑡)𝐂(𝑡)𝐆𝜻_𝑖(𝑡)]𝑑𝑡, which inherits the near-periodicity from the characteristic frequencies and depends only on the controlled elliptic equilibrium. Therefore, 𝐌 inherited from the longest period of the characteristic frequency 𝜃₁, that is, 𝑡_𝑓=(2𝜋/𝜃₁)𝑇_𝑑, can serve as the foundational measuring matrix that can be obtained offline (or on the ground) from the linearized dynamics for the fast and approximative estimation. The consumption function during [0,(2𝜋/𝜃₁)𝑇_𝑑] has the form as: 𝐽()=1𝐺,𝜶2𝜶^𝑇𝐌(𝐺)𝜶, (38) where 𝐌 is symmetric and positive-definite and 𝜶=[𝛼₁𝛼₂𝛼₃𝛼₄𝛼₅𝛼₆]^𝑇. Generally, there exists a unitary matrix 𝐔 from Singular Value Decomposition (SVD) to transform the function as 𝐽=(1/2)∑⁶_𝑖=1𝜆_𝑖𝛽²_𝑖, where 𝜆_𝑖 represents the eigenvalues of 𝐌 and 𝛽_𝑖 is the 𝑖th column element of the intermediate variable 𝜷=𝐔𝜶.

The identified gain can be derived to minimize 𝐽 if we know the initial relative position and velocity. The specific gain satisfying the minimal trace of 𝐌(𝐺), that is, tr𝐌=∑⁶_𝑖=1𝜆_𝑖, is however preferable to minimize 𝐽 for the maximum likelihood, for an unknown case on the initial conditions.

4. Application to Nonlinear Astrodynamics

4.1. Stable Lissajous Orbits of Solar Sail: Application to Second-Order Hamiltonian System

A solar sail is a new type of spacecraft with no fuel, and its orbital control is realized by the solar radiation pressure by means of orientating its attitude relative to the Sun. Therefore, we can neglect the control consumption and implement the complicated control law, which is incomparable to the low-trust ion propulsion or the impulse propulsion.

We noticed the fact that there exists quasiperiodic or periodic orbits near the sail’s equilibria, and several periodic orbits have been proposed in the solar sail restricted three-body problem. A survey [13] aimed to summarize the combined wealth of literature concerned with the dynamics, stability and control of highly non-Keplerian orbits for various low thrust propulsion devices, and to demonstrate some of these potential applications, like Capturing near earth objects [14]. McInnes [15] is the first person who applied the classical libration point orbit theory to generate the halo orbit around the on-axis sail’s equilibrium. Baoyin and McInnes [16] achieved two different halo orbits types around the on-axis sail’s equilibrium. McInnes [15] has tried to generate lissajous orbits near the equilibrium as well, but failed to keep a bounded trajectory for a long time. Waters and McInnes [17] applied the classical theory to generate some halo orbits that is only available for the off-axis sail equilibria in some specified regions. What is more, all the interesting orbits presented by researchers are unstable and have one-dim unstable manifolds, hence the station-keeping strategies becomes quite necessary.

The sail associated with the Sun and Earth-Moon system, (where the Earth and Moon are regarded as a whole celestial body, denoted as Earth/Moon), is regarded as the solar sail restricted three-body system. In this section we suppose that the Sun-Earth/Moon system revolves around a circular orbit.

The structure-preserving controller is proposed to generate stable lissajous orbits and Lyapunov orbits for nonresonance shown in Figure 7, and resonant periodic orbits shown in Figure 8 (1 : 2, 1 : 3, 1 : 4, 2 : 3, resp.).

Figure 7 

Stable Lissajous orbit and Lyapunov orbits near the sail equilibrium: 𝐺1=20,𝐺2=10,𝐺3=10.

Figure 8 

Resonant periodic orbits near the sail equilibrium.

The feedback gains are initially chosen with 𝐺₁=20,𝐺₂=10,𝐺₃=10 and Δ=0 in view of the local optimization, and then the gains will be optimized so as to obtain the same poles. The gains refined by optimization are 𝐺₁=13.0147,𝐺₂=13.0147,𝐺₃=12.0063, and Δ=0, with the optimization index improved from 600 to 483 (80.5%) for the position feedback. While for full feedback, the gains refined by optimization are 𝐺₁=11.3253,𝐺₂=11.3254,𝐺₃=14.4533, and Δ=−0.5956, with the optimization index improved from 2.0058×10³ to 1.5587×10³ (77.71%).

4.2. Bounded Orbits by Degenerated Solar Sail: Application to Degenerated Hamiltonian System

McInnes [18] has researched the displaced non-Keplerian orbits above the Sun. Similar steps have been taken by Bookless and McInnes [19, 20] to the displaced non-Keplerian orbits above the planet. Nevertheless, just as the stability analysis for motions around the equilibrium, Bookless’ analysis for the dynamics and control is linear and local. There presented the nonlinear dynamical analysis for a displaced orbit above a planet, and there also investigated the motions around the equilibria for the nonresonance case with the help of the Birkhoff normal form and dynamical system techniques. M. Xu and S. Xu [1] has obtained two important contributions to derive necessary and sufficiency conditions for the motion stability around the equilibria, and to illustrate numerically that the non-transition critical KAM torus is filled with the Lyapunov (1,1)-homoclinic orbits. ℎ_𝑧=ℎ_𝑧^max defines the degradation case for the two body problem of solar sails, which is the saddle-node bifurcation point.

For the degenerated system, there exists only one double equilibrium, and the frequencies for the liberalized motion are, respectively, 𝜔₁=0 and 𝜔₂>0. What is more, it is demonstrated from the variation equation around the equilibrium that the equilibrium has no stable or unstable manifolds but centre manifolds. Owing to Theorem 1, it guarantees that the stabilization for the degenerated ℎ_𝑧^max can be achieved by the equilibrium centre manifolds.

Similar to the nondegradation case, the suitable gains can be chosen to guarantee the modified elliptic equilibrium is the only equilibrium for the controlled system. So the controller’s stabilizing region is all the phase space. The quasiperiodic orbits generated by the controller are shown in Figure 9, and their positions in 3D physical space are shown in Figure 10.

Figure 9 

Quasiperiodic orbits by the controller near the degenerated equilibrium: 𝐺1=0,𝐺2=0,𝐺3=0.3,Δ=0.0005.

(a)

(b)

(c)

(d)

Figure 10 

Quasiperiodic orbits by the controller in the 3D physical space: 𝐺1=0,𝐺2=0,𝐺3=0.3,Δ=0.0005.

4.3. Quasiperiodic Formation Flying: Application to Time-Periodic Hamiltonian System

Roughly, we can classify the relative motion control approaches into two branches: continuous control that using low-thrust electric propulsion and impulsive control that relying on chemical thrusters. While the former approaches used mostly full-state feedback of the Cartesian relative position and velocity to develop high-accuracy tracking control laws, thus capable to maintain small steady-state errors under a myriad of orbital perturbations [21–25]; however, the maintenance of the later assumes that the tracking reference orbital elements are known and attempts to generate control commands that will match the osculating orbital elements to the expected values [12, 26, 27].

The initial relative position of 𝑥₀=𝑦₀=𝑧₀=1km and the relative velocity generated from the classical 𝐶-𝑊 equations 𝑦̇₀=−2𝜛𝑥₀,𝑥̇₀=𝜛𝑦₀/2,𝑧̇₀=0 for the formation flying on a 𝐽₂-perturbed mean circular orbit, are applied to generate the quasiperiodic relative trajectories with the preferred gain 0.31 illustrated in Figure 11. More rapid changes in configuration are obtained by the controller developed in this paper, rather than by the change in 𝐽₂ invariant relative orbits caused by the secular term of perigee argument, from the projection of the configuration on the 𝑥-𝑧 or 𝑦-𝑧 plane. Therefore, the mission that requires rapid changes in formation configuration, like on-orbit inspection and repair could potentially use the addressed controller. Furthermore, the quasiperiodic relative trajectories which comprise six foundational motions with different frequencies have no reference relative configuration to track for the continuous thrust, make the trajectory prediction difficult by other nonallied spacecrafts.

(a)

(b)

(c)

(d)

Figure 11 

Quasiperiodic relative trajectories generated by the optimal gain of 0.31 during 30 days: presented is the relative orbit during the first orbital period in every day; the black thicker line indicates the initial formation configuration.

Figure 12 illustrates the topological type of the equilibrium of the controlled system, which is elliptic during the period, as confirmed by zero real part of the eigenvalues, for the mean circular orbit whose mean inclination is 80°, mean semi-major axis is 7178.137 km, and initial argument of latitude is 0°.

(a)

(b)

Figure 12 

Topological type of the equilibrium in the controlled system: (a) the time history of real part of the equilibrium during the orbital period; (b) the time history of imaginary part of the equilibrium during the orbital period; on the mean circular orbit with 𝑎=7178.137km,𝑖 = 80°, and an initial argument of latitude of 0°.

𝐺₊ should share its value with 𝐺₋ to minimize the controller outputs from view of local optimization, according to the contour curves of 𝜅 illustrated in Figure 13. The specific gain owning the minimal trace of 𝐌(𝐺), that is, tr𝐌=∑⁶_𝑖=1𝜆_𝑖, is preferable to minimize 𝐽 for the maximum likelihood if the initial conditions are unknown. The optimizing result obtained from local optimization is illustrated in Figure 14 that the optimal gain of 0.31 has a smaller trace than the critical gain of 0.27.

Figure 13 

Stable regions of the controller in the (𝐺+-𝐺−) space: the hyperbolic-like curve 𝐶1 is the boundary keeping the topological type of the equilibrium elliptic; the line-like curve 𝐶2 is the boundary keeping the Floquet multipliers lying on the unit circle in the complex plane; the closed curves labeled by numbers are contour curves of the averaging values of 𝜅; the circle 𝐶3 is the approximation of the contour curve.

Figure 14 

Trace of the measuring matrix: compared with the trace of the critical gain of 0.27 obtained from local optimization, the gain of 0.31 can achieve the minimal trace.

4.4. Controlled Frozen Orbits: Application to First-Order Hamiltonian System

The high area-to-mass spacecraft’s orbital evolution behaves strangely under the influence of the solar radiation pressure and the perturbation due to the oblateness of the Earth. Since 1960, there have been observed long-term oscillations of the eccentricity in the satellites’ orbital behaviours, like Vanguard and the ECHO balloon. The two variables, the osculating orbital eccentricity and the solar angle 𝜙 between the orbit pericentre and the Sun-line [10, 28], can be used to describe this dynamical system with low inclinations through a Hamiltonian function.

All spacecrafts around the hyperbolic eccentricity cannot be maintained on such an orbit; on the contrary, they will move on the unstable manifold towards higher eccentricities, or on the unstable manifold towards lower eccentricities based on the quite small change of the orbital elements. Controlled frozen trajectories can be obtained in the phase space around the hyperbolic eccentricity by means of applying the Hamiltonian structure-preserving controller, such as the one illustrated by red dotted line in Figure 15.

(a)

(b)

Figure 15 

Controlled frozen orbits generated by the controller in 𝑒-𝜙 space.

The mission illustrated in Figure 15(a) comprises three phases: (i) transferring (blue line) from a circular orbit to the hyperbolic equilibrium, (ii) stabilizing near the hyperbolic equilibrium (red line), and (iii) transferring (dotted blue line) beyond the critical eccentricity. While the mission illustrated in Figure 15(b) comprises three phases as well: (i) transferring from a circular orbit to the hyperbolic equilibrium (blue line), (ii) stabilizing (red line) near the hyperbolic equilibrium, and (iii) transferring (the dotted blue line) back to the circular orbit.

5. Conclusion

Bounded motions near hyperbolic equilibria have been broadly applied to various astronautical missions. The Hamiltonian structure-preserving controller can change the hyperbolic equilibrium (saddle) to an elliptic one (center), and then some quasiperiodic bounded trajectories will emerge naturally from the KAM tori generated by the controller.

In this paper, a systematic investigation is presented on the general form for full feedback and position-only feedback modes, Lyapunov and Floquet stability analysis, and control gain selection of the structure-preserving controller, in the context of second-order, first-order, time-periodic, and degenerated Hamiltonian system, respectively. One type of periodic trajectories is achieved by the resonant conditions of control gains for the resonant periodic orbits, and another type is making judicious choice in the foundational motions with different frequencies. Therefore, the controller has potential applications in stable lissajous orbits for solar sail’s three-body problem and degenerated two-body problem, quasiperiodic formation flying on a 𝐽₂-perturbed mean circular orbit, or the controlled frozen orbits for a spacecraft with high area-to-mass ratio.

Furthermore, some further researches can be implemented, including how to use unfixed dimensional center manifolds to feedback the controller, the further investigation for the global Floquet stable controller gain for quasiperiodic Hamiltonian system, and the relationships between the different frequencies and the critical gain for Floquet stability or the optimal gain.

Acknowledgments

The research is supported by the National Natural Science Foundation of China (11172020), the National High Technology Research and Development Program of China (863 Program: 2012AA120601), Talent Foundation supported by the Fundamental Research Funds for the Central Universities, Aerospace Science and Technology Innovation Foundation of China Aerospace Science Corporation, and Innovation Fund of China Academy of Space Technology.