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: March 21st, 2019 Reviewed: June 10th, 2019 Published: January 15th, 2020

DOI: 10.5772/intechopen.87954

Chapter metrics overview

716 Chapter Downloads

View Full Metrics


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.


  • 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 Izpoints along the radius vector, Iypoints along the normal to the orbit plane, and Ixis along the tangent to the orbit in the LVLH frame (Hill frame) as shown in Figure 1a. In addition, the condition Iy>Ix>Izmust 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):


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


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




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


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., ẋt=Atxt,s.tAt=At+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 xt=Qtztreduces the original nonautonomous linear differential system to an autonomous one with the form żt=Rztwhere Ris an n×nconstant matrix. In the LFT matrix, Qtcan 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 Qijtwith truncated Fourier representations as shown in Eq. (8). Qij1thas similar Fourier representation:


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.2and σ=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 ecosf<1and 1ecosf1, then the magnitude of the terms in the binomial series progressively become smaller. Therefore, the binomial expansion of the term 1+ecosf1can be approximated as 1ecosf. Eq. (7) becomes




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


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


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.5radians, 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 01.

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):


where k=1,2,3,4, also series coefficient ak0, 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 a5=6.7791×104, while the 11th order term’s coefficient is a6=1.98412×105.Substituting the expanded trigonometric product in Eq. (7), we obtain


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


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Θ/=x1and d2Θ/dζ2=x2. Further, Θ0=0Θ0=0constitutes the initial conditions of the represented motion which correspond to pitch librations at position 1 of Figure 1b:


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


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 Zare the normalized principal period; hence, the computed FTM=ΦZis


The computed Floquet multipliers are critical since they lie on the unit circle with values of 0.9462±0.3236ias 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.3295iare 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 ϵ=1012to 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


where δytis 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 λ1which measures expansion in one direction is equal to the magnitude of λ2which measures contraction in another direction. Since λ1>0always, 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 e01and σ01. 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 eand σ. Constructing stability charts which partition the e-σplane into stable and unstable regions enables scrutiny of motion stability as eand σ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, ρk1. 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 TraceFTM=±2or ρ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 fxζ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) σ00.2and increasing values of eplus (ii) e00.3and increasing values of σ. Essentially, a spacecraft with mass distribution such that σ>0.2is 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 regioneσFloquet multipliers (ρ1,ρ2)ρ1ρ2

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ζ+fxζ, or xζ=fxζ+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 Q1ζ. On the other hand, periodic elements of the forcing matrix are likewise multiplied by Q1ζ. 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 eis 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 Θ=Myto Eq. (12) and obtain Eq. (20), where Θ=Θ1Θ2qpT:


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


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


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., ±i0.948683).

Moreover, after obtaining the straightforward solutions of v1fand v2fand then substituting these values in the equations for v3fand v4f, the last two equations of the normal form can be expressed as


where C1 and C2 are the integration constants originating from the analytical solutions of v1fand v2f, 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 μ1represents 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 μ1and the original system bifurcation parameter e. Incorporating μ1in Eq. (22), we obtain Eq. (23):


where ϖ=i0.94868330.4002C1C2.

By defining small increments on the bifurcation parameter as η, we can write ek=ec+ηkto represent the kdisparate sets of bifurcation parameter in the neighborhood of the critical parameter ec=0.2. We employ the least squares, curve fitting technique proposed by [42] to obtain the relationship between μ1and η, as μ1=1.47476+i0.301628η1.82052+i0.414608η2. NB values of C1and C2in 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 v1fand v2fin the versal deformation in Eq. (23) are straightforward. To obtain v3fand v4f, we introduce the complex changes of variable v3f=u1iu2and v4f=u1+iu2followed by the polar coordinates u1=Rcosθand u2=Rsinθ. The last two equations in Eq. (23) become


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


Similarly, C3and C4are the integration constants originating from the analytical solutions of v3fand v4f, 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 v1fand v2fconstitute the augmented states given in Figure 3c. Moreover, μ10but is generally small in the order of magnitude 104. The integration constants on evaluation are imaginary whose magnitudes are close to identity. Consequently, from Eq. (25), back transformation of the sinusoidal v1fand v2fwill 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 Ras R=eiiμ1±C. Since μ10, 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.9896i) and purely imaginary Floquet exponents, ±1.4268i. The computed FTM and Rmatrices are given in Eq. (26):


The computed periodic L-F transformation Qζmatrix and Q1ζ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):


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.


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+ecos2πζ1has been approximated by the binomial expansion equivalent, i.e., 1+ep11ep.

Apart from raising the order of nonlinearity, state augmentation further introduces periodic linear terms, 4πeQ121qand 4πeQ221qin 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:


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


The computed second periodic L-F transformation matrix, Qζ, and its inverse, Q1ζ, 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=Qw. Here, z=z1z2pqT. After applying the second L-F transformation to the state augmented periodic system, we obtain


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






J, y, and ytake 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 v1ζand v2ζare constants. Variables v1ζand v2ζin the v3and v4differential 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 v1and v2.

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


where λ3=0.0106945+i2.12186andλ4=0.0106945i2.12186.

After defining small increments on the bifurcation parameter again as η, we express the kdisparate sets of bifurcation parameter in the neighborhood of the critical parameter ec=0.1as, ek=ec+η. The relationship between μ2and ηis evaluated via the procedure previously stated in Section 3.3.1 to yield μ2=0.132784i0.842528η3.04196i7.13545η2.

The closed-form analytical solutions for v1ζand v2ζof the versal deformation normal form are straightforward. To obtain v3ζand v4ζ, we introduce the complex changes of variable, v3ζ=u1iu2and v4ζ=u1+iu2followed by the polar coordinates u1=Rcosθand u2=Rsinθ. The last two equations in (35) become


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


where Γ=2.1219ζ+0.00028e0.02138+2μ2ζReC32μ20.01069+ReC4.

Ci(i=1,2,3,4) are the respective constants of integration whose value is evaluated from the initial conditions specified in the original coordinates. C1and C2are real, whereas C3and C4are complex. μ2is similarly small in the order of magnitude 104. v1ζand v2ζ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μ20.0106945. When μ20, the solution of v3ζand v4ζresults in locally stable limit cycle with amplitude corresponding to R=Reμ20.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.2and σ=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, Q1ζmatrices (for the e=0.2,σ=0.3case), we computed the FTM and Rmatrices over the interval ζ01via 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 Q1ζ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 ζ=f2πand Zare the normalized principal period; hence, FTM=ΦZ:


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


where k=Q11z1+Q12z2.

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 utin Eq. (1) as shown below:


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


The control action ufwill 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):


where fΘfconstitutes the nonlinear terms, Ais the linear matrix, and Θ=Θ1Θ2qpTis the state vector. To synthesize the parameter-invariant linear state feedback controller, Eq. (42) becomes


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


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, pqin 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 G2=11Tin L-F transformed coordinates. Back transformation of the G2uζ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). Ris full rank and the linear pair RG2is controllable. The L-F transformed Eq. (42) will be shown in Eq. (46):


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


We initially place poles at 10. Then we evaluate the corresponding matrix K=K1K2to 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 p1=0.1p2=0.2produce a response for a duration of slightly beyond 1.5 cycles before the states abruptly become indeterminate at about ζ1.74as 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ζ=Q1ζ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 Q1ζ) 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. 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 s1=0, Θ10as Θ20:


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


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


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




A sigmoid function, s1s1+ε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=12s12to be the Lyapunov function. Hence, V=s1s1. 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 s1and ufin the equation of V:


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 Θ1and Θ2are adequately confined to zero. The augmented states pand qremain unaffected. 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). G3=01Tis the control input scaling vector:


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


where k=Q11z1+Q12z2. We define the sliding function according to Eq. (56) to ensure that when s2=0, z10as z20. The sliding surface represents the reference pitch angle error. The controller attempts to maintain a zero error throughout, i.e., s2=0,f>0:


After obtaining the derivative of the sliding function, we substitute for z1and z2from Eq. (55) to obtain Eq. (57). Moreover, from Eq. (38), R11=R22=0:


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




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=12s22as the Lyapunov function. Asymptotic stability will be guaranteed if the sliding function derivative is negative definite. Hence, the switching function derivative is V=s2s2. Substituting for s2with the control input likewise substituted, we obtain the stability-criteria satisfying relationship below:


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. 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 v1fand v2fin Eq. (23) is straightforward. On the other hand to evaluate v3fand v4f, we introduce the complex changes of variable, v3=u1iu2and v4=u1+iu2followed by the polar coordinates u1=Rcosθand u2=Rsinθ. The last two equations in Eq. (23) become


where C1 and C2 are the integration constants obtained when solving for v1fand v2f, 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 v1f,v2f,v3f,and v4fare 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):


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


Back transformation of the G4uproduct 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 K1=5and K2=10. γ1=1is 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. 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 eand σ. Therefore, e=0.1and σ=0.2is once again considered in this section. L-F transformation analysis of the attitude dynamics associated with these values of eand σ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 v1ζand v2ζare constants. Variables v1and v2in the v3and v4differential 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 v1and v2.

To obtain v3ζand v4ζ, we introduce the complex changes of variable v3=u1iu2and v4=u1+iu2followed by the polar coordinates u1=Rcosθand u2=Rsinθ. The last two equations in Eq. (34) become


Results from Eq. (64) which is easier to solve are then used to obtain the closed-form analytical solutions to Eq. (34). Then, v1ζv2ζv3ζv4ζTare 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:


Let the scaling matrix and control input be of the form


The proportional gains are custom-tuned to K1=K2=2and γ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 andx2 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.1and σ=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 ecof the order 104<η<103trigger 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.



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.






true anomaly, radians


orbit radius, m


center of mass


identity matrix


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




semilatus rectum


Lyapunov-Floquet transformation matrix (LFT)


time-dependent normal forms


time-independent normal forms


earth gravitational parameter, 3.986004418×1014m3s2




spacecraft angular velocity, rad/s


spacecraft pitch angle, radians


spacecraft roll angle, radians


spacecraft yaw angle, radians


state transition matrix (STM)


Floquet transition matrix (FTM)


  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: //
  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:
  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:
  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; 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:
  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: //
  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


  • These authors contributed equally.

Written By

Peter M.B. Waswa and Sangram Redkar

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