Advances in Cislunar Periodic Solutions via Taylor Polynomial Maps
Mohammed Atallah Simone Servadio








 Advances in Cislunar Periodic Solutions via Taylor Polynomial Maps
     Mohammed AtallahPhD Student, Department of Aerospace Engineering, Iowa State University, IA 50011, USA. email: matallah@iastate.edu 
 and Simone ServadioAssistant Professor, Department of Aerospace Engineering, Iowa State University, IA 50011, USA. email: servadio@iastate.edu

    Received xxxx ; accepted xxxx 
==========================================================================================================================================================================================================================================================================================


§ ABSTRACT 
In this paper, novel approaches are developed to explore the dynamics of motion in periodic orbits near libration points in cislunar space using the Differential Algebra (DA) framework. The Circular Restricted Three-Body Problem (CR3BP) models the motion, with initial states derived numerically via differential correction. Periodic orbit families are computed using the Pseudo-Arclength Continuation (PAC) method and fitted. Two newly developed polynomial regression models (PRMs) express initial states as functions of predefined parameters and are used in the DA framework to evaluate propagated states. The initial states, expressed via PRM, are propagated in the DA framework using the fourth-order Runge-Kutta (RK4) method. The resultant polynomials of both PRM and DA are employed to develop a control law that shows significantly reduced control effort compared to the traditional tracking control law, demonstrating their potential for cislunar space applications, particularly those requiring computationally inexpensive low-energy transfers.




 


§ INTRODUCTION

With the increasing potential for deep-space exploration missions, there is significant interest in cislunar space due to its role in designing low-energy trajectories [koon2000dynamical]. However, this region is recognized as a chaotic system because of its multi-body gravitational environment. These characteristics have drawn attention to the bounded motion represented by the periodic and quasi-periodic orbits near libration points [wilmer2024preliminary]. In the past few decades, numerous studies have investigated motion in cislunar space. The Circular Restricted Three-Body Problem (CR3BP) is one of the simplified mathematical models commonly used to find solutions for bounded trajectories. This model is linearized around a libration point to obtain the trajectory of the periodic orbit, which introduces inaccuracy. This inaccuracy is compensated using a high-order differential correction scheme that obtains a more accurate trajectory [richardson1980halo]. However, even the differential correction scheme cannot provide a high-fidelity solution for the trajectory due to unmodeled dynamics and external disturbances. Therefore, there is a need for a real-time correction scheme that retains the satellite in a periodic orbit by leveraging sensor measurements.


In recent years, several missions have utilized periodic orbits near the libration points in cislunar space. For instance, the first stationkeeping operations around ℒ_1 and ℒ_2 in cislunar space were performed by the ARTEMIS mission [woodard2009artemis]. In light of these advancements, the Lunar Orbiter Platform-Gateway (LOP-G) is one of the largest international cooperative space programs, aiming to assemble a space station around the Moon [merri2018lunar]. Additionally, NASA's journey to Mars will utilize cislunar space to conduct advanced operations [national2016nasa]. Moreover, the next decade will witness over thirty missions being launched in the cislunar region [baker2024comprehensive]. These missions require staging locations to conduct various activities, where the periodic orbits near the libration points are being investigated as potential choices for that role [whitley2016options].


The CR3BP is the most commonly used framework for transfers in cislunar space, where approximate trajectory solutions can be obtained analytically [richardson1980analytic]. In [pritchett2018impulsive], a methodology for conducting low-energy transfers between periodic orbits is developed using CR3BP. In [singh2020low], the Halo orbit (HO) family near ℒ_1 is utilized to design a low-thrust transfer of a small spacecraft to a low-altitude lunar orbit. In [davis2017stationkeeping], stationkeeping in cislunar space is investigated, with the CR3BP being employed to generate the periodic orbit families, while the n-body dynamical model simulates higher-fidelity trajectories. In [van2016tadpole], the Pseudo-Arclength Continuation (PAC) method is developed to compute the members across each periodic orbit family based on the CR3BP. In [wilmer2021cislunar], the benefits of ℒ_1 and ℒ_2 HOs for orbit maintenance are investigated, where the CR3BP is used to model the constellations. In [fay2024investigation], search and rescue operations are investigated, and the response times are compared for rescuer spacecraft located in distant retrograde orbits and ℒ_1/ℒ_2 Lyapunov orbit (LO) families. More accurate trajectory solutions can be obtained using the bicircular restricted four-body problem (BCR4BP), as presented in [negri2020generalizing, oshima2022multiple, wilmer2021lagrangian]; however, these solutions are more computationally expensive and cannot be obtained analytically. Therefore, this study employs the CR3BP to develop a methodology for representing motion in a periodic orbit family leveraging Differential Algebra (DA).


DA is a computationally efficient tool based on Taylor expansion, that can be employed to represent differentiable and continuous dynamic models as high-order polynomials [chao2023handbook, hawkes1999modern]. Several tools are supplied in the DA framework to obtain the derivatives and integrals of the models in low-level computation environments, such as FORTRAN [berz1987differential], and C/C++ [massari2018differential]. In addition, DA has been proven to be a reliable tool for numerical integration of Ordinary Differential Equations (ODE) carried out by an arbitrary integration scheme. Several applications have leveraged DA framework, such as describing beam dynamics [berz1988differential], and high-order nonlinear filtering [valli2014nonlinear, cavenago2018based, servadio2021differential]. The DA tool is fundamentally based on expressing a continuous differentiable function as an infinite series expanded at a predetermined operating point [servadio2022maximum]. For a small deviation of this operating point, the series returns a precise value of the function using a finite number of the terms. This introduces the concept of Truncated Power Series (TPS), which is computationally reliable and can be used for applying arithmetic and calculus operations. In this study, the DA is employed to represent the initial states of periodic orbit families as functions of predetermined parameters, then these TPS are propagated according to ODE of the CR3BP using fourth-order Runge-Kutta scheme.



This paper aims to investigate and analyze the periodic orbit families near libration points in cislunar space within the framework of DA. First, the general translational motion in cislunar space is expressed using the CR3BP. Then, the initial states of an arbitrary periodic orbit in a given family are obtained using the linearized model around the nearest libration point. These approximate states are refined using a high-order differential correction scheme to obtain more accurate ones. Next, the members across the family are computed using the Pseudo-Arclength Continuation (PAC) method. After that, these members are employed to fit a Polynomial Regression Model (PRM) using the Least-Squares Error (LSE) method. The resultant polynomials of the periodic orbit initial states are then propagated to specific times using the fourth-order Runge-Kutta scheme within the DA framework. The first propagation process uses absolute time, while the second process uses normalized time, in which each orbit in the family is propagated for a fractional amount of its period. Finally, numerical simulations are conducted to demonstrate the reliability and accuracy of DA in representing different periodic orbit families near libration points in cislunar space. Additionally, a Proportional-Derivative (PD) control law is developed using the proposed method and compared to the traditional tracking control law to demonstrate the optimality of the proposed approach.


The rest of the paper is organized as follows: Section 2 presents the mathematical model of the CR3BP. Section 3 introduces the basics of DA. Section 4 shows the evaluation of the periodic orbits near ℒ_1 and ℒ_2 and the application of DA in that process. Section 5 presents and discusses the results of the numerical simulations and demonstrates the applicability of the proposed methodology. Section 6 concludes the paper.



§ MATHEMATICAL MODEL OF THE CR3BP

The translational motion in cislunar space can be approximated by an autonomous dynamic model by applying the following assumptions: 

    
  * The Earth and the Moon are treated as mass points.
    
  * The Moon moves in a circular orbit around the Earth.
    
  * The gravity of the Earth and the Moon is the only source of force influencing the motion, while all other perturbations are neglected.


Conventionally, the parameters and states in the CR3BP model are dimensionless, and the motion is expressed in rotating coordinates centered at the barycenter of the Earth-Moon system. The X axis is in the direction of the vector between the Earth and the Moon. Equation (<ref>) presents the mathematical model of the CR3BP according to the aforementioned assumptions and conventions.

    ẍ=2 ẏ+x-(1-μ)(x+μ)r_1^3-μ[x-(1-μ)]r_2^3
       ÿ=-2 ẋ+y-(1-μ) yr_1^3-μ yr_2^3
       z̈=-(1-μ) zr_1^3-μ zr_2^3

Here, x, y, and z denote the components of the dimensionless position vector of the satellite, where x points to the Moon, y is in the direction of the relative motion of the Moon with respect to the Earth, and z completes the set according to the right-hand rule. μ = 0.01215 is the dimensionless mass of the Moon. r_1 and r_2 are the relative distances between the satellite and the Earth, and the satellite and the Moon, respectively.



§ PERIODIC ORBITS NEAR ℒ_1 AND ℒ_2

In the CR3BP model, there are five equilibrium points, known as libration points, where the gravitational forces exerted by the Earth and the Moon on a satellite are balanced. The first two points, ℒ_1 and ℒ_2, have special characteristics due to their symmetry relative to the Moon. There are two common families of periodic orbits near these points: LOs, which exist in two-dimensional space [henon1969numerical], and HOs, which exist in three-dimensional space [breakwell1979halo]. Computing these periodic orbits requires a series of iterative steps due to the chaotic behavior of the CR3BP and the absence of an analytical solution for the model, as follows:


    
  * The mathematical equations are linearized at the libration point, and approximate initial states are computed using this linear model.
    
  * An iterative high-order differential correction scheme is employed to determine the period and refine the initial states of the orbit.
    
  * The computed member of the orbit family is used to generate other members in the family using the PAC method.




 §.§ Linearized Equations of Motion

The detailed steps of linearizing the model are found in [richardson1980analytic, richardson1980halo]. The resultant linear model is derived as follows:

    ẍ = 2 ẏ+(1+2 c_2) x
       ÿ = -2 ẏ-(c_2-1) y 
       z̈ = -c_2 z

where  

    c_n=1γ_L^3[( ± 1)^n μ+(-1)^n (1-μ) γ_L^n+1(1 ∓γ_L)^n+1],   (L_1  or  L_2)

Here, μ_E=G M_E, G is the gravitational constant, M_E is the Earth mass, γ_L=r_E / a, r_E is the Earth mean radius, and a is the astronomical unit.



 §.§ Differential Correction for Computing Initial States

The linearized model in Equation (<ref>) has a closed-form analytical solution, from which an analytical formula for the initial states of periodic orbits can be derived. However, the accuracy of these approximate initial states is insufficient due to the chaotic behavior of the system, and the high non-linearity of the real-time system. Therefore, further correction is required to achieve the desired accuracy. The differential correction method, commonly used for this purpose, implements an iterative algorithm using the nonlinear model to modify the initial states. The differential correction scheme used in this study was first proposed in [richardson1980halo]. The procedure of this scheme is as follows: First, the initial states of the periodic orbit are obtained using the closed-form analytical solution of Equation (<ref>). These states are then propagated using the CR3BP model until they intersect the x-z plane for HOs or the y-axis for LOs. Due to the symmetry of the periodic orbit, the states ẋ and ż must equal zero at the intersection point. Next, the state transition matrix is evaluated at this half-period. Using the states and the state transition matrix at the half-period point, the next iteration of the corrected initial states is computed as follows:

    [[     Δ x_0;      Δẏ_0; Δ T_1 / 2 ]]=-Φ^-1[[ ẋ(x_0, ẏ_0, T_1 / 2); ż(x_0, ẏ_0, T_1 / 2); y(x_0, ẏ_0, T_1 / 2) ]]

where Φ is the matrix of partials that is defined as follows:

    Φ={[      ∂ẋ/∂ x_0       ∂ẋ/∂ẏ_0  ∂ẋ/∂ T_1 / 2;      ∂ż/∂ x_0       ∂ż/∂ẏ_0  ∂ż/∂ T_1 / 2;     ∂ y/∂ x_0      ∂ y/∂ẏ_0 ∂ y/∂ T_1 / 2 ]}_t=T_1 / 2

This process repeats until the desired state error is achieved. 




 §.§ Computing Members across the Orbit Families

The differential correction scheme is used obtain the initial states of a single periodic orbit. It might be used to obtain the other members of the periodic orbit family, however, it has a limited range of these members and it is not the most computationally efficient method for that purpose [servadio2022dynamics]. The PAC method is one of the most reliable methods that is developed to compute a wide range of members across the family using a single predetermined member [van2016tadpole, patel2023low]. It computes the members with a step Δ s, that is predetermined depending on the desired number of members, in the tangent direction to the solution manifold. In this paper, the method is employed for both Lyapunov and HOs, though it can be used for all periodic and quasi-periodic orbit families. 

Assume 𝐱_i is the initial state vector of an arbitrary member that satisfies the constraints F(𝐱_i) = 0 of the family. In case of LO, the constraints are y|_t=T/2 = ẋ|_t=T/2 = 0, while ż|_t=T/2 = 0 is added for HOs. Shifting this member by a step size Δ s yields the next member 𝐱_i+1, which also satisfied the constraints F(𝐱_i+1) = 0. To guarantee that the step size equals Δ s in the tangent direction, an additional constraint is added as follows: 

    G(𝐱_i+1)=[[               F(𝐱_i+1); (𝐱_i+1-𝐱_i)^T Δ𝐱_i-Δ s ]]=0

where G(·) is the augmented constraints, and Δ𝐱_i is the null vector of the Jacobian matrix for 𝐱_i, which is defined as Δ𝐱_i=𝔑(D_ F(𝐱_i)). Here, 𝔑(·) denotes the null vector, and D_(·) is the Jacobian matrix of (·). The new Jacobian matrix of the augmented constraints is defined as follows:

    D_G(𝐱_i+1)=[[ D_F(𝐱_i+1);     Δ𝐱_i^T ]]

In order to obtain 𝐱_i+1 that satisfies the constraints in Equation (<ref>), an initial guess of the solution ^0𝐱_i+1 is selected. Then, Newton's method is employed as follows:

    ^k+1𝐱_i+1=^k𝐱_i+1-[D_G(^k𝐱_i+1)]^-1 G(^k𝐱_i+1)

where ^k𝐱_i+1 is the k^th iteration. This process is repeated iteratively until the constraints in Equation (<ref>) are satisfied within a certain tolerance.

Starting from an arbitrary member in the family, the other members can be computed in both directions ±Δ𝐱. In this study, the PAC method is employed to compute the members across LO families, HO, and NRHO families near ℒ_1 and ℒ_2.



  §.§.§ LO Families

Figure <ref> shows several LOs near ℒ_1 and ℒ_2. A key advantage of this method is its ability to obtain orbits that are close to the Moon.




  §.§.§ HO and NRHO Families
 
Similarly, Figures <ref> and <ref> show numerous HOs and NRHOs near ℒ_1 and ℒ_2. Specifically, Figure <ref> displays the ℒ_1 families, separated by the bold blue orbit, while Figure <ref> shows the ℒ_2 families.







§ OVERVIEW OF DIFFERENTIAL ALGEBRA

The basic idea of treating numbers and implementing various operations on them using computers is to represent these numbers with a finite amount of information. Generally, numbers can be irrational and ideally represented by infinite digits, which makes it impractical for computers to handle their ideal representation. Therefore, only a finite amount of relevant information is extracted. These approximations are known as floating-point numbers. This approximated form of the numbers allows operations on real numbers by transforming these numbers into floating-point numbers and implementing the operations as depicted in Figure <ref>. Here, a and b are real numbers, a̅ and b̅ are floating-point numbers, and  denotes an arbitrary operation.

In a similar manner, the DA technique extends the concept of floating-point numbers to encompass differentiable functions [servadio2020recursive]. According to the Taylor expansion, any differentiable function at a certain point can be represented by an infinite series expanded at that point. This brings an analogy between real numbers and differentiable functions, in which both are represented by an infinite amount of information (i.e., digits of real numbers and series coefficients of differentiable functions). Similar to floating-point numbers, DA extracts a finite number of terms to represent the function in an approximated way that can be handled by computers. This approximation is used to implement various operations on these functions in a computationally efficient manner. Figure <ref> depicts the equivalent DA approximation of functions in computer environments [servadio2020nonlinear]. Here, f and g are differentiable functions, while F and G are finite series that represent these functions in the DA framework, defined by their coefficients.


For any differentiable function 𝐲 = f(𝐱), the DA mapping is represented as follows: 

    𝐲(δ𝐱) = _Nℱ^𝐱̂(δ𝐱)

where 𝐱̂ denotes the operating point of the series expansion, δ𝐱 denotes the deviation of the 𝐱 defined as δ𝐱 = 𝐱 - 𝐱̂, and N is the highest order of the series with nonzero coefficient. This representation can be used to derive a highly accurate approximation of the solution for a dynamical system [servadio2022nonlinear]. For a given dynamic system 𝐱̇ = f(t, 𝐱), the DA framework can be employed to evaluate the states at a certain time t_j, as depicted in Figure <ref>.

Here, ℳ^𝐱̂_t_0 → t_j(δ𝐱) denotes the State Transition Propagation Matrix (STPM) of the states propagated from the initial time t_0 to a given time t_j, expressed as a function of the deviation δ𝐱 with respect to the operating point of the expansion 𝐱̂. This approach allows for the propagation of the neighborhood of a given state 𝐱 to multiple times by propagating the series instead of propagating each point individually [servadio2024likelihood]. This method is computationally efficient in any application that requires computing multiple states at different times. Additionally, the accuracy can be balanced with computation time by tuning the order N of the series. 




§ POLYNOMIAL REGRESSION MODEL

The computed members in each periodic orbit family are used to construct PRMs. Two different approaches are implemented: The first is the global PRM, where the domain of the predefined parameter κ is divided into multiple regions. The mean of each region is selected as the operating point κ̂, and the polynomial for that region is fitted. The deviations of the predefined parameter δκ for the members act as the query points, while the initial states of the members serve as the fitted values at these query points. The second approach is the local PRM, which uses the parameter κ at the designed periodic orbit as the operating point for polynomial fitting. However, this approach only uses the neighbors of the designed member to fit the polynomial, in order to avoid fitting issues.



 §.§ Global Polynomial Regression Model

In this study, the x component of the initial states of the members is used as the parameter κ for both Lyapunov and Halo families. The x component is chosen because it is unique for each member, unlike the y component, which is always zero, and the z component, which is zero in Lyapunov families and not unique in Halo families. Figure <ref> illustrates the concept of dividing the domain into multiple regions and fitting a polynomial in each, where κ_il and κ_iu represent the lower and upper bounds of the ith region of the domain, respectively, while 𝒫^κ̂_i_i is the polynomial of the members in the ith region, expanded at the operating point κ̂_i. For any given point κ = x_0, the deviation is calculated with respect to the operating point of each region to evaluate each polynomial. The red line in Figure <ref> represents the deviation of x_0 with respect to the operating point of an arbitrary region, δκ_m.


For any arbitrary parameter κ = x_0, the initial states of a family member 𝐱 are evaluated as follows:

    𝐱(κ = x_0) = ∑_i = 1^M a_i(κ) 𝒫_i^κ̂_i(δκ_i)

where δκ_i = κ - κ̂_i, as visualized in Figure <ref>, and a_i(κ) is the activation function of 𝒫_i, defined as follows:

    a_i(κ) = 
        1    if κ_il≤κ < κ_iu, 
    
        0    otherwise




 §.§ Local Polynomial Regression Model

The local PRM uses the parameter κ of the designed member as the operating point for the polynomial. This approach yields a single polynomial, which is more computationally efficient than the global method and provides a more precise computation of members that are close to the designed member (i.e., when δκ≪ 1). However, this model becomes less effective for larger values of δκ. The initial states of the member are obtained using this model as follows:

    𝐱(κ) =  𝒫^κ̂_d(δκ)

In some cislunar applications, such as station-keeping, only the neighbors of the designed orbit are required, making this local PRM an optimal choice for these cases.



§ PERIODIC ORBIT REPRESENTATION USING DIFFERENTIAL ALGEBRA

As mentioned earlier, a Polynomial Regression Model (PRM) is developed to represent the initial states of each periodic orbit family as a polynomial function of a parameter κ. Using this model, the initial state vector is expressed as a function of δκ, as follows: 𝐱_0(δκ) = _Nℱ^κ̂(δκ). Here, ℱ denotes the series of the initial states. Instead of propagating the initial states of the periodic orbits individually, it is more efficient to propagate the series ℱ. In this approach, the resultant series represents the propagated state vector as a function of δκ and can be used to obtain the state vector at different values of κ in a less computationally expensive manner. The DA technique is employed to represent the propagated state vector 𝐱_ij at a certain time t_j, starting from an initial state vector 𝐱_i0, as a series expanded at a given parameter vector κ̂ to N order, where each parameter vector κ_i can be mapped to a certain initial state vector 𝐱_i0, as follows:

    𝐱_i0(δκ_i) = _Nℳ^κ̂_t_0 → t_0(δκ_i) 
       𝐱_ij(δκ_i)
            = _Nℳ^κ̂_t_0 → t_j(δκ_i)

where ℳ denotes the STPM. This STPM is obtained by propagating the CR3BP as a function of the PRM of the periodic orbit family using the fourth-order Runge-Kutta scheme. Therefore, the STPM returns a mapped state vector that must be in the subspace of the periodic orbit family for any arbitrary κ is the domain. This mapping is initially performed at discretized times; however, it is used to obtain the states at any randomly selected time by interpolating the states of the surrounding points.



 §.§ Propagation with Respect to Normalized Time

Equation (<ref>) represents the mapping of the propagated states at a given time. In this context, any deviation in κ might lead to a significant deviation in 𝐱 due to the different time periods of the deviated members. This large deviation increases the control effort required to transfer to the deviated member. To minimize the deviation of the propagated states, the STPM can be computed at a certain dimensionless normalized time η, where η = t/T_p. In this case, the period T_p is a function of δκ as follows: 

    T_pi(δκ_i) = _N𝒯^δκ̂(δκ_i)

where _N𝒯^δκ̂(δκ_i) is the map function of the time period. 
In this approach, the number of time steps is fixed, while the sampling time is variable and is determined as a ratio of the time period function, as follows:

    T_s(δκ) = 1N_s_N𝒯^κ̂(δκ)

where N_s denotes the fixed number of time steps per period. In this case, the resultant state vector at any normalized time can be obtained by adjusting the number of time steps (e.g., N_s|_10% = 1/10N_s). In this case, the mapping from time to normalized time is performed as follow:

    η(δκ_i) = t_N𝒯^κ̂(δκ)

The propagated states at a certain η is represented as follows:

    𝐱_i0(δκ_i) = _Nℳ^κ̂_η_0→η_0(δκ_i) 
       𝐱_ij(δκ_i) =  _Nℳ^κ̂_t_0 → t_j(δκ_i, T_s(δκ_i)) = _Nℳ^κ̂_η_0→η_j(δκ_i)


where ℳ denotes the STPM using normalized time. Figure <ref> depicts the difference between mapping to time and normalized time, with the horizontal solid lines representing points at the same normalized time η, while the dashed lines representing points at the same time t. 





§ NUMERICAL SIMULATIONS

The proposed PRMs and the DA representation of periodic orbit families are verified through a series of numerical simulations.



 §.§ Precision of the Derived  Polynomial Regression Models

The purpose of the PRM is to obtain the initial states of any member in the periodic orbit family in a computationally efficient manner. In this numerical simulation, the accuracy of the developed model is assessed by measuring the state error after propagation for multiple orbits.
To evaluate the accuracy of the proposed global PRM, random states are generated at random times and propagated using the RK4 method for a finite number of orbits. Both random points and times are generated using a uniform distribution in a MATLAB environment. For the LOs near ℒ_2, the domain is divided into eight regions, with polynomials of order thirty. Figure <ref> shows the propagation of these random initial points over time. Figure <ref> demonstrates that the states maintained the periodic orbit over three orbits with only insignificant deviations. However, after three orbits, Figure <ref> reveals that most of these states fail to maintain the periodic orbits. The accuracy would vary if the number of regions or the polynomial order were changed.  Figure <ref> shows the Root Mean Square Error (RMSE) of both position and velocity after different numbers of orbital revolutions. In this analysis, initial states are randomly generated using the developed PRM. The initial position vector is randomly selected, and the corresponding velocity vector is then computed using the PRM to satisfy the periodic orbit conditions. The figure indicates that both position and velocity errors follow the same trend and are of similar magnitude for all number of samples. The errors start from significantly small values and settle to an order of magnitude of one. This analysis demonstrates that the developed PRM can compute accurate initial states for periodic orbits, maintaining a significantly small error even after up to three orbits.







 §.§ Periodic Orbit Families Using the Derived Polynomial Regression Models

The developed global PRM is applied to compute the members of the Lyapunov and HO families near the libration points ℒ_1 and ℒ_2. The polynomial approximations of these orbits are obtained, representing the initial states of the periodic orbit as a function of x_0. Then, these polynomials are propagated over time using the RK4 method within a DA framework. To validate the method's accuracy and demonstrate its effectiveness, the polynomials are propagated over various time and normalized time intervals. This approach provides a comprehensive view of the orbits' trajectories and their stability characteristics. The results, showcasing the generated polynomials, are presented in this section, highlighting the efficiency and precision of the global and local PRMs in orbit computation and propagation near the libration points.



 §.§ Lyapunov Orbits Family

In LOs, the initial condition on the Earth-Moon rotating frame is governed by two states: x_0 and ẏ_0, as they are planar orbits. Since the predefined parameter of the PRM is x_0, the state ẏ_0 is expressed as a polynomial. However, after propagating the states to a certain time, each of the four planar states will have a different polynomial expressed as a function of δκ or δ x_0. Figure <ref> shows the states of the LOs near ℒ_1 after propagating the initial states to ten different times. Each solid line represents the locus of points that share the same time. Figure <ref> displays the propagation of the states to ten different normalized times, with each solid line representing the locus of points that share the same normalized time. It is worth noting that normalized time propagation covers the entire domain of the family. Additionally, points with the same normalized time are closer to each other than points with the same time, especially in long-term propagation. This demonstrates the superiority of normalized time propagation over time propagation in proximity operations. Figure <ref> illustrates the variation of the position states with x_0 and η, and also demonstrates how η varies with x_0, with each solid line representing normalized times of the family at a given time. Figure <ref> shows the variation of the velocity states with x_0 and η in a similar manner. Here, x_0 den






 §.§ Halo Orbit Families

In HOs, the initial conditions are governed by three states: x_0, z_0, and ẏ_0, as these are 3-D orbits. The states z_0 and ẏ_0 are defined as functions of δ x_0. After propagation, each of the six states will have a polynomial representation at each time or normalized time. Figure <ref> shows the states of the HOs near ℒ_2 after propagating the initial states to ten different times, with each solid line representing the locus of points that share the same time. Figure <ref> illustrates the propagation of the states to ten different normalized times, with each solid line representing the locus of points that share the same normalized time. Figure <ref> shows the variation of the position states with x_0 and η, while Figure <ref> shows the variation of the velocity states.







 §.§ Comparative Analysis of Global and Local Polynomial Regression Models

The strength of the global PRM lies in its ability to cover the entire domain of the orbit family. However, locating the orbit in the correct region can be computationally expensive for some applications. Conversely, if an application only requires information about the neighborhood of a specific design orbit rather than the entire family, the local PRM becomes the optimal choice, since it uses a single polynomial, making it computationally efficient. Figure <ref> compares the accuracy of the global and local PRMs, represented by the error between the final and initial states after one time period. An arbitrary point is selected for comparison, and its neighboring points are propagated using the three methods. The comparison shows that the global model maintains consistent accuracy across the domain, whereas the local models exhibit better accuracy near the designed operating point. However, as expected, the local model fails to provide accurate results when the states deviate significantly from the designed operating point. It is worth noting that the global PRM requires evaluating multiple polynomials to maintain consistent error, while the local PRM requires evaluating only a single polynomial.




 §.§ Proportional-Derivative Controller

A PD control law is designed to verify the applicability of the proposed models in real-world missions. The controller gains are manually tuned until the system is stabilized. The objective of the controller is to transfer the satellite to one of the nearest periodic orbits while keeping it within the same family. The control law is implemented as follows:

    
  * The position of the satellite is measured to identify the closest periodic orbit, as determined by the proposed model.
    
  * The nearest periodic orbit parameter, κ_0, and its normalized time, η_0, are determined using the normalized time mapping, as depicted in Figure <ref> and Figure <ref>.
    
  * The normalized transfer time, η_t, is selected; then, the reference state is obtained as follows: 
    𝐱_r = _Nℳ^κ̂_η_0→η_0 + η_t(δκ_0)

    
  * The thrust is evaluated as follows: 𝐮 = - [𝐊] ×(𝐱 - 𝐱_r),
    where [𝐊] = Diag(k_p 𝐈_3,  k_d 𝐈_3). Here, 𝐈_3 denotes the identity matrix, and k_p and k_d are the controller gains.
    
  * The equivalent transfer time is computed as follows: t_t = η_t×_N𝒯^κ̂(δκ_0).
    
  * The velocity impulse Δ𝐕 is evaluated as follows: Δ𝐕 = 𝐮  t_t, which is a valid approximation for small t_t.
    
  * These steps are repeated frequently after each transfer.


The control law is tested on the LO family near ℒ_2, where a random periodic orbit is selected, and the satellite starts from a random point on that orbit with a disturbed velocity. The control law is then implemented to retain the satellite within the family by performing multiple transfers. The transfer time is set to η_t = 0.05. Figure <ref> shows the transfers over ten orbits until the satellite is retained in a specific periodic orbit. Figure <ref> illustrates the time history of the velocity impulses, indicating that the satellite rapidly converges to a steady-state periodic orbit within two revolutions. Although the final periodic orbit differs from the initial one, it maintains the same ground track on the Moon. These results emphasize the potential of the proposed method for various cislunar applications, including low-energy transfers.



The proposed control law is compared to a traditional controller, which tracks the states of a predetermined orbit without using the STPM. As previously mentioned, the proposed method converges to a random orbit within the family. To efficiently compare the two methods, the same simulation parameters are used, and the final orbit obtained from the proposed method is set as the target orbit for the traditional method. Figure <ref> presents the control efforts of both methods, demonstrating a significant reduction of 13.54% in the impulses required by the proposed method.




§ CONCLUSIONS

In conclusion, this study explores the dynamics of motion in periodic orbits near libration points in cislunar space using the DA framework. The CR3BP model is employed to describe the motion in this environment. Initial states of the periodic orbits are numerically generated using a differential correction scheme that leverages an analytical solution, while the members of each orbit family are computed using the PAC method. These computed members are then used to fit PRMs for the orbit families, with the initial states expressed as functions of predefined parameters. These regression models are incorporated into the DA framework to evaluate propagated states at a given time as functions of deviations in these predefined parameters. The accuracy of the computed states at various times is assessed, and the execution time for computing these states is compared with traditional propagation methods using the Runge-Kutta method. The analysis demonstrates the effectiveness of using DA for representing motion in periodic orbits in cislunar space. Additionally, it shows significantly reduced control effort compared to the
traditional tracking control law.




§ ACKNOWLEDGEMENT

The authors wish to thank Erin Ashley and Batu Candan for their help in reviewing and improving this paper.

ieeetr


