Open access peer-reviewed chapter

Analysis and Control of Nonlinear Attitude Motion of Gravity-Gradient Stabilized Spacecraft via Lyapunov-Floquet Transformation and Normal Forms

Written By

Peter M.B. Waswa and Sangram Redkar

Submitted: 21 March 2019 Reviewed: 10 June 2019 Published: 15 January 2020

DOI: 10.5772/intechopen.87954

From the Edited Volume

Advances in Spacecraft Attitude Control

Edited by Timothy Sands

Chapter metrics overview

987 Chapter Downloads

View Full Metrics

Abstract

This chapter demonstrates analysis and control of the attitude motion of a gravity-gradient stabilized spacecraft in eccentric orbit. The attitude motion is modeled by nonlinear planar pitch dynamics with periodic coefficients and additionally subjected to external periodic excitation. Consequently, using system state augmentation, Lyapunov-Floquet (L-F) transformation, and normal form simplification, we convert the unwieldy attitude dynamics into relatively more amenable schemes for motion analysis and control law development. We analyze the dynamical system’s periodicity, stability, resonance, and chaos via numerous nonlinear dynamic theory techniques facilitated by intuitive system state augmentation and Lyapunov-Floquet transformation. Versal deformation of the normal forms is constructed to investigate the bifurcation behavior of the dynamical system. Outcome from the analysis indicates that the motion is quasi-periodic, chaotic, librational, and undergoing a Hopf bifurcation in the small neighborhood of the critical point-engendering locally stable limit cycles. Consequently, we demonstrate the implementation of linear and nonlinear control laws (i.e., bifurcation and sliding mode control laws) on the relatively acquiescent transformed attitude dynamics. By employing a two-pronged approach, the quasiperiodic planar motion is independently shown to be stabilizable via the nonlinear control approaches.

Keywords

  • gravity gradient
  • nonlinear attitude control
  • sliding mode
  • Lyapunov-Floquet transformation
  • normal forms

1. Introduction

Ever since the launch of Sputnik—the first artificial satellite put into earth’s orbit in 1957—mankind has increasingly become dependent on space-based technology in many areas of our daily lives. For instance, space-based technology performs an indispensable role in telecommunications, navigation, personal entertainment, weather forecasting, farming, security, defense, scientific exploration, research, innovation, etc. Undoubtedly, the prominence of space technology in shaping humanity’s future is unequivocal. The success or failure of a given space mission is largely contingent upon the complex system analysis and design methodologies exerted in converting the initial idea into an elaborate functioning enterprise [1]. It is for this reason that reliable and efficacious methodologies and tools are consistently utilized in space mission formulation and implementation. Thus, there is a need to continuously examine the effectiveness of prevailing space mission analysis and design methodologies. This is in order to improve prevailing tools and approaches that shall expedite relatively simpler, more reliable, and accurate mission modeling and analysis.

Space systems are required to function nominally in their designated orbital locations, maintain appropriate orientation, and conform to planned trajectories despite the ambient perturbing space environment. Strict mission pointing requirements normally constrain spacecraft in orbit around a large body to maintain a fixed stable orbital position and orientation during operation. However, perturbing space-environment torques act to dislodge positioned spacecraft and disorient stabilized ones [2, 3].

Modeling, analyzing, and controlling dynamics of space systems are therefore a crucial component of space mission design. The quest for relatively simpler, more accurate, and more reliable analytical methodologies and tools to represent, scrutinize, and manipulate the dynamics of space systems is therefore a worthwhile undertaking.

Inopportunely, dynamics of space systems tend to be commonly represented by coupled analytical models that possess complex structures encompassing nonlinearity, parameter-variant coefficients, and periodic external excitation terms [4, 5, 6, 7, 8, 9]. The requisite analysis essential to fathom such motion is not a trivial undertaking—except for few special cases, the general solution for such dynamical systems cannot be found. The complex structures of the motion’s analytical models characteristically point to nondeterministic and potentially chaotic systems over a range of initial conditions and system parameters. Therefore, to analyze dynamical space systems, we often have to be content with nonautonomous, nonlinear, and periodic differential equations [10, 11]. This presents an immense analysis challenge. For instance, time-varying eigenvalues of the periodic linear system matrix cannot determine the system stability. Consequently, methods such as linearization [12], averaging [13], and perturbation techniques [14, 15] have been consistently used to analyze such complex nonlinear, periodic motion. However, the two latter approaches tend to be limited to minimally excited systems (parameter multiplying the periodic terms is small), while linearization is restricted to small domains about the operating point. Further, such methods are inclined to be relatively cumbersome and normally augmented with numerical approaches to analyze dynamical systems and accomplish real-life applications [16, 17, 18].

The presence of perturbing torques in the ambient space environment tends to disorient an already stabilized spacecraft and further alters the orbital motion [2, 3]. This is contrary to the prevailing strict pointing requirements that constrain the satellite in orbit to maintain a stable orbital motion and fixed orientation during mission operation. A number of strategies are employed to stabilize spacecraft attitude motion and maintain a desired orientation despite the presence of perturbing torques in the space environment. The most common attitude control and stabilization approaches are three-axis stabilization, spin stabilization, and gravity-gradient stabilization. To provide the determined control input required to offset undesired attitude deviations, these methods employ either active control systems (e.g., thrusters, magnetic torquers, reaction wheels) or passive control systems (e.g., booms). Unlike passive controllers, active controllers utilize an external source of energy to drive the attitude control actuators [17, 19].

Among the stated attitude stabilization methods, gravity-gradient stabilization of spacecraft attitude is attractive due to its relatively intrinsic simplicity, reliability, and low cost [20]. However, it is mostly feasible in low earth orbit due to its principle of operation as discussed in Section 2.

The motion about COM of a rigid gravity-gradient stabilized spacecraft is libratory about the pitch axis. This axis is normal to the orbital plane in an inverse-square gravity field. The satellite will oscillate about a position of stable relative equilibrium if the work done by external perturbing forces is greater than the rotational kinetic energy. The sufficient conditions for stability of relative equilibrium are explained in Section 2. The complete formulation of COM motion for a gravity-gradient stabilized satellite in eccentric orbit consists of six coupled, nonlinear second-order differential equations. This system of equations is considered analytically unsolvable in closed form [5, 7, 21, 22, 23].

This chapter aims to first investigate the periodicity, quasi-periodicity, and chaotic behavior of the gravity-gradient stabilized attitude motion. Moreover, the motion stability, resonances and bifurcation behavior will also be examined. Subsequently in Section 4, we synthesize suitable controllers to adequately offset the attitude perturbations experienced by the gravity-gradient stabilized spacecraft in an eccentric orbit. Requisite assumptions made to facilitate the attitude motion analyses will be explicitly stated and qualified.

To model, analyze, and control the nonlinear motion with parameter-variant coefficients and periodic forcing terms, we intend to use approaches based on:

  • System state augmentation

  • Lyapunov-Floquet (L-F) transformations

  • Normal form (NF) theory

The fitting use of the aforementioned transformations and techniques enables dynamical system analysis and control law development in transformed, parameter-invariant, and more tractable coordinates that preserve the original system Lyapunov stability properties [24, 25]. Consequently, we intend to exploit this propitious attribute in our investigation. Applications of L-F transformations in spacecraft dynamics have been previously investigated by authors such as [8, 26]. The former demonstrates how L-F theory enhances the representation of relative spacecraft dynamics in elliptical orbits, while the latter further proposes an orbit control law based on L-F theory. On the other hand, this chapter focuses on the dynamics of a rigid body about its COM while in elliptical orbit around a central large mass.

2. Gravity-gradient attitude stabilization in eccentric orbit

A gravity-gradient stabilized spacecraft attains a state of stable relative equilibrium when its I z points along the radius vector, I y points along the normal to the orbit plane, and I x is along the tangent to the orbit in the LVLH frame (Hill frame) as shown in Figure 1a. In addition, the condition I y > I x > I z must similarly be satisfied.

Figure 1.

Gravity-gradient stabilization. (a) Geometry of orbit and attiutde parameters, and (b) Pitch angle librations.

If the work done by external perturbing torques is greater than the rotational kinetic energy of the spacecraft about its COM, motion of the spacecraft in an elliptical orbit will be libratory as illustrated in Figure 1b. Equations representing the spacecraft orbital motion are identical with those of a point mass in an inverse-square law force field. To analyze the attitude dynamics, the spacecraft orbital motion (motion of COM) can be reasonably assumed to be independent of the spacecraft attitude motion (motion about COM). This assumption is justifiable because the satellite is small compared to the dimensions of the orbit. Under this assumption, the spacecraft’s orbital motion can hence transfer energy to the attitude motion, but the converse is assumed not to be possible. Thus, orbital parameters as determined functions of time are considered in analyzing attitude motion [7, 27].

When the spacecraft is considered as a rigid body in an inverse-square force field along an elliptical orbit, a complete formulation of equations of motion can be derived [5, 7, 21]. The resulting six second-order differential equations of motion are nonlinear and coupled. These equations of motion cannot be solved analytically in this exact form.

Ignoring other torques such as aerodynamic, magnetic, thermal bending, and solar radiation pressure, we can derive the equations of spacecraft attitude motion under the influence of inverse-square force field in an elliptical orbit. Additional assumptions are an ideal, perfect sphere earth without oblateness; largest spacecraft dimensions are extremely small compared to the orbit radius, and the spacecraft mass is negligible compared to the mass of the central body [27].

We further assume that the exact equations of motion can be linearized in small-angle motion characterization. Subsequently, the attitude dynamic models may be considered to consist of two equations with coupled roll-yaw angles and a third uncoupled equation describing the pitch angle dynamics. The pitch motion equation is hence independent of roll-yaw motion. The coupled roll-yaw ( Ψ - Ω ) equations are homogeneous and can be solved for Ψ = Ψ ̇ = Ω = Ω ̇ = 0 [5, 7, 21, 27].

Consequently, the exact problem is reduced to the equation of pitch motion with orbital parameters as functions of time and spacecraft mass parameters shown in Eq. (1):

Θ ¨ + 3 μ r 3 σ sin Θ cos Θ = ω ̇ E1

where 0 σ 1 is a dimensionless ratio of the spacecraft’s principal moments of inertia given by

σ = I x I z I y = I roll I yaw I pitch E2

To analyze attitude motion in eccentric orbit, we substitute time with true anomaly, where f is the independent variable. Moreover, the COM will obey the following Keplerian orbit relations:

r = P 1 + e cos f E3
ω = df dt = μP r 2 = μP P 2 1 + e cos f 2 E4

Thus

ω ̇ = 2 μ r 3 e sin f E5
Θ ¨ = 1 + e cos f μ r 3 d 2 Θ df 2 2 μ r 3 e sin f d Θ df E6

Substituting Eqs. (5) and (6) into Eq. (1) yields

1 + e cos f Θ 2 e sin f Θ + 3 σ sin Θ cos Θ = 2 e sin f E7

This is the well-known equation of plane pitch angle libratory motion in elliptical orbit [5, 7, 21, 27]. The primes indicate differentiation with respect to f. The planar pitch attitude motion equation is hence nonlinear with periodic coefficients in f. Analysis of this motion and subsequent synthesis of a fitting controller is not a trivial task. We hence intend to analyze this motion and synthesize suitable controllers to stabilize the system.

3. Attitude motion analysis

In general, L-F transformation techniques facilitate obtaining solutions of dynamical systems with periodic coefficients, evaluate periodically forced responses, and design feedback control laws. We shall augment these capabilities with normal form techniques (simplifies nonlinearity) and state augmentation (converts nonautonomous to autonomous system). Subsequently, the emanating synergies serve to accomplish the objectives of this chapter.

Floquet theory enables stability and response analysis of linear systems with periodic coefficients, i.e., x ̇ t = A t x t , s . t A t = A t + T . Floquet theory utilizes knowledge of characteristic exponents of the state transition matrix (STM) to infer that if the solution of a system is obtained over a full principle period, then the solution is known for all time [28, 29].

Further, the the Lyapunov-Floquet (L-F) transformation x t = Q t z t reduces the original nonautonomous linear differential system to an autonomous one with the form z ̇ t = Rz t where R is an n × n constant matrix. In the LFT matrix, Q t can be approximated via the methodology described by [25, 30] using shifted Chebyshev polynomials of the first kind. Chebyshev polynomials are chosen because they produce better approximation and convergence than other special orthogonal polynomials [31]. The approximation of LFT matrix via Chebyshev polynomials contains elements Q ij t with truncated Fourier representations as shown in Eq. (8). Q ij 1 t has similar Fourier representation:

Q ij t a 0 2 + n = 1 q a n cos πnt T + n = 1 q b n sin πnt T E8

Time-independent normal form (TINF) simplification facilitates construction of relatively lesser complex but qualitatively equivalent models of the original nonlinear dynamical systems. On the other hand, time-dependent normal form (TDNF) simplification considers original nonlinear dynamical systems with periodic coefficients by utilizing Lyapunov-Floquet transformation. This simplification is generally implemented on equations arising from Taylor series expansion via successive application of nonlinear near-identity transformations. Such a transformation entails preservation of the original system’s stability and bifurcation characteristics by the transformed models. The fundamental concept behind normal forms methodology is to simplify the system by eliminating as many nonlinear terms as possible. This is accomplished via application of successive series of near-identity transformations on the original system. The near-identity coordinate transformations are nonlinear and local. The reader is directed to the well-documented literature on normal form theory found in works by authors such as [32, 33, 34, 35].

To normalize nonlinear systems subjected to external periodic excitation, several authors such as [34, 35, 36, 37] either utilize approaches that introduce equation variables and/or detuning parameters or incorporate a bookkeeping parameter in their methodology. However, the augmenting parameters involved seemingly lack a uniform explicit connection to the terms in the dynamic equations under consideration. Consequently, we shall utilize a relatively more straightforward and intuitive approach that involves augmenting the system states by converting the periodic external excitation into a system state. The state augmentation approach here is intuitive because the augmented states directly emanate from the periodic forcing term(s); hence, they are neither ad hoc nor arbitrary. Moreover, neither a detuning parameter nor a bookkeeping parameter is required.

Strictly speaking the behavior of the attitude motion as shown in Eq. (7) is characterized in terms variation in true anomaly, f. However, the true anomaly similarly varies with time; hence, the reference as implicit time history is preferred. To demonstrate implicit time history behavior, we initially select e = 0.2 and σ = 0.3 . Later on in Section 3.3, we shall analyze the impact of different e - σ pair values on the motion.

Both the motion in original coordinates, Eq. (7) and the corresponding state augmented system, will be scrutinized. Figure 2 shows the implicit time history behavior of the motion in original coordinates.

Figure 2.

Motion history behavior in original coordinates for seven complete orbits. (a) Θ and Θ′, (b) Phase portrait.

Similarly, we shall scrutinize the history behavior of the system in Eq. (7) after augmenting the system states. In accordance with binomial expansion theorem [38], since e cos f < 1 and 1 e cos f 1 , then the magnitude of the terms in the binomial series progressively become smaller. Therefore, the binomial expansion of the term 1 + e cos f 1 can be approximated as 1 e cos f . Eq. (7) becomes

Θ = 1 e cos f 2 e sin f Θ 3 σ sin Θ cos Θ + 2 e sin f E9

Let

p = cos f p = sinf = q q = cos f = p E10

After further substituting the trigonometric product term with its series approximation given in Eq. (13) to the 7th order, the motion in Eq. (9) can be expressed as

Θ = 1 ep 2 eq Θ 3 σ Θ 2 3 Θ 3 + 2 15 Θ 5 4 315 Θ 7 + 2 eq E11

Therefore, the augmented system state-space representation with Θ 1 ' = Θ 2 Θ 2 ' = Θ 1 ' is shown in Eq. (12):

Θ 1 Θ 2 q p = 0 1 0 0 3 σ 0 2 e 0 0 0 0 1 0 0 1 0 Θ 1 Θ 2 q p + 0 2 eq ep 1 Θ 2 + Θ 2 3 σep Θ 1 3 σ 1 + ep 2 3 Θ 1 3 + 2 15 Θ 1 5 4 315 Θ 1 7 0 0 E12

Figure 3 shows the augmented state system history behavior. Allowing for the expected minor discrepancies due to series and binomial expansion approximations, the system state response is comparable to that of the system in original coordinates shown in Figure 2.

Figure 3.

Motion history behavior of the augmented state system for seven complete orbits. (a) Θ and Θ′ response, (b) Θ and Θ′ phase portrait, (c) q(f) and p(f) response.

From Figures 2 and 3aandb, the attitude motion is quasiperiodic as characterized by the absence of closed trajectory attractors in the phase space. The pitch angle librates roughly between 1.5 < Θ < + 1.5 radians , while the pitch angle rate of change varies between 1.0 < Θ < + 1.5 . The orbits generally appear to follow a “heart-shaped” path starting at the origin with two non-closing lobes on either side. Conversely, the augmented states are periodic as characterized by the closed circular limit cycle attractor centered at the origin in Figure 3d. The motion behavior discussed here indicates that Θ and Θ are susceptible to instability, unpredictability, and chaos. Consequently, these aspects of the motion flow will be investigated next. Neither the eigenvalues of the linear periodic term of Eq. (7) nor of the state augmented system Eq. (12) can be used to determine stability. Consequently, we’ll have to construct the Floquet transformation matrix (FTM) to analyze the dynamical system’s stability.

3.1 Stability and chaos

Stability analysis is preceded by computation of the dynamical system’s state transition matrix, Φ t . We first prepare the motion in Eq. (7) for expansion via shifted Chebyshev polynomials of the first kind by normalizing the principal period. This is because shifted Chebyshev polynomials are only valid for the period interval 0 1 .

The pitch angle trigonometric product term in Eq. (7) can be represented as a product of the respective Taylor series expansion of sine and cosine as shown in Eq. (13):

sin Θ cos Θ = Θ 2 3 Θ 3 + 2 15 Θ 5 4 315 Θ 7 + + a k Θ 2 k 1 E13

where k = 1,2,3,4 , also series coefficient a k 0 , as k .

We can ignore the terms of order greater than 7 in Eq. (13) without significant loss of accuracy because the follow-on terms have relatively small successive coefficients that rapidly approach zero. For instance, the 9th order term has the coefficient a 5 = 6.7791 × 10 4 , while the 11th order term’s coefficient is a 6 = 1.98412 × 10 5 . Substituting the expanded trigonometric product in Eq. (7), we obtain

Θ = 1 1 + e cos f 2 e sin f Θ 3 σ Θ 2 3 Θ 3 + 2 15 Θ 5 4 315 Θ 7 + 2 e sin f E14

To normalize the principal period, let f = 2 πζ . Equivalently, ζ 0 1 represents duration within the principal period. Let Z represent the principal period, and then this implies that for a periodic term, A ζ = A ζ + Z . It follows

Θ = d Θ 1 2 π , Θ = 1 4 π 2 d 2 Θ d ζ 2 E15

After substituting Θ and Θ from Eq. (15) into Eq. (14), the obtained state-space representation of the normalized attitude motion is given in Eq. (16), where d Θ / = x 1 and d 2 Θ / d ζ 2 = x 2 . Further, Θ 0 = 0 Θ 0 = 0 constitutes the initial conditions of the represented motion which correspond to pitch librations at position 1 of Figure 1b:

x 1 x 2 = 0 1 12 π 2 σ 1 + e cos 2 πζ 4 πe sin 2 πζ 1 + e cos 2 πζ x 1 x 2 + 12 π 2 σ 1 + e cos 2 πζ 0 2 3 x 1 3 2 15 x 1 5 + 4 315 x 1 7 + 0 8 π 2 e sin 2 πζ 1 + e cos 2 πζ E16

It is clear that Eq. (16) is of the form

x ζ = A ζ x ζ + f ζ x + F ζ E17

3.1.1 Floquet multipliers and exponents

To facilitate computation of STM, FTM, and L-F transformation matrices using Chebyshev polynomials, we utilized the Chebfun software package on MATLAB™ [39]. Summarily, Chebfun applies piecewise Chebyshev polynomial interpolation to construct smooth functions over the interval 1 + 1 . Recall that ζ = f / 2 π and Z are the normalized principal period; hence, the computed FTM = Φ Z is

Φ Z = 0.9462 0.0529 1.9796 0.9462 E18

The computed Floquet multipliers are critical since they lie on the unit circle with values of 0.9462 ± 0.3236 i as shown in Figure 4. Consequently, this reveals a marginally stable system for the chosen e - σ pair. It follows that the corresponding Floquet exponents 0 ± 0.3295 i are purely imaginary. This is consistent with the quasiperiodic system phase portraits that illustrate a librational motion.

Figure 4.

Floquet multiplier location.

The motion is hence stable in the sense of Lyapunov, but the inherent oscillations are disruptively significant to jeopardize nominal execution of the spacecraft mission.

3.1.2 Poincaré map

Figure 5 shows the constructed Poincaré section of the flow. There is a discernible main cluster of points in close proximity to the origin but restricted to the positive side of Θ . Relatively scanty, isolated discrete points occupy the lower bottom half of the plot bound by 1 < Θ < 2 .

Figure 5.

Poincaré section of the attitude motion.

The Poincaré section composition suggests two possible flow behaviors. The groupings suggest a quasiperiodic trajectory. On the other hand, the scanty random points devoid of clustering could be due to transient behavior or chaos. A chaotic motion can briefly dwell on a near periodic trajectory before changing to a disparate trajectory with a period that is a multiple of the preceding motion. Consequently, it is needed to further investigate the presence of chaos in the attitude motion dynamics.

3.1.3 Chaos

We define chaos as a bounded aperiodic steady-state motion behavior that is not in equilibrium and is sensitive to initial conditions. A minuscule divergence in the input rapidly grows to spawn an overwhelming difference in the system response. We begin investigating chaos in the attitude motion by plotting the system implicit time history with a minute divergence to the initial condition of the state, Θ . We set the first initial condition to zero; the divergent second initial condition is obtained by adding ϵ = 10 12 to the zero initial condition. The obtained implicit time history of the two curves reveal the onset and progression of an overwhelming difference in the system response as illustrated in Figure 6. The slightly divergent initial condition results in an overwhelming difference in response that begins in the second half of the 4th orbit and then rapidly grows in the subsequent orbits. The attitude motion is hence chaotic.

Figure 6.

Chaos: attitude motion sensitivity to initial conditions.

3.1.4 Lyapunov exponents

To determine the average rate of divergence between the initially neighboring trajectories defined locally in the state space, we shall scrutinize the dynamic behavior of the motion’s Lyapunov exponents. Lyapunov exponent stability analysis affords a means of quantifiable expression for initial conditions sensitivity and dependence (i.e., chaos), by describing the exponential rate of growth or decay of a perturbation to a trajectory as time progresses [40]. Lyapunov exponent λ is numerically expressed as

λ = lim t 1 t ln δy t δy 0 E19

where δy t is the tiny separating perturbation vector between the trajectories. The value of Lyapunov exponent will distinguish the nature of the trajectory according to the following criteria: (i) λ < 0 : Trajectory is stable and the motion is asymptotically stable. (ii) λ = 0 : Trajectory is neutral and the motion is characterized by some sort of steady-state. (iii) λ > 0 : Trajectory is unstable and chaotic. The Lyapunov exponent behavior for the motion given Eq. (16) is illustrated in Figure 7.

Figure 7.

Dynamics of Lyapunov exponents for the attitude motion.

The computed Lyapunov exponents are equal in magnitude but opposite in sign with increasing periods because the flow in Eq. (16) is nonautonomous Hamiltonian. Hamiltonian systems are conservative. Therefore, the magnitude of λ 1 which measures expansion in one direction is equal to the magnitude of λ 2 which measures contraction in another direction. Since λ 1 > 0 always, then as prescribed by the above distinction for λ , the attitude motion is chaotic. This outcome is consistent with the preceding chaos analyses that scrutinized motion sensitivity to initial conditions.

3.1.5 Stability charts

The orbit eccentricity and spacecraft’s ratios of principal moments of inertia are, respectively, defined as e 0 1 and σ 0 1 . On the other hand, our stability analysis so far has been illustrated in the courtesy of arbitrarily set values of e = 0.2 σ = 0.3 . Consequently is it essential to holistically scrutinize the motion behavior for all possible values of e and σ . Constructing stability charts which partition the e - σ plane into stable and unstable regions enables scrutiny of motion stability as e and σ vary simultaneously. Transition curves in stability charts constitute frontiers that separate stable regions from unstable regions. We can derive transition curves in closed form via the FTM. Floquet theory requires a stable system to have a Floquet multiplier of magnitude, ρ k 1 . It can hence be proven that the transition from stability to instability occurs when both Floquet multipliers are equal to 1 or both are equal to −1 (see Figure 4). Therefore, the transition curves in the e - σ plane where the solution to the linear periodic term of Eq. (17) changes from stable to unstable (or vice versa) are determined by the conditions Trace FTM = ± 2 or ρ k = 1 .

Even though we can only construct transition curves associated with the linear term of the attitude motion equation, the outcome provides insightful perusal into the behavior of the whole equilibrium solution. The complete solution behavior can be arrived at by augmenting the evaluated linear periodic stability behavior with the combined effect of the nonlinear term f x ζ and forcing term F ζ in Eq. (17). For instance, if a given e - σ pair is initially located in the unstable region of the FTM-dependent stability chart, then the nonlinear and forcing terms will tend to exacerbate this instability, rendering the complete solution unstable for that particular e σ combination. A similar argument can be made for an e - σ pair located in a stable region. The constructed stability chart of the attitude motion is shown in Figure 8.

Figure 8.

Stability chart of the attitude motion.

The darker regions constitute stable points while the lighter regions are unstable. The stability chart has a slightly larger stable region than the unstable region. Unstable solutions appear to be dominated by two regions approximately defined by (i) σ 0 0.2 and increasing values of e plus (ii) e 0 0.3 and increasing values of σ . Essentially, a spacecraft with mass distribution such that σ > 0.2 is amenable to a wider range of eccentricity values above 0.2 to achieve an intended stable pitch angle motion.

For instance, we have previously considered the pair e = 0.2 σ = 0.3 ; this pair is located in the stable region corresponding to critical Floquet multipliers confirmed in Figure 4. The motion at this location is stable in the sense of Lyapunov. The stability chart further accords a means of scrutinizing generalized stability behavior trends or commonalities between disparate e - σ pairs. By picking representative e - σ pairs from different regions, we tabulate the resultant illustrative Floquet multipliers as shown in Table 1.

Stability region e σ Floquet multipliers ( ρ 1 , ρ 2 ) ρ 1 ρ 2
0.82 0.6 1.06522 , 0.93877 1.06 0.94
Unstable 0.6 0.1 0.20126 , 4.96854 0.20 4.96
0.8 0.14 0.10212 , 9.79244 0.10 9.79
Marginal 0.2 0.5 0.13463 ± 0.99099 i 1 1
0.4 0.4 0.79713 ± 0.60382 i 1 1
Stable 0.6 0.5 0.1117 ± 0.99373 i 1 1
0.8 0.9 0.71565 ± 0.69844 i 1 1

Table 1.

Stability of representative e - σ pairs.

From Table 1, we note that both selected marginal and stable regions are associated with critical Floquet multipliers. This implies that the pertinent e - σ pair characterizes a motion stable in the sense of Lyapunov; i.e., the pitch angle librations are bound by Θ π + π . This is consistent with the dynamics presented in Figures 1 and 3. However, in the unstable regions, Floquet multipliers have magnitudes ρ k > 1 , implying that the pitch angle wanders beyond ± π .

3.2 Resonance

The attitude dynamics are dominated by the linear and forcing terms delineated in Eq. (17). This is because if we consider the motion composed of only the linear and forcing terms, i.e., x ζ = A ζ x ζ + F ζ , the numerical solution is unbounded as shown in Figure 9. This is not the case if we consider motion composed of any of the following term combinations: x ζ = A ζ x ζ , x ζ = A ζ x ζ + f x ζ , or x ζ = f x ζ + F ζ .

Figure 9.

Resonance in linear and forcing terms.

Moreover, if we consider the L-F transformed motion, Floquet exponents (eigenvalues of R ) represent frequencies associated with the linear term. Periodic elements of the nonlinear matrix are a product of the truncated Fourier series matrices Q ζ and Q 1 ζ . On the other hand, periodic elements of the forcing matrix are likewise multiplied by Q 1 ζ . Subsequently, resonance between the Floquet exponents and any of the periodic terms in the forcing matrix elements will trigger instability in the motion as well.

Bifurcation triggers the system’s equilibrium solutions to transition between the disparate regions of the stability chart. The orbit eccentricity, e, is the bifurcation parameter for the attitude motion (see Eq. (16)). This is because the general on-orbit spacecraft mass properties represented by σ typically tend to be constant for a gravity-gradient stabilized spacecraft. Therefore, it is essential to analyze the equilibrium solution dynamics as small increments are applied on the bifurcation parameter. We develop the normal form of our dynamics in the next section to facilitate bifurcation behavior analysis. Normal forms are not unique, consequently near identity transformation for the state-augmented system, and the L-F transformed system will be undertaken separately.

3.3 Versal deformation of the normal form and bifurcation analysis

Versal deformation refers to embedding the system in a parameterized family of systems containing all possible dynamics that can occur near the bifurcation point. Moreover, the family of systems should be transverse to the bifurcation surface with the number of parameters equal to the codimension of the bifurcation [41]. The attitude motion undergoes a codimension-one bifurcation because only one parameter e is responsible for the loss of stability (for gravity-gradient stabilization to be maintained, σ has to remain fixed). Because our critical Floquet multipliers are complex and lie on the unit circle, this system will experience a Hopf bifurcation. Further, it is “well-known” that Hopf bifurcation is a codimension-one bifurcation. Firstly, we shall formulate the normal forms of nonlinearities up to the cubic order in Eqs. (16) and (12) to demonstrate the intended approach. Normalization of dynamics with higher order nonlinearities can be accomplished via the same techniques.

3.3.1 State-augmented system

To normalize the state-augmented system, we first apply the modal transformation Θ = My to Eq. (12) and obtain Eq. (20), where Θ = Θ 1 Θ 2 q p T :

y = Jy + M 1 0 2 eq ep 1 Θ 2 + Θ 2 3 σep Θ 1 3 σ 1 + ep 2 3 Θ 1 3 0 0 E20

This system now possesses 4th order nonlinearity. J is in the Jordan canonical form. The normal form is evaluated by successive application of near identity transformation of the form

y = v + h 4 v f E21

where h r v is an n × 1 homogeneous vector of monomials in v of degree r. The state-augmented system is independent of periodic coefficients; hence, we solely obtain the TINF given in Eq. (22):

v 1 v 2 v 3 v 4 = iv 1 iv 2 i 0.948683 v 3 + i 30.4002 v 1 v 2 v 3 + i 1.05409 v 3 2 v 4 i 0.948683 v 4 i 30.4002 v 1 v 2 v 4 i 1.05409 v 3 v 4 2 E22

When the external forcing term is augmented as a system state, the magnitude of the external forcing frequency appears as solitary, linear, imaginary conjugate coefficients in the normal form, i.e., 1 (see Eq. (9)). Eigenvalues of the linear matrix in Eq. (12) constitute the conjugate coefficients in the linear terms of the reduced nonlinearity normal forms (i.e., ± i 0.948683 ).

Moreover, after obtaining the straightforward solutions of v 1 f and v 2 f and then substituting these values in the equations for v 3 f and v 4 f , the last two equations of the normal form can be expressed as

v 3 = i 0.948683 + 30.4002 C 1 C 2 v 3 + i 1.05409 v 3 2 v 4 v 4 = i 0.948683 30.4002 C 1 C 2 v 4 i 1.05409 v 3 v 4 2

where C1 and C2 are the integration constants originating from the analytical solutions of v 1 f and v 2 f , respectively.

We shall investigate the bifurcation of Eq. (12) via its normal form given above. Because the periodic system maintains the same general structure, we may treat the respective limit cycles as equilibria and study their bifurcations. We utilize the versal deformation of the normal form to investigate the change in the stability structure of the dynamics in the neighborhood of the critical point of the bifurcation parameter. Essentially, construction of versal deformation of the normal form facilitates characterization of system dynamics at the critical point and its small neighborhood. Therefore, we handily gain complete understanding of the qualitative phase space dynamics of the dynamical system in the neighborhood of the critical point.

We define the normal form versal deformation parameter as μ 1 . The parameter μ 1 represents a small change in the eigenvalues of the normal form corresponding to a small change in the bifurcation parameter in the original system coordinates. It is a prerequisite condition to obtain a relationship between the versal deformation parameter μ 1 and the original system bifurcation parameter e . Incorporating μ 1 in Eq. (22), we obtain Eq. (23):

v 1 v 2 v 3 v 4 = μ 1 i 0 0 0 0 μ 1 + i 0 0 0 0 μ 1 ϖ 0 0 0 0 μ 1 + ϖ v 1 v 2 v 3 v 4 + 0 0 i 1.05409 v 3 2 v 4 i 1.05409 v 3 v 4 2 E23

where ϖ = i 0.948683 30.4002 C 1 C 2 .

By defining small increments on the bifurcation parameter as η , we can write e k = e c + η k to represent the k disparate sets of bifurcation parameter in the neighborhood of the critical parameter e c = 0.2 . We employ the least squares, curve fitting technique proposed by [42] to obtain the relationship between μ 1 and η , as μ 1 = 1.47476 + i 0.301628 η 1.82052 + i 0.414608 η 2 . NB values of C 1 and C 2 in Eq. (23) were evaluated by forward action transformations of initial conditions in the original coordinates, i.e., Θ 1 = Θ 2 = 0 , p = 1 , q = 0 .

The closed-form analytical solutions for v 1 f and v 2 f in the versal deformation in Eq. (23) are straightforward. To obtain v 3 f and v 4 f , we introduce the complex changes of variable v 3 f = u 1 iu 2 and v 4 f = u 1 + iu 2 followed by the polar coordinates u 1 = R cos θ and u 2 = Rsin θ . The last two equations in Eq. (23) become

R = Re μ 1 R θ = 0.948683 30.4002 C 1 C 2 1.05409 R 2 E24

After solving Eq. (24), we utilize the results to complete the closed-form analytical solution of Eq. (23) as given in Eq. (24):

v 1 f = e μ 1 f C 1 exp if v 2 f = e μ 1 f C 2 exp if v 3 f = e μ 1 f C 3 exp 0.948683 30.4 C 1 C 2 f 0.52705 e 2 μ 1 t C 3 2 μ 1 + C 4 i v 4 f = e μ 1 f C 3 exp 0.948683 30.4 C 1 C 2 f 0.52705 e 2 μ 1 f C 3 2 μ 1 + C 4 i E25

Similarly, C 3 and C 4 are the integration constants originating from the analytical solutions of v 3 f and v 4 f , respectively. The values of these integration constants are evaluated from the initial conditions specified in the original coordinates. After back transformation of the normal form closed-form analytical solutions above, we obtain the motion in the original coordinates.

The back transformed v 1 f and v 2 f constitute the augmented states given in Figure 3c. Moreover, μ 1 0 but is generally small in the order of magnitude 10 4 . The integration constants on evaluation are imaginary whose magnitudes are close to identity. Consequently, from Eq. (25), back transformation of the sinusoidal v 1 f and v 2 f will result in the trigonometric augmented states whose amplitude is determined by the magnitude of the integration constants. From Eq. (24), we can express the transient solution of R as R = e i i μ 1 ± C . Since μ 1 0 , the motion is characterized by a locally stable limit cycle in the neighborhood of the bifurcation point. The limit cycle is stable in the sense of Lyapunov but not asymptotically stable. Post-bifurcation attractors that transform into quasiperiodic attractors portraying a limit cycle in the original coordinates are obtained via back transformation and subsequently shown in Figure 10 ( η = 0.00001 ).

Figure 10.

Poincaré map of state augmented motion post-bifurcation behavior.

3.3.2 L-F transformed system

Here, we also demonstrate analysis of bifurcation behavior subject to different values of e - σ pair. Consequently, by utilizing the developed stability chart (Figure 8), we select e = 0.1 σ = 0.2 . This e - σ pair lies on a transition curve; hence, the system is guaranteed to be bifurcating. Again, by considering up to the cubic nonlinearity, the history behavior from Eq. (16) is shown in Figure 11.

Figure 11.

L-F transformed system implicit time history response for seven complete orbits. (a) x1 and x2, and (b) Phase portrait.

The system likewise possesses critical Floquet multipliers that lie on the unit circle of values ( 0.1435 ± 0.9896 i ) and purely imaginary Floquet exponents, ± 1.4268 i . The computed FTM and R matrices are given in Eq. (26):

Φ Z = 0.1435 0.1914 5.1161 0.1435 , R = 0 0.276 7.376 0 E26

The computed periodic L-F transformation Q ζ matrix and Q 1 ζ are plotted in Figure 12.

Figure 12.

Plot of elements of the LFT matrix and its inverse. (a) Qij (ζ), (b) Qij−1 (ζ).

We consider the L-F transformed dynamics in Eq. (27) up to the cubic nonlinearity. After applying the L-F transformation given in Eq. (19) on the attitude motion in Eq. (16), we obtain the system in Eq. (27):

z = Rz + Q 1 0 12 π 2 σ 1 + e cos 2 πζ 2 3 Q 11 z 1 + Q 12 z 2 3 + Q 1 0 8 π 2 e sin 2 πζ 1 + e cos 2 πζ E27

To normalize this externally excited motion, the system states are augmented to convert the system from nonautonomous to autonomous. We define additional states shown in Eq. (28) and plotted in Figure 13:

Figure 13.

Augmented system states behavior. (a) p(ζ) and q(ζ), (b) Phase portrait.

p = cos 2 πζ p = 2 πsin 2 πζ = q q = 4 π 2 cos 2 πζ = 4 π 2 p E28

After substituting the above augmented states into Eq. (27), we obtain the system shown in Eq. (29)—whose order of nonlinearity increases to four. The transformed denominator term 1 + e cos 2 πζ 1 has been approximated by the binomial expansion equivalent, i.e., 1 + ep 1 1 ep .

Apart from raising the order of nonlinearity, state augmentation further introduces periodic linear terms, 4 π eQ 12 1 q and 4 π eQ 22 1 q in Eq. (29). Consequently, the augmented dynamics with linear parameter-variant coefficients necessitate a second L-F transformation to convert the linear term to parameter invariant:

z 1 z 2 p q = 0 0.276 0 0 7.376 0 0 0 0 0 0 1 0 0 4 π 2 0 z 1 z 2 p q + 8 π 2 σ 1 ep Q 11 z 1 + Q 12 z 2 3 Q 12 1 1 ep Q 11 z 1 + Q 12 z 2 3 Q 22 1 0 0 + 4 πe q 1 ep Q 12 1 q 1 ep Q 22 1 0 0 E29

The computed parameters corresponding to the second L-F transformation are as follows: critical Floquet multipliers, 0.1435 ± 0.9897 i , and purely imaginary Floquet exponents, 0 ± 1.4268 i . The computed second FTM and constant R matrices are given in Eqs. (30) and (31), respectively:

Φ Z = 0.1435 0.1914 0.4343 0.054 5.1161 0.1435 2.0844 0.307 0 0 1 1 0 0 0 1 E30
R = 0 0.276 0.6456 0.0022 7.376 0 0.1165 0.4524 0 0 0 0 0 0 0 0 E31

The computed second periodic L-F transformation matrix, Q ζ , and its inverse, Q 1 ζ , are similarly presented as plots in Figure 14.

Figure 14.

Plot of elements of the second LFT matrix and its inverse. (a) Qij* (ζ), (b) Qij*−1 (ζ).

We designate the second L-F transformation as z = Q w . Here, z = z 1 z 2 p q T . After applying the second L-F transformation to the state augmented periodic system, we obtain

w = R w + Q 1 8 π 2 σ 1 ep Q 11 z 1 + Q 12 z 2 3 Q 12 1 1 ep Q 11 z 1 + Q 12 z 2 3 Q 22 1 0 0 Q 1 4 π e 2 qpQ 12 1 qpQ 22 1 0 0 E32

Applying the modal transformation w = M y to Eq. (32) transmutes this system into Eq. (33):

y = Jy + M 1 Q 1 8 π 2 σ 1 ep Q 11 z 1 + Q 12 z 2 3 Q 12 1 1 ep Q 11 z 1 + Q 12 z 2 3 Q 22 1 0 0 M 1 Q 1 4 π e 2 qpQ 12 1 qpQ 22 1 0 0 E33

where

z 1 = i = 1 4 Q 1 i w i , z 2 = i = 1 4 Q 2 i w i , p = i = 1 4 Q 3 i w i , q = i = 1 4 Q 4 i w i

and

w 1 = i = 1 4 M 1 i y i , w 2 = i = 1 4 M 2 i y i , w 3 = i = 1 4 M 3 i y i , w 4 = i = 1 4 M 4 i y i

J , y , and y take the form previously described in Eqs. (20) and (21).

We first evaluate the TDNF and then average out the periodic terms to obtain the simplified TINF. The closed-form analytical solutions for v 1 ζ and v 2 ζ are constants. Variables v 1 ζ and v 2 ζ in the v 3 and v 4 differential equations are substituted by their respective computed constants. This computation is carried out through the forward action transform of the L-F, modal, and near-identity transformations of the initial conditions declared in the original coordinates:

v 1 v 2 v 3 v 4 = 0 0 0.0106945 i 2.12186 v 3 + i 0.0005599 v 3 2 v 4 0.0106945 + i 2.12186 v 4 i 0.0005599 v 3 v 4 2 E34

The Floquet exponents are conjugate coefficients in the linear terms of the normal forms before being multiplied by the substituted constant values equal to v 1 and v 2 .

Similarly, we define the normal form versal deformation parameter as μ 2 and incorporate it into Eq. (34):

v 1 v 2 v 3 v 4 = 0 0 0 0 0 0 0 0 0 0 μ 2 λ 3 0 0 0 0 μ 2 λ 4 v 1 v 2 v 3 v 4 + 0 0 i 0.0005599 v 3 2 v 4 i 0.0005599 v 3 v 4 2 E35

where λ 3 = 0.0106945 + i 2.12186 and λ 4 = 0.0106945 i 2.12186 .

After defining small increments on the bifurcation parameter again as η , we express the k disparate sets of bifurcation parameter in the neighborhood of the critical parameter e c = 0.1 as, e k = e c + η . The relationship between μ 2 and η is evaluated via the procedure previously stated in Section 3.3.1 to yield μ 2 = 0.132784 i 0.842528 η 3.04196 i 7.13545 η 2 .

The closed-form analytical solutions for v 1 ζ and v 2 ζ of the versal deformation normal form are straightforward. To obtain v 3 ζ and v 4 ζ , we introduce the complex changes of variable, v 3 ζ = u 1 iu 2 and v 4 ζ = u 1 + iu 2 followed by the polar coordinates u 1 = R cos θ and u 2 = Rsin θ . The last two equations in (35) become

R = Re μ 2 0.0106945 R θ = 2.12187 + 0.0005599 R 2 E36

We solve Eq. (36) and use the results to obtain the remaining analytical solutions of Eq. (35) as shown in Eq. (37):

v 1 ζ = μ 2 ζ + C 1 v 2 ζ = μ 2 ζ + C 2 v 3 ζ = e 0.01069 + μ 2 ζ C 3 e Γ i v 4 ζ = e 0.01069 + μ 2 ζ C 3 e Γ i E37

where Γ = 2.1219 ζ + 0.00028 e 0.02138 + 2 μ 2 ζ Re C 3 2 μ 2 0.01069 + Re C 4 .

C i ( i = 1,2,3,4 ) are the respective constants of integration whose value is evaluated from the initial conditions specified in the original coordinates. C 1 and C 2 are real, whereas C 3 and C 4 are complex. μ 2 is similarly small in the order of magnitude 10 4 . v 1 ζ and v 2 ζ back transformation via inverse near-identity, modal, and single L-F transformations forms the augmented states given in Eq. (28) and is plotted in Figure 13. Eq. (36) yields a steady-state solution of the limit cycle amplitude as R = Re μ 2 0.0106945 . When μ 2 0 , the solution of v 3 ζ and v 4 ζ results in locally stable limit cycle with amplitude corresponding to R = Re μ 2 0.0106945 . Consequently, the quasiperiodic attractors in the original coordinates delineating a limit cycle are obtained after back transformation as shown in Figure 15 ( η = 0.0001 ).

Figure 15.

Poincaré map of L-F transformed motion post-bifurcation behavior.

Solutions of the versal deformation equations enable investigation of the post-bifurcation steady-state behavior in the small neighborhood of the bifurcation point. However, as observed by [42], this method is only useful for local analysis. This is because minor errors introduced by back transformation close to the bifurcation points significantly grow as you move further away.

4. Attitude motion feedback control

After setting e = 0.2 and σ = 0.3 , the motion in Eqs. (9) and (16) is then numerically integrated to obtain the uncontrolled responses shown in Figure 16a and b, respectively (NB Figure 16a is the same as Figure 2a). The slight difference in the long-term motion behavior in Figure 16b may be attributed to the approximation of the trigonometric product term by a truncated series in Eq. (13) and possibly fidelity of the numerical integrator used in the computation. As established in Section 3.1, the attitude motion for the considered ( e - σ ) pair is quasiperiodic, marginally stable, and chaotic. Despite the system being stable in the sense of Lyapunov, the inherent oscillations are disruptively significant and require stabilization if the spacecraft is expected to successfully conduct its mission.

Figure 16.

Uncontrolled attitude motion.

Motion controller design shall be undertaken on augmented state system, L-F transformed, and near-identity transformed coordinates. We shall hence first transform the system dynamics into these more amenable but topologically equivalent dynamical structures that retain the Lyapunov stability and bifurcation properties of the original system. Augmentation of the attitude dynamic states has been conducted in Section 3. Control law development will be first considered in L-F transformed coordinates and then followed by near-identity transformed coordinates of Eq. (16).

4.1 Lyapunov-Floquet transformation

Prior to computing the L-F transformation matrix Q ζ and its inverse, Q 1 ζ matrices (for the e = 0.2 , σ = 0.3 case), we computed the FTM and R matrices over the interval ζ 0 1 via shifted Chebyshev polynomials of the first kind for the system in Eq. (16). The evaluated aforementioned matrices are shown below. Here we alternatively present elements of Q ζ and Q 1 ζ as truncated Fourier series described in Eq. (8). Previously in Figures 12 and 14, we have presented the periodic plots of these series. Further, recall that ζ = f 2 π and Z are the normalized principal period; hence, FTM = Φ Z :

Φ Z = 0.9462 0.0529 1.9796 0.9462 , R = 0 0.0539 2.0159 0 E38
Q ζ = Q 11 Q 12 Q 21 Q 22 , Q 1 ζ = Q 11 1 Q 12 1 Q 21 1 Q 22 1 E39
where Q 11 = 0.138896 + 1.24968 cos 2 πζ 0.121824 cos 4 πζ + 0.0121464 cos 4 πζ 0.00121854 cos 8 πζ + 0.000122567 cos 10 πζ Q 12 = 0.201902 sin 2 πζ 0.0196812 sin 4 πζ + 0.00196227 sin 6 πt 0.000196857 sin 8 πζ Q 21 = 7.44496 sin 2 πζ + 1.49121 sin 4 πζ 0.224998 sin 6 πζ + 0.0302284 sin 8 πζ Q 22 = 0.0074811 + 1.20128 cos 2 πζ 0.24076 cos 4 πζ + 0.0363337 cos 6 πζ 0.00488192 cos 8 πζ + 0.000615458 cos 10 πζ
Q 11 1 = 0.178684 + 0.831759 cos 2 πζ 0.00416167 cos 4 πζ 0.00201274 cos 6 πζ 0.001088 cos 8 πζ 0.000706808 cos 10 πζ 0.000523275 cos 12 πζ 0.000466034 cos 14 πζ 0.0258801 sin 2 πζ + 0.000259231 sin 4 πζ + 0.000188365 sin 6 πζ + 0.000136072 sin 8 πζ + 0.000110821 sin 10 πζ + 0.000103106 sin 14 πζ
Q 12 1 = 0.0043031 cos 2 πζ 0.000772303 cos 4 πζ + 0.00012908 cos 6 πζ + 0.000115479 cos 8 πζ + 0.000107006 cos 10 πζ + 0.000110097 cos 12 πζ + 0.000112704 cos 14 πζ 0.138297 sin 2 πζ 0.0123985 sin 4 πζ + 0.00137926 sin 6 πζ + 0.000923347 sin 8 πζ + 0.00068248 sin 10 πζ + 0.000583058 sin 12 πζ + 0.000509416 sin 14 πζ
Q 21 1 = 0.155206 cos 2 πζ 0.00365144 cos 4 πζ 0.00363534 cos 6 πζ 0.00338305 cos 8 πζ 0.00343186 cos 10 πζ 0.00327608 cos 12 πζ 0.00335147 cos 14 πζ + 4.98815 sin 2 πζ 0.0586199 sin 4 πζ 0.0388447 sin 6 πζ 0.0270502 sin 8 πζ 0.0218882 sin 10 πζ 0.0173497 sin 12 πζ 0.0151485 sin 14 πζ
Q 22 1 = 0.0835005 + 0.844365 cos 2 πζ + 0.0815938 cos 4 πζ 0.00321618 cos 6 πζ 0.00173261 cos 8 πζ 0.0010434 cos 10 πζ 0.000742632 cos 12 πζ 0.00062635 cos 14 πζ 0.0262723 sin 2 πζ 0.00508249 sin 4 πζ + 0.000300991 sin 6 πζ + 0.00021669 sin 8 πζ + 0.000163595 sin 10 πζ + 0.000140229 sin 12 πζ + 0.000138575 sin 14 πζ

After applying the L-F transformation, x ζ = Q ζ z ζ , to the attitude motion in Eq. (16), it becomes

z = Rz + Q 1 0 12 π 2 σ 1 + e cos 2 πζ 2 3 k 3 2 15 k 5 + 4 315 k 7 + Q 1 0 8 π 2 e sin 2 πζ 1 + e cos 2 πζ E40

where k = Q 11 z 1 + Q 12 z 2 .

The Lyapunov stability properties are preserved in the new coordinates after the system is transformed by the L-F transformation matrix. The L-F transformation theory guarantees that a suitable controller realized in the L-F transformed coordinates will be correspondingly efficacious after back transformation into the original system coordinates. Consequently, we shall endeavor to systematically synthesize suitable controllers to stabilize the motion in the transformed coordinates. Our control synthesis strategy will first consider linear control laws before exploring nonlinear control strategies.

In order to formulate appropriate control laws that would stabilize the quasiperiodic motion analyzed in Section 3, we introduce a control input u t in Eq. (1) as shown below:

Θ ¨ = 3 μ r 3 σ sin Θ cos Θ ω ̇ + u t E41

Using Eqs. (5) and (6), we perform a change of independent variable from time (t) to true anomaly (f). The closed-loop attitude dynamics hence will be

Θ = 1 1 + e cos f 2 e sin f Θ 3 σ sin Θ cos Θ + 2 e sin f + u f E42

The control action u f will generally represent torque per unit moment of inertia as a function of true anomaly. Eq. (42) is first used to synthesize linear control laws followed by nonlinear control law development.

4.2 Linear control

Though linear control law principles are conventionally intended for controlling linear parameter invariant systems [43], we initially consider them to control our nonlinear dynamics as an initial analysis step. Since most linear control methods tend to be relatively simpler to analyze and implement compared to nonlinear control methods, it is prudent to ascertain the suitability of linear control prior to embarking on relatively more complicated techniques. To implement linear control, we shall consider pole-placement approach to determine the negative feedback gain required to stabilize the system.

4.2.1 State-augmented system

The autonomous state-augmented system in Eq. (12) can be represented in abbreviated form as shown in Eq. (43):

Θ = A Θ f + f Θ f E43

where f Θ f constitutes the nonlinear terms, A is the linear matrix, and Θ = Θ 1 Θ 2 q p T is the state vector. To synthesize the parameter-invariant linear state feedback controller, Eq. (42) becomes

Θ = A Θ f + f Θ f + G 1 u f E44

The linear state feedback controller is of the form u = K Θ f , and the control input scaling vector is G 1 = 0 1 0 0 T . Though A is full rank, the linear pair A G 1 is not controllable. This is because the system controllability matrix, CM , shown in Eq. (45) has a rank of 2, instead of rank 4:

CM = 0 1 0 0.9 1 0 0.9 0 0 0 0 0 0 0 0 0 E45

Consequently, a linear state feedback controller cannot stabilize the system dynamics associated with the e - σ pair considered. It is to be noted that the two states, p q in Eq. (44), are virtual states serving to simplify the system but are not accessible in the actual system dynamics. This is illustrated by the fact that controllability matrix does not have full rank.

4.2.2 L-F transformed system

In this case, the parameter-invariant linear state feedback controller is similarly of the form u ζ = Kz ζ . The control input is scaled by the matrix G 2 = 1 1 T in L-F transformed coordinates. Back transformation of the G 2 u ζ product via inverse L-F transformation will guarantee a single control input in the system original coordinates as will be demonstrated in Eqs. (54) and (55). R is full rank and the linear pair R G 2 is controllable. The L-F transformed Eq. (42) will be shown in Eq. (46):

z = Rz ζ + Q 1 f z ζ + Q 1 F ζ + G 2 u ζ E46

Therefore, the system closed-loop dynamics subjected to a linear control law will be of the form

z = R G 2 K z ζ + Q 1 f z ζ + Q 1 F ζ E47

We initially place poles at 1 0 . Then we evaluate the corresponding matrix K = K 1 K 2 to realize this pole placement. Several stable double pole locations with a decreasing factor of 10 (i.e., −1, −10, −100, −1000…) were considered. None of these pole-placement locations demonstrated notable success in stabilizing the system. For instance, poles placed at p 1 = 0.1 p 2 = 0.2 produce a response for a duration of slightly beyond 1.5 cycles before the states abruptly become indeterminate at about ζ 1.74 as shown in Figure 17.

Figure 17.

Linear control of L-F transformed state.

In this analysis, the system response in the original coordinates is realized via the back transformation z ζ = Q 1 ζ x ζ . Therefore, similar to the state augmentation case, the L-F transformed nonlinear system demonstrates inability to be stabilized by a linear control law. The presence of periodic coefficients (elements of Q 1 ζ ) associated with the nonlinear and forcing terms renders the system untenable to be controlled via LTI system control approaches. In general, the “region of application” of linear control for nonlinear systems is dependent on magnitude of nonlinearity and initial conditions. Many times, linear control may stabilize nonlinear systems locally, but this is not guaranteed.

4.3 Nonlinear control

Nonlinear control appears more suitable than linear control to stabilize the attitude motion. However, conventional nonlinear techniques are not readily amenable to dynamics with periodic coefficients and periodic external excitation. Hence, requisite system state augmentation, L-F, or near-identity transformations will be undertaken prior to controller design. We shall first consider sliding mode control (SMC), and then bifurcation control will be implemented on the marginally stable system to stabilize post-bifurcation response.

4.3.1 Sliding mode control

Sliding mode control is a robust nonlinear feedback control methodology that is suitable for achieving accurate tracking for a class of nonlinear systems. SMC methodology is based on variable structure control law that results in the state trajectory “sliding” along a discontinuity surface in the state space [44, 45]. Though SMC is deterministic, nonlinear, and robust, its implementation is prone to undesirable “chattering” along the sliding surface [46]. Design of SMC involves (i) selection of the switching function (stable hyperplane in the state space on which the dynamics will be restricted) and (ii) control law synthesis.

4.3.1.1 State-augmented system

Here, we implement a SMC that tracks a desired null pitch angle via a negative rate of growth. Dynamics in the original coordinates possess periodic coefficients rendering the dynamics unwieldy and unfavorable to synthesize a SMC. Therefore, we develop the SMC law based on the augmented state dynamics—which are liberated from periodic coefficients. To design a sliding mode controller for the state augmented systems, we designate the switching function as given in Eq. (48). The switching function represents the actual system state (i.e., attitude pitch angle) reference error (difference between desired and actual pitch angle) that the controller desires to maintain at zero. Therefore, when s 1 = 0 , Θ 1 0 as Θ 2 0 :

s 1 = β Θ 1 + Θ 2 2 E48

Subsequently, the closed-loop dynamics of the controlled system are similar to Eq. (42) as shown in Eq. (49):

Θ 1 = Θ 2 Θ 2 = 2 eq + 2 e 2 pq 1 Θ 2 + 2 eq Θ 2 3 Θ 1 σ 3 ep Θ 1 σ 3 1 + ep 2 3 Θ 1 3 + 2 15 Θ 1 5 4 315 Θ 1 7 σ + u f q = p p = q E49

where u f represents the control input. Derivative of the switching function with Θ 1 and Θ 2 substituted from Eq. (49) is

s 1 = 2 105 β Θ 1 + Θ 2 1 + ep 210 eq + σ Θ 1 315 + 210 Θ 1 2 42 Θ 1 4 + 4 Θ 1 6 + 105 β 2 e 1 + ep q Θ 2 + u f E50

Setting β = 1 , we derive the following control input:

u f = 2 105 Θ 1 + Θ 2 1 + ep 210 eq + σ Θ 1 315 + 210 Θ 1 2 42 Θ 1 4 + 4 Θ 1 6 105 1 2 e 1 + ep q Θ 2 ρ sgn s 1 E51

where

ρ Θ > 2 105 β Θ 1 + Θ 2 1 + ep 210 eq + σ Θ 1 315 + 210 Θ 1 2 42 Θ 1 4 + 4 Θ 1 6 105 β 2 e 1 + ep q Θ 2 E52

A sigmoid function, s 1 s 1 + ε is preferred instead of the signum function to reduce chattering around the sliding surface typical of sliding mode controllers. ε is generally small. Employing a direct Lyapunov approach, stability of the sliding mode controller applied here is ascertained by setting V = 1 2 s 1 2 to be the Lyapunov function. Hence, V = s 1 s 1 . The guaranteed negative definiteness of the Lyapunov function derivative demonstrated in Eq. (53) points to a stable controller. This equation is obtained after subsequent substitution for s 1 and u f in the equation of V :

V = 2 s 1 3 2 ϱ s 1 s 1 + ϵ < 0 s 1 0 E53

Figure 18 shows the sliding mode controlled system response.

Figure 18.

Sliding mode controlled actual and augmented states. (a) Θ1 and Θ2, (b) q and p.

It is observable that the sliding mode controller in the state-augmented system achieves stabilization of the motion throughout any number of orbits. Both Θ 1 and Θ 2 are adequately confined to zero. The augmented states p and q remain unaffected.

4.3.1.2 L-F transformed system

In this case, we similarly assume a control input u ζ first applied to Eq. (16), prior to L-F transformation as shown in Eq. (54). G 3 = 0 1 T is the control input scaling vector:

x 1 x 2 = 0 1 12 π 2 σ 1 + e cos 2 πζ 4 πe sin 2 πζ 1 + e cos 2 πζ x 1 x 2 + 12 π 2 σ 1 + e cos 2 πζ 0 2 3 x 1 3 2 15 x 1 5 + 4 315 x 1 7 + 0 8 π 2 e sin 2 πζ 1 + e cos 2 πζ + G 3 u ζ E54

It then follows from Eq. (40) that the controlled L-F transformed system is as shown in Eq. (55):

z = R 11 z 1 + R 12 z 2 + Q 12 1 12 π 2 σ 1 + e cos 2 πζ 2 3 k 3 2 15 k 5 + 4 315 k 7 + 8 π 2 e sin 2 πζ 1 + e cos 2 πζ + u ζ z = R 21 z 1 + R 22 z 2 + Q 22 1 12 π 2 σ 1 + e cos 2 πζ 2 3 k 3 2 15 k 5 + 4 315 k 7 + 8 π 2 e sin 2 πζ 1 + e cos 2 πζ + u ζ E55

where k = Q 11 z 1 + Q 12 z 2 . We define the sliding function according to Eq. (56) to ensure that when s 2 = 0 , z 1 0 as z 2 0 . The sliding surface represents the reference pitch angle error. The controller attempts to maintain a zero error throughout, i.e., s 2 = 0 , f > 0 :

s 2 = α z 1 + z 2 E56

After obtaining the derivative of the sliding function, we substitute for z 1 and z 2 from Eq. (55) to obtain Eq. (57). Moreover, from Eq. (38), R 11 = R 22 = 0 :

s 2 = R 21 z 1 + α R 12 z 2 + α Q 12 1 + Q 22 1 12 π 2 σ 1 + e cos 2 πζ 2 3 k 3 2 15 k 5 + 4 315 k 7 + 8 π 2 e sin 2 πζ 1 + e cos 2 πζ + u ζ E57

From Eq. (57), we set α = 1 and derive the following sliding mode control law:

u ζ = R 21 z 1 R 12 z 2 Q 12 1 + Q 22 1 12 π 2 σ 1 + e cos 2 πζ 2 3 k 3 2 15 k 5 + 4 315 k 7 8 π 2 e sin 2 πζ 1 + e cos 2 πζ 1 Q 12 1 + Q 22 1 ρ sgn s 2 E58

where

ρ z > R 21 z 1 R 12 z 2 Q 12 1 + Q 22 1 12 π 2 σ 1 + e cos 2 πζ 2 3 k 3 2 15 k 5 + 4 315 k 7 8 π 2 e sin 2 πζ 1 + e cos 2 πζ 1 Q 12 1 + Q 22 1 ρ sgn s 2 E59

To reduce chattering around the sliding surface typical of sliding mode controllers due to fast switching of the signum function, a sigmoid function is similarly preferred. We again apply the direct Lyapunov approach to analyze the sliding mode controller stability by selecting V = 1 2 s 2 2 as the Lyapunov function. Asymptotic stability will be guaranteed if the sliding function derivative is negative definite. Hence, the switching function derivative is V = s 2 s 2 . Substituting for s 2 with the control input likewise substituted, we obtain the stability-criteria satisfying relationship below:

V = s 2 ρ s 2 s 2 + ϵ < 0 , s 2 0 E60

Figure 19 shows the sliding mode controlled system in L-F transformed coordinates. The response in Figure 19 is back transformed via the inverse L-F transformation resulting in controlled states in the original coordinates shown in Figure 20.

Figure 19.

Sliding mode controlled states in L-F transformed coordinates.

Figure 20.

Sliding mode controlled states in original coordinates. (a) x1 state, (b) x2 state.

We observe that, similar to the state augmented case, the sliding mode controller stabilizes the L-F transformed motion as well by invariably confining the states to zero as desired. Though specific values for the e - σ were used to demonstrate this technique, stabilization of the planar pitch motion by SMC approach is independent of the assigned e - σ values. However, the possibility of a synthesized sliding mode controller being impractical to implement exists if the required control effort is colossally prohibitive.

4.3.2 Bifurcation control

The critical Floquet multipliers corresponding to purely imaginary Floquet exponents indicate that the system is in the stability boundary. Consequently, it is essential to stabilize the system post bifurcation apart from modifying other motion characteristics such as rate of growth. To achieve these objectives, we engage nonlinear bifurcation control with full state feedback. Synthesis of such a controller is facilitated by the normalized dynamics which are relatively more tractable compared to the dynamics as represented in the original coordinates. Periodic coefficients and complexity in structure of the dynamic equations in the original coordinates drastically convolute synthesis of bifurcation control law. Dynamics of the states in the original coordinates will eventually be obtained via back transformation of the normal form, modal, and L-F transformations. Location of the complex Floquet multipliers on the unit circle (Figure 4) indicates that the pitch attitude motion is undergoing a Hopf bifurcation with a limit cycle attractor of controllable radius. Therefore, the structure of the normal form will also verify a Hopf bifurcation occurring in the neighborhood of the critical point of the bifurcation parameter (i.e., orbit eccentricity).

To illustrate the intended approach, we shall formulate the normal forms of nonlinearities up to the cubic order in Eqs. (12) and (16). Normalization of dynamics with higher order nonlinearities can be accomplished through the same techniques. Similar to the preceding cases, we’ll consider the augmented states and L-F transformed systems separately.

4.3.2.1 State-augmented system

In Section 3.1.1, we demonstrated how to obtain the TINF of the state-augmented system—shown in Eq. (23). Obtaining the closed-form analytical solutions for v 1 f and v 2 f in Eq. (23) is straightforward. On the other hand to evaluate v 3 f and v 4 f , we introduce the complex changes of variable, v 3 = u 1 iu 2 and v 4 = u 1 + iu 2 followed by the polar coordinates u 1 = Rcosθ and u 2 = Rsin θ . The last two equations in Eq. (23) become

R = 0 θ = 0.948683 30.4002 C 1 C 2 1.05409 R 2 E61

where C1 and C2 are the integration constants obtained when solving for v 1 f and v 2 f , respectively.

We solve Eq. (61) and then utilize the results to complete the closed-form analytical solution of Eq. (23). The closed-form solutions of v 1 f , v 2 f , v 3 f , and v 4 f are then back transformed to the original coordinates producing the uncontrolled motion behavior shown in Figure 20. The system response in Figure 21 is a cognate approximation of the originally obtained numerical solution in Figures 3a and 16. Again, a quasiperiodic motion is characterized by non-closed curves and is observed in the corresponding phase portrait. Moreover, a codimension-one Hopf Bifurcation is verified by the normal form structure.

Figure 21.

Uncontrolled dynamics of the normalized state augmented system for seven complete orbits. (a) Θ and Θ′ response, (b) Phase portrait.

To synthesize a bifurcation control law of the normal form, we first add a control input in Eq. (62):

v 1 v 2 v 3 v 4 = i 0 0 0 0 i 0 0 0 0 i 0.948683 0 0 0 0 i 0.948683 v 1 v 2 v 3 v 4 + 0 0 i 30.4002 v 1 v 2 v 3 + i 1.05409 v 3 2 v 4 i 30.4002 v 1 v 2 v 4 i 1.05409 v 3 v 4 2 + G 4 u E62

Let the scaling matrix and control input, respectively, be of the form

G 4 = 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 , u = γ 1 0 0 K 1 v 1 v 2 v 3 + K 2 v 3 2 v 4 K 1 v 1 v 2 v 3 + K 2 v 3 v 2 4 E63

Back transformation of the G 4 u product via inverse normal form and modal transformations will guarantee a single control input in the system original coordinates as demonstrated in Eqs. (54) and (55). The proportional gains are custom tuned to K 1 = 5 and K 2 = 10 . γ 1 = 1 is a scalable parameter meant to suppress strange trajectory behavior according to Poincaré-Bendixson theorem in the system phase space. The resulting response of the bifurcation-controlled augmented state system is shown in Figure 22. The augmented states remain unaffected as previously shown in Figure 18b.

Figure 22.

Dynamics of the bifurcation-controlled state-augmented system. (a) Θ and Θ′ response, (b) Phase portrait.

The libratory amplitude of the quasiperiodic pitch angle motion is tremendously stabilized and confined to a significantly diminished limit cycle attractor as illustrated in Figure 21.

4.3.2.2 L-F transformed system

As already indicated in Section 3.3.2, in addition to synthesizing bifurcation control law via L-F transformed dynamics, we shall also demonstrate analysis of the spacecraft attitude dynamics due to different values of e and σ . Therefore, e = 0.1 and σ = 0.2 is once again considered in this section. L-F transformation analysis of the attitude dynamics associated with these values of e and σ has been comprehensively conducted in Section 3.3.2. Subsequently, the normalized TINF system was obtained in Eq. (34).

In Eq. (34), the closed-form analytical solutions for v 1 ζ and v 2 ζ are constants. Variables v 1 and v 2 in the v 3 and v 4 differential equations are substituted by their respective computed constants. This computation is carried out through the forward action transform of the L-F, modal, and near-identity transformations of the initial conditions declared in the original coordinates.

The Floquet exponents are conjugate coefficients in the linear terms of the normal forms before being multiplied by the substituted constant values equal to v 1 and v 2 .

To obtain v 3 ζ and v 4 ζ , we introduce the complex changes of variable v 3 = u 1 iu 2 and v 4 = u 1 + iu 2 followed by the polar coordinates u 1 = R cos θ and u 2 = Rsin θ . The last two equations in Eq. (34) become

R = 0.0106945 R θ = 2.12186 + 0.0005599 R 2 E64

Results from Eq. (64) which is easier to solve are then used to obtain the closed-form analytical solutions to Eq. (34). Then, v 1 ζ v 2 ζ v 3 ζ v 4 ζ T are then back transformed to the original coordinates, producing the uncontrolled motion shown in Figure 23. The system response in Figure 23 (with nonzero initial conditions) is a cognate approximation of the originally obtained numerical integration solution in Figure 11. The back transformed augmented states are similarly shown in Figure 24 corresponding to Eq. (28) where the amplitude of q ζ is 2 π times that of p.

Figure 23.

Behavior of the normalized L-F transformed system states for seven complete orbits. (a) Original back transformed normalized states, (b) Phase portrait.

Figure 24.

Behavior of back-transformed normalized augmented states for seven complete orbits.

The normal form in Eq. (34) verifies that this is a system undergoing a codimension-one Hopf bifurcation. To synthesize a bifurcation control law, a control input is added to Eq. (34) as shown below:

v 1 v 2 v 3 v 4 = 0 0 0 0 0 0 0 0 0 0 0.0106945 i 2.12186 0 0 0 0 0.0106945 + i 2.12186 v 1 v 2 v 3 v 4 + 0 0 i 0.0005599 v 3 2 v 4 i 0.0005599 v 3 v 4 2 + G 5 u E65

Let the scaling matrix and control input be of the form

G 5 = 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 , u = γ 2 0 0 K 1 v 3 2 v 4 K 2 v 3 v 4 2 E66

The proportional gains are custom-tuned to K 1 = K 2 = 2 and γ 2 = 1 . Figure 25 shows dynamic behavior of the implemented bifurcation control in original coordinates with nonzero initial conditions.

Figure 25.

State response of the bifurcation-controlled L-F transformed system. (a) x1 and x2 response, (b) Phase portrait.

The oscillating amplitude of the quasiperiodic pitch angle motion is tremendously stabilized relative to the initial behavior illustrated in Figure 23. This hence demonstrates successful control of the post-bifurcation attitude dynamics about the spacecraft center of mass.

Bifurcation control is a nonlinear control technique that affects the behavior of the closed-loop system by modifying nonlinearity and post-bifurcation behavior. Therefore, the location of Floquet multipliers (exponents) is generally preserved post-bifurcation control. Figure 26 shows the preserved locations of the Floquet multipliers before and after bifurcation control ( e = 0.1 and σ = 0.2 ). This location of Floquet multipliers is consistent with the limit cycle shown in Figure 25b corresponding to a simply stable system with relatively subdued librations.

Figure 26.

Preserved locations of Floquet multipliers after bifurcation control.

5. Conclusions

In this chapter, we illustrated techniques for analyzing and stabilizing the attitude motion of a gravity-gradient stabilized spacecraft. The motion dynamics are shown to be nonlinear with periodic coefficients and subjected to external periodic excitation. Methodologies employed here utilize state augmentation, Lyapunov-Floquet transformation theory, and normal forms to realize relatively more tractable dynamical systems that are amenable to conventional controller synthesis techniques. Floquet theory was used to investigate system stability. State augmentation facilitated analysis via normal forms by transforming the dynamical system from nonautonomous to an autonomous one.

Outcome from the analysis showed that the attitude motion is quasiperiodic, chaotic, and stable in the sense of Lyapunov for the particular e - σ pairs considered. Subsequently, the motion stability chart that was constructed facilitated prediction of e - σ combination leading to stable or unstable dynamics. The stable regions of the stability curves were found to predict marginal and not asymptotically stable dynamics. However, the emanating librations need to be stabilized for nominal mission operations to be realized. Conversely, the e - σ combinations located in the unstable regions resulted in aperiodic unstable dynamics. The computed Lyapunov exponents indicate that the chaotic dynamics also depend on initial values of Θ Θ pair, not just on the magnitudes of e - σ pairs.

Both outcomes of the twofold versal deformation analyses (disparate values of e - σ pairs considered) indicate establishment of locally stable limit cycles by the quasiperiodic flow post bifurcation. Since the eccentricity varies as 0 < e < 1 , relatively small deviations from the critical point e c of the order 10 4 < η < 10 3 trigger a significant topological change in the structure of the motion flow.

The quasiperiodic, nonlinear, and periodically forced pitch attitude motion is challenging to control. The synthesized linear controller served as starting point for developing more adept control laws. Not surprisingly, the linear controller failed to stabilize the complexly structured nonlinear dynamical system. As stated, in general the “region of application” of linear control for nonlinear systems is dependent on magnitude of nonlinearity and initial conditions. Many times, linear control may stabilize nonlinear systems locally, but this is not guaranteed.

On-orbit perturbations cause disturbing torques that bifurcate the attitude motion; it is hence imperative to stabilize the system attitude dynamics in the small neighborhood of the bifurcation parameter’s critical point. Local nonlinear bifurcation control law implemented on the attitude motion undergoing a Hopf bifurcation was shown to stabilize the attitude motion. The bifurcation controller which modifies the nonlinearity and post-bifurcation behavior further prevents the attitude motion from becoming chaotic because bifurcation is the path to chaos. Implemented in the TINF, the bifurcation control law would subsequently stabilize the secular and periodic attitude perturbations experienced by a spacecraft in elliptical orbit about its nominal operating point.

Sliding mode control law was based on driving both system states to zero on the sliding surface when the sliding surface reference error is equal to zero. The SMC law was similarly shown to be successful by invariably restricting the pitch angle to zero.

Future work will consider torques generated by sources such as magnetism and oblateness of the earth, atmospheric drag, solar radiation pressure, thermal bending, etc. Further, nonlinearities beyond the cubic term in the L-F transformation and TDNF case of near identity transformation would also be analyzed. In addition, physical implementation of the controllers and derivation of TDNF-based control laws requires future scrutiny. As demonstrated, all the control effort inputs are single torques per unit moment of inertia which for instance can be implemented via thrusters. Consequently, sizing and implementation of the control effort are an essential task.

Acknowledgments

Financial support for this research was partially provided by the Interplanetary Initiative of Arizona State University.

Conflict of interest

The authors declare no conflict of interest.

Nomenclature

e

eccentricity

f

true anomaly, radians

r

orbit radius, m

COM

center of mass

I

identity matrix

I x , I y , I z

principal moment of inertia about the roll, pitch, and yaw axes, respectively, kgm2

LVLH

local-vertical/local-horizontal

P

semilatus rectum

Q t

Lyapunov-Floquet transformation matrix (LFT)

TDNF

time-dependent normal forms

TINF

time-independent normal forms

μ

earth gravitational parameter, 3.986004418 × 10 14 m 3 s 2

σ

I x I z I y

ω

spacecraft angular velocity, rad/s

Θ

spacecraft pitch angle, radians

Ψ

spacecraft roll angle, radians

Ω

spacecraft yaw angle, radians

Φ t

state transition matrix (STM)

Φ T

Floquet transition matrix (FTM)

References

  1. 1. Waswa PMB, Redkar S. A survey of space mission architecture and system actualisation methodologies. International Journal of Space Science and Engineering. 2017;4(3):234-252
  2. 2. Larson WJ, Wertz JR, editors. Space Missions Analysis and Design. El Segundo, CA: Microcosm Press; 2006
  3. 3. Wie B. Space Vehicle Dynamics and Control. Reston, VA: American Institute of Aeronautics and Astronautics; 2008
  4. 4. Suleman A. Multibody dynamics and nonlinear control of flexible space structures. Modal Analysis. 2004;10(11):1639-1661. Available from: https: //doi.org/10.1177/1077546304042049
  5. 5. Beletskii VV. Librations on an eccentric orbit: NASA Technical Translation F-8504. Washington, DC, USA: National Aeronautics and Space Administration; 1963
  6. 6. Fujii HA, Ichiki W. Nonlinear dynamics of the tethered subsatellite system in the station keeping phase. Journal of Guidance, Control, and Dynamics. 2011;20(2):403-406
  7. 7. Beletskii VV. Motion of an artificial satellite about its center of mass: NASA Technical Translation F-429. U.S. Department Commerce: Springfield, VA, USA; 1963
  8. 8. Sherrill RE, Sinclair AJ, Sinha SC, Lovell TA. Lyapunov-Floquet control of satellite relative motion in elliptic orbits. IEEE Transactions on Aerospace and Electronic Systems. 2015;51(4):2800-2810
  9. 9. Mingori DL. Effects of energy dissipation on the attitude stability of dual-spinsatellites. AIAA Journal. 1969;7(1):20-27
  10. 10. Kojima H, Iwasaki M, Fujii HA, Blanksby C, Trivailo P. Nonlinear control of librational motion of tethered satellites in elliptic orbits. Journal of Guidance, Control, and Dynamics. 2004;27(2):229-239. American Institute of Aeronautics and Astronautics
  11. 11. Sinha S, Wu D, Juneja V, Joseph P. Analysis of dynamic systems with periodically varying parameters via Chebyshev polynomials. Journal of Vibration and Acoustics. 1993;115(1):96-102
  12. 12. Vincent TL, Grantham WJ. Nonlinear and Optimal Control Systems. New York: John Wiley & Sons; 1997
  13. 13. Yakubovich VA, Starzhinskii VM. Linear Differential Equations with Periodic Coefficients. Vol. I and II. New York: John Wiley & Sons; 1975
  14. 14. Johnson W. Perturbation solutions for the influence of forward flight on helicopter rotor flapping stability. NASA TM X-62; 1974. p. 165
  15. 15. Jordan DW, Smith P. Nonlinear Ordinary Differential Equations: An Introduction for Scientists and Engineers. 4th ed. Oxford: Oxford University Press; 2007
  16. 16. Lovera M, Marchi ED, Bittanti S. Periodic attitude control techniques for small satellites with magnetic actuators. IEEE Transactions on Control Systems Technology. 2002;10(1):90-95
  17. 17. Rossa FD, Dercole F, Lovera M. Attitude stability analysis for an earth pointing magnetically controlled spacecraft. IFAC Proceedings Volumes. 2013;46(19):518-523. 19th IFAC Symposium on Automatic Control in Aerospace
  18. 18. Sheng J, Elbeyli O, Sun JQ. Stability and optimal feedback controls for time-delayed linear periodic systems. American Institute of Aeronautics and Astronautics Journal. 2004;42(5):908-911
  19. 19. Wisniewski R, Blanke M. Fully magnetic attitude control for spacecraft subject to gravity gradient. Automatica. 1999;35(7):1201-1214
  20. 20. Arduini C, Baiocco P. Active magnetic damping attitude control for gravity gradient stabilized spacecraft. Journal of Guidance, Control, and Dynamics. 1997;20(1):117-122. American Institute of Aeronautics and Astronautics
  21. 21. Chernous’ko FL. On the motion of a satellite about its center of mass under the action of gravitational moments. Journal of Applied Mathematics and Mechanics. 1963;27(3):708-722
  22. 22. Kane T. Attitude stability of earth-pointing satellites. AIAA Journal. 1965 Apr;3:726-731
  23. 23. Breakwell J, Pringle Jr R. Nonlinear resonance affecting gravity-gradient stability. In: 16th International Astronautical Congress; Astrodynamics. Athens, Greece; 1965. pp. 305-325
  24. 24. D'Angelo H. Linear Time-Varying Systems: Analysis and Synthesis. Boston: Allyn and Bacon; 1970
  25. 25. Sinha S, Joseph P. Control of general dynamic systems with periodically varying parameters via Lyapunov-Floquet transformation. Transanctions-American Society of Mechanical Engineers Journal of Dynamic Systems Measurement and Control. 1994;116:650-650
  26. 26. Sherrill RE, Sinclair AJ, Sinha SC, Lovell TA. Time-varying transformations for Hill-Clohessy-Wiltshire solutions in elliptic orbits. Celestial Mechanics and Dynamical Astronomy. 2014;119(1):55-73. DOI: 10.1007/s10569-014-9543-x
  27. 27. Fleig A Jr. On the libration of a gravity gradient stabilized spacecraft in an eccentric orbit. NASA TECHNICAL REPORT R-350; Washington, DC, USA: National Aeronautics and Space Administration; 1970
  28. 28. Coddington A, Levinson N. Theory of Ordinary Differential Equations. International Series in Pure and Applied Mathematics. New York, NY: McGraw-Hill; 1955
  29. 29. Floquet G. Sur les équations différentielles linéaires á coefficients périodiques. Annales scientifiques de l'École Normale Supérieure. 1883;12:47-88. Available from: http://eudml.org/doc/80895
  30. 30. Sinha S. On the analysis of time-periodic nonlinear dynamical systems. Sadhana. 1997;22:411-434
  31. 31. Mason JC, Handscomb DC. Chebyshev Polynomials. Boca Raton, Florida: CRC Press; 2003
  32. 32. Poincaré H. Les méthodes Nouvelles de la mécanique céleste. Vol. 1. Paris: Gauthier-Villars; 1892
  33. 33. Arnold VI, Vogtmann K, Weinstein A. Mathematical Methods of Classical Mechanics. Graduate Texts in Mathematics. New York: Springer; 2013. Available from: https://books.google.com/books?id=UOQlBQAAQBAJ
  34. 34. Nayfeh AH. Method of Normal Forms. Wiley Series in Nonlinear Science. New York: Wiley; 1993
  35. 35. Jezequel L, Lamarque CH. Analysis of non-linear dynamical systems by the normal form theory. Journal of Sound and Vibration. 1991;149(3):429-459
  36. 36. Gabale AP, Sinha SC. A direct analysis of nonlinear systems with external periodic excitations via normal forms. Nonlinear Dynamics. 2009;55(1-2):79-93
  37. 37. Neild SA, Champneys AR, Wagg DJ, Hill TL, Cammarano A. The use of normal forms for analysing nonlinear mechanical vibrations. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. 2015;373(2051):1-22
  38. 38. Herman EA. Convergence, Approximation, and Differential Equations. Hoboken, NJ: Wiley; 1986
  39. 39. Driscoll T A NH, Trefethen LN, editors. Chebfun Guide: For Chebfun version 5. www.chebfun.org; 2014
  40. 40. Nayfeh AH, Balachandran B. Applied Nonlinear Dynamics: Analytical, Computational and Experimental Methods. Wiley Series in Nonlinear Science. Wiley & Sons Inc: New York, NY; 1995
  41. 41. Wiggins S. Introduction to Applied Nonlinear Dynamical Systems and Chaos. Texts in Applied Mathematics. New York: Springer; 2003. Available from: https://books.google.com/books?id=RSI4RGdwnU4C
  42. 42. Dávid A, Sinha SC. Versal deformation and local bifurcation analysis of time-periodic nonlinear systems. Nonlinear Dynamics. 2000;21(4):317-336. Available from: https: //doi.org/10.1023/A:1008330023291.
  43. 43. Franklin GF, Powell J, Emami-Naeini A. Feedback Control of Dynamic Systems. 3rd ed. Reading, MA: Addison-Wesley; 1994
  44. 44. Utkin VI. Sliding Modes and their Application in Variable Structure Systems. Rev. from the 1974 ed. Moscow: Mir Publishers; 1978
  45. 45. Slotine JJ, Sastry SS. Tracking control of non-linear systems using sliding surfaces with application to robot manipulators. In: 1983 American Control Conference; 1983. pp. 132-135
  46. 46. Kachroo P, Tomizuka M. Chattering reduction and error convergence in the sliding-mode control of a class of nonlinear systems. IEEE Transactions on Automatic Control. 1996;41(7):1063-1068

Notes

  • These authors contributed equally.

Written By

Peter M.B. Waswa and Sangram Redkar

Submitted: 21 March 2019 Reviewed: 10 June 2019 Published: 15 January 2020