A Functional Interpolation Approach to Compute Periodic Orbits in the Circular-Restricted Three-Body Problem

Abstract

In this paper, we develop a method to solve for periodic orbits, i.e., Lyapunov and Halo orbits, using a functional interpolation scheme called the Theory of Functional Connections (TFC). Using this technique, a periodic constraint is analytically embedded into the TFC constrained expression. By doing this, the system of differential equations governing the three-body problem is transformed into an unconstrained optimization problem where simple numerical schemes can be used to find a solution, e.g., nonlinear least-squares is used. This allows for a simpler numerical implementation with comparable accuracy and speed to the traditional differential corrector method.

1. Introduction

According to Poincaré, “periodic orbits” provide the only gateway into the otherwise impenetrable domain of nonlinear dynamics. With the advent of space exploration, periodic orbits have become an indispensable part of missions in space. The amazing fish-like Apollo orbit was the first three-body orbit used for space missions. The second three-body orbit used for space missions was the Halo orbit, discovered by Robert Farquhar in his Ph.D. thesis [1] under John Breakwell [2]. In 1978, Farquhar convinced NASA and led the International Sun-Earth Explorer 3 mission (ISEE3) to study the Sun from a Halo orbit around the Earth’s L1 Lagrange point. Farquhar’s original idea was to place a satellite in Halo orbit around the Lunar L2 for telecommunication support for the backside of the Moon. Today, this is indeed part of NASA’s planned return of humans to the Moon in the next few years.

Computing Halo orbits, and finding nonlinear periodic orbits in general, has become an important part of modern mission design. This paper introduces a novel, simple, and efficient method for finding periodic orbits using the Theory of Functional Connections (TFC) [3,4]. First, we review perhaps the most popular standard method for computing periodic orbits: the differential correction method (also called shooting method), such as presented by Kathleen Howell [5]. One begins with an approximate solution obtained typically from normal form expansions. Using the variational equation, the guess solution is iteratively corrected for periodicity. Assuming that the initial guess is in a reasonable basis of attraction to a periodic orbit, the process converges to a periodic orbit. In Hamiltonian systems, periodic orbits typically occur in one-parameter families. Often, there are multiple families nearby. Hence, the convergence may not always lead to the desired orbit. Moreover, control over the specific features of the periodic orbit, such as its period or energy, requires additional work—for example, using continuation methods to reach the exact orbit desired. Using TFC, the functional interpolation theory, a simpler formulation and more efficient algorithm for finding periodic orbits is possible.

With this work, we removed the need for a shooting method by analytically embedding the boundary conditions in the solution space using TFC [3,4]. These expressions, which will be explicitly derived in the following sections, transform a differential equation subject to constraints into an unconstrained problem, reducing the whole solution function space to only the subspace satisfying the constraints. This method drastically simplifies the problem and allows for simple numerical implementation; in our case, we utilize a nonlinear least-squares technique.

2. Circular-Restricted Three-Body Problem

The circular-restricted three-body problem (CR3BP) is a dynamical model used to describe the motion of a particle $\mathit{r}={\{x,y,z\}}^{\mathrm{T}}$ of negligible mass under the influence of a primary body of mass $m_1$ and the secondary body of mass $m_2$. Furthermore, the orbits of $m_1$ and $m_2$ are subject to circular motion about the system’s barycenter and lie in the x-y plane as depicted in Figure 1. Following this, the system can be non-dimensionalized by the following scaled units: unit mass is defined as $m_1+m_2$; unit length is taken as the separation between $m_1$ and $m_2$; the unit time is chosen such that the orbital periods of $m_1$ and $m_2$ about the system’s barycenter are $2\pi$. By following these steps, the system can be reduced to a single parameter called the mass parameter, $\mu$, where,

$$\mu ={\displaystyle \frac{m_2}{m_1+m_2}}$$

(1)

From this, we define the terms ${\mu }_1$ and ${\mu }_2$ based on Equation (1) as

$${\mu }_1=1-\mu \hspace{1em}\mathrm{and}\hspace{1em}{\mu }_2=\mu .$$

Using this definition of the system, the equations of motion can be derived in the rotating frame, leading to the following system of equations:

$$\begin{array}{c}\hfill \ddot{x}-2\dot{y}={\displaystyle \frac{\partial \mathsf{\Omega }}{\partial x}}\\ \hfill \ddot{y}+2\dot{x}={\displaystyle \frac{\partial \mathsf{\Omega }}{\partial y}}\\ \hfill \ddot{z}={\displaystyle \frac{\partial \mathsf{\Omega }}{\partial z}}\end{array}$$

(2)

such that the shorthand dot notation is used to signify a derivative with respect to time (e.g., $\dot{x}:={\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}}$). Additionally, $\mathsf{\Omega }$ is defined as

$$\begin{array}{c}\hfill \mathsf{\Omega }(x,y,z):={\displaystyle \frac{1}{2}}(x^2+y^2)+{\displaystyle \frac{1-\mu }{R_1}}+{\displaystyle \frac{\mu }{R_2}}+{\displaystyle \frac{1}{2}}(1-\mu )\mu \end{array}$$

(3)

where $R_1=\sqrt{{(x+\mu )}^2+y^2+z^2}$ and $R_2=\sqrt{{(x+\mu -1)}^2+y^2+z^2}$ are the distances to the primaries. Furthermore, the equations of motion are Hamiltonian and independent of time and, thanks to Noether’s theorem, have an energy integral of motion $\mathcal{E}$, where, in the celestial mechanics community, the Jacobi constant is used, which is $C:=-2\mathcal{E}$ and given as

$$C=2\mathsf{\Omega }-(\dot{x}+\dot{y}+\dot{z})=(x^2+y^2)+2{\displaystyle \frac{1-\mu }{R_1}}+2{\displaystyle \frac{\mu }{R_2}}+(1-\mu )\mu -(\dot{x}+\dot{y}+\dot{z}).$$

(4)

Moving forward, we solve the dynamics defined by the system of equations in Equation (2) such that the orbit is at a fixed energy level (or rather Jacobi constant) using Equation (4). For our implementation, it is useful to define the residuals of these equations:

$$\begin{array}{c}\hfill 0=F_x:=\ddot{x}-2\dot{y}-{\displaystyle \frac{\partial \mathsf{\Omega }}{\partial x}}\end{array}$$

(5)

$$\begin{array}{c}\hfill 0=F_y:=\ddot{y}+2\dot{x}-{\displaystyle \frac{\partial \mathsf{\Omega }}{\partial y}}\end{array}$$

(6)

$$\begin{array}{c}\hfill 0=F_z:=\ddot{z}-{\displaystyle \frac{\partial \mathsf{\Omega }}{\partial z}}\end{array}$$

(7)

$$\begin{array}{c}\hfill 0=F_c:=(x^2+y^2)+2{\displaystyle \frac{1-\mu }{R_1}}+2{\displaystyle \frac{\mu }{R_2}}+(1-\mu )\mu -(\dot{x}+\dot{y}+\dot{z})-C\end{array}$$

(8)

Next, we generate analytical expressions for the states to guarantee a periodic orbit.

However, before moving forward, it is necessary to summarize the major steps of the differential corrector method [5], since this technique is used as a baseline to compare results with the proposed method. In general, the differential correction method can be classified as a shooting method since it relies on transforming the boundary-value problem into an initial-value problem where the initial conditions are iteratively adjusted. First, unlike the the proposed technique, the differential correction method can take advantage of the symmetry of the orbits about the x-z plane, which leads to the following conditions:

$$X\left(0\right)={\left\{\begin{array}{cccccc}{x}_0,& 0,& z_0,& 0,& {\dot{y}}_0,& 0\end{array}\right\}}^{\mathrm{T}}$$

(9)

and

$$X(T/2)={\left\{\begin{array}{cccccc}x,& 0,& z,& 0,& \dot{y},& 0\end{array}\right\}}^{\mathrm{T}}$$

(10)

where Equations (9) and (10) represent the initial state and the state at half of the period, T, where the orbit crosses the x-z plane (the sign of the y variable). Therefore, the equations of motion only need to be propagated to this sign change. At this point, if the values of $|\dot{x}|$ and $|\dot{z}|$ are zeros within an acceptable tolerance (Ref. [5] uses ${10}^{-8}$), then the orbit is considered periodic. If they are not, the initial conditions are updated by determining

$$\delta X\left(0\right)={\left\{\begin{array}{cccccc}\delta x_0,& 0,& \delta z_0,& 0,& \delta {\dot{y}}_0,& 0\end{array}\right\}}^{\mathrm{T}}$$

(11)

and adding the result to Equation (9). The values in Equation (11) are determined by the state transitions matrix calculated in the propagation step and the fact that the only changes to state should be the variables $\dot{x}$ and $\dot{z}$, which should be equal to the negative of their values at $T/2$, i.e., $\delta \dot{x}$ should be $-\dot{x}$ (see Refs. [2,5] for the detailed equations). After updating the initial conditions, the equations of motion are propagated again and the process continues until the periodic conditions are met.

3. Solving the Problem Using the Theory of Functional Connections Framework

The numerical technique to solve differential equations based on the Theory of Functional Connections [4] has produced machine error accuracy within milliseconds for both linear and nonlinear [6] differential equations for a wide array of constraint cases. In fact, several works have already adopted the TFC framework to solve mathematical problems including hybrid systems [7], optimal control problems [8,9], quadratic and nonlinear programming [10], and other applications [11,12].

The foundation of this approach is the ability to derive a class of functionals, called constrained expressions, which have the characteristic of always analytically satisfying assigned linear constraints. Using these constrained expressions, the differential equation can be transformed into an algebraic expression that can ultimately be solved via a variety of optimization techniques. In general, the TFC methodology provides functional interpolation by embedding a set of k linear constraints. Let $g\left(t\right)$ be a general function in our unconstrained function space, G. We define next the subspace, $\mathsf{\Gamma }\subset G$, of constrained functions, $x(t,g(t\left)\right)$, subject to these k linear constraints written as the following functional:

$$x(t,g\left(t\right))=g\left(t\right)+{\displaystyle \sum _{j=1}^k}{\varphi }_j\left(t\right) {\rho }_j(t,g\left(t\right)),$$

(12)

where $g\left(t\right)$ is a free function, and ${\varphi }_j\left(t\right)$ are switching functions—by definition, a function is equal to 1 when evaluated at the constraint it is referencing, and equal to 0 when evaluated at all other constraints. In our case, the switching functions are composed of a set of mutually linearly independent functions called support functions, $s_k\left(t\right)$, with unknown coefficients ${\alpha }_{ij}$, such that ${\varphi }_j\left(t\right)=s_i\left(t\right) {\alpha }_{ij}$; however, the switching functions can be very general depending on the nature of the original function space, G.

Furthermore, ${\rho }_j(t,g\left(t\right))$ are called projection functionals, since they are derived by setting the constraint function equal to zero and replacing $x\left(t\right)$ with $g\left(t\right)$. The following step-by-step procedure can be used to derive a constrained expression from Equation (12):

1.

Choose k linearly independent support functions, $s_k\left(t\right)$.

2.

Write each switching function as a linear combination of the support functions with k unknown coefficients.

3.

Based on the switching function definition, write a system of equations to solve for the unknown ${\alpha }_{ij}$ coefficients.

3.1. Derivation of Constrained Expression

First, since, in our problem, the orbital period is unknown, let us normalize the time domain from $t\in [0,T]$, where $t_0=0, t_f=T$, to $\tau \in [-1,+1]$, where ${\tau }_0=-1, {\tau }_f=+1$. The simplest mapping is a linear one as follows:

$$\tau \left(t\right)={\tau }_0+{\displaystyle \frac{{\tau }_f-{\tau }_0}{t_f-t_0}}(t-t_0)\hspace{1em}⟷\hspace{1em}t=t_0+{\displaystyle \frac{t_f-t_0}{{\tau }_f-{\tau }_0}}(\tau -{\tau }_0).$$

(13)

where we defined the following coefficient term of the derivative of $\tau$ with respect to t. (Note: here, we define the coefficient term as the square of b in order to ensure that the ratio is always positive, i.e., the problem domain is always consistent. This is important in numerical implementation when solving for the value of b.)

$$b^2:={\displaystyle \frac{\mathrm{d}\tau }{\mathrm{d}t}}={\displaystyle \frac{{\tau }_f-{\tau }_0}{t_f-t_0}},$$

(14)

Now, we write the constrained expression in terms of this new variable $\tau$. Let us define the position vector of the spacecraft as $\mathit{r}\left(\tau \right)={\{r_x\left(\tau \right),r_y\left(\tau \right),r_z\left(\tau \right)\}}^{\mathrm{T}}={\{x\left(\tau \right),y\left(\tau \right),z\left(\tau \right)\}}^{\mathrm{T}}$. Since the problem presented in this paper is to find periodic orbits in the CR3BP, we can use the TFC expression to satisfy the following constraints:

$$\begin{array}{c}\hfill r_i\left({\tau }_0\right)=r_i\left({\tau }_f\right)={\alpha }_i\hspace{1em}\mathrm{and}\hspace{1em}{r}_i^{\prime }\left({\tau }_0\right)=r_i^{\prime }\left({\tau }_f\right)={\displaystyle \frac{{\beta }_i}{b^2}}\end{array}$$

or, in other words, the trajectory must return to the initial state at some period. Since the constraints for each component of the vector $r_i\left(\tau \right)$ are of the same form, Equation (12) can be easily adapted to the vector form:

$$r_i(\tau ,g_i\left(\tau \right))=g_i\left(\tau \right)+{\displaystyle \sum _{j=1}^k}{\varphi }_j\left(\tau \right) {\rho }_j(t,g_i\left(\tau \right)),$$

(15)

Following the process detailed in the prior section, one must first choose the support functions $s_k\left(\tau \right)$. For these specific constraints, the support functions are defined in the same manner as Ref. [9], where $s_k\left(\tau \right)={\tau }^{(k-1)}$. Now, the definition of the switching functions is used to produce a set of equations. For example, the first switching function has the four following equations:

$${\varphi }_1\left({\tau }_0\right)=1,\hspace{1em}{\varphi }_1\left({\tau }_f\right)=0,\hspace{1em}{\varphi }_1^{\prime }\left({\tau }_0\right)=0,\hspace{1em}\mathrm{and}\hspace{1em}{\varphi }_1^{\prime }\left({\tau }_f\right)=0,$$

These equations are expanded in terms of the support functions:

$$\begin{array}{cccccccccc}\hfill {\varphi }_1\left({\tau }_0\right)& =\left(1\right)\hfill & \hfill ·& {\alpha }_{11}+\left({\tau }_0\right)\hfill & \hfill ·& {\alpha }_{21}+\left({\tau }_0^2\right)\hfill & \hfill ·& {\alpha }_{31}+\left({\tau }_0^3\right)\hfill & \hfill ·& {\alpha }_{41}=1\hfill \\ \hfill {\varphi }_1\left({\tau }_f\right)& =\left(1\right)\hfill & \hfill ·& {\alpha }_{11}+\left({\tau }_f\right)\hfill & \hfill ·& {\alpha }_{21}+\left({\tau }_f^2\right)\hfill & \hfill ·& {\alpha }_{31}+\left({\tau }_f^3\right)\hfill & \hfill ·& {\alpha }_{41}=0\hfill \\ \hfill {\varphi }_1^{\prime }\left({\tau }_0\right)& =\left(0\right)\hfill & \hfill ·& {\alpha }_{11}+\left(1\right)\hfill & \hfill ·& {\alpha }_{21}+\left(2{\tau }_0\right)\hfill & \hfill ·& {\alpha }_{31}+\left(3{\tau }_0^2\right)\hfill & \hfill ·& {\alpha }_{41}=0\hfill \\ \hfill {\varphi }_1^{\prime }\left({\tau }_f\right)& =\left(0\right)\hfill & \hfill ·& {\alpha }_{11}+\left(1\right)\hfill & \hfill ·& {\alpha }_{21}+\left(2{\tau }_f\right)\hfill & \hfill ·& {\alpha }_{31}+\left(3{\tau }_f^2\right)\hfill & \hfill ·& {\alpha }_{41}=0\hfill \end{array}$$

which can be written in matrix form as

$$\left[\begin{array}{cccc}1& {\tau }_0& {\tau }_0^2& {\tau }_0^3\\ 1& {\tau }_f& {\tau }_f^2& {\tau }_f^3\\ 0& 1& 2{\tau }_0& 3{\tau }_0^2\\ 0& 1& 2{\tau }_f& 3{\tau }_f^2\end{array}\right]\left\{\begin{array}{c}{\alpha }_{11}\\ {\alpha }_{21}\\ {\alpha }_{31}\\ {\alpha }_{41}\end{array}\right\}=\left\{\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right\}.$$

The same is done for the other three switching functions to produce a set of equations that can be solved by matrix inversion, where ${\tau }_0=-1$ and ${\tau }_f=1$.

$$\left[\begin{array}{cccc}1& {\tau }_0& {\tau }_0^2& {\tau }_0^3\\ 1& {\tau }_f& {\tau }_f^2& {\tau }_f^3\\ 0& 1& 2{\tau }_0& 3{\tau }_0^2\\ 0& 1& 2{\tau }_f& 3{\tau }_f^2\end{array}\right]\left[\begin{array}{cccc}{\alpha }_{11}& {\alpha }_{12}& {\alpha }_{13}& {\alpha }_{14}\\ {\alpha }_{21}& {\alpha }_{22}& {\alpha }_{23}& {\alpha }_{24}\\ {\alpha }_{31}& {\alpha }_{32}& {\alpha }_{33}& {\alpha }_{34}\\ {\alpha }_{41}& {\alpha }_{42}& {\alpha }_{43}& {\alpha }_{44}\end{array}\right]=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right].$$

that is,

$$\left[\begin{array}{cccc}{\alpha }_{11}& {\alpha }_{12}& {\alpha }_{13}& {\alpha }_{14}\\ {\alpha }_{21}& {\alpha }_{22}& {\alpha }_{23}& {\alpha }_{24}\\ {\alpha }_{31}& {\alpha }_{32}& {\alpha }_{33}& {\alpha }_{34}\\ {\alpha }_{41}& {\alpha }_{42}& {\alpha }_{43}& {\alpha }_{44}\end{array}\right]={\left[\begin{array}{cccc}1& {\tau }_0& {\tau }_0^2& {\tau }_0^3\\ 1& {\tau }_f& {\tau }_f^2& {\tau }_f^3\\ 0& 1& 2{\tau }_0& 3{\tau }_0^2\\ 0& 1& 2{\tau }_f& 3{\tau }_f^2\end{array}\right]}^{-1}={\displaystyle \frac{1}{4}}\left[\begin{array}{cccc}2& 2& 1& -1\\ -3& 3& -1& -1\\ 0& 0& -1& 1\\ 1& -1& 1& 1\end{array}\right].$$

This leads to the switching functions,

$$\begin{array}{ccc}{\varphi }_1\left(\tau \right)={\displaystyle \frac{1}{4}}\left(2-3\tau +{\tau }^3\right)\hfill & \hfill \hspace{1em}& {\varphi }_2\left(\tau \right)={\displaystyle \frac{1}{4}}\left(2+3\tau -{\tau }^3\right)\hfill \\ {\varphi }_3\left(\tau \right)={\displaystyle \frac{1}{4}}\left(1-\tau -{\tau }^2+{\tau }^3\right)\hfill & \hfill \hspace{1em}& {\varphi }_4\left(\tau \right)={\displaystyle \frac{1}{4}}\left(-1-\tau +{\tau }^2+{\tau }^3\right).\hfill \end{array}$$

By their definition, the projection functionals are

$$\begin{array}{ccc}{\rho }_1={\alpha }_i-g_i\left({\tau }_0\right)\hfill & \hfill \hspace{1em}& {\rho }_2={\alpha }_i-g_i\left({\tau }_f\right)\hfill \\ {\rho }_3={\displaystyle \frac{{\beta }_i}{b^2}}-g_i^{\prime }\left({\tau }_0\right)\hfill & \hfill \hspace{1em}& {\rho }_4={\displaystyle \frac{{\beta }_i}{b^2}}-g_i^{\prime }\left({\tau }_f\right)\hfill \end{array}$$

Collecting these terms in the form of Equation (15), the constrained expression and its derivatives can be expressed as

$$\begin{array}{cc}\hfill r_i(\tau ,g_i\left(\tau \right))=g_i\left(\tau \right)& +{\varphi }_1\left(\tau \right)\left({\alpha }_i-g\left({\tau }_0\right)\right)+{\varphi }_2\left(\tau \right)\left({\alpha }_i-g\left({\tau }_f\right)\right)\hfill \\ \hfill & +{\varphi }_3\left(\tau \right)\left({\displaystyle \frac{{\beta }_i}{b^2}}-g^{\prime }\left({\tau }_0\right)\right)+{\varphi }_4\left(\tau \right)\left({\displaystyle \frac{{\beta }_i}{b^2}}-g^{\prime }\left({\tau }_f\right)\right)\hfill \end{array}$$

(16)

In order to use this expression in our differential equation, we need to utilize the mapping parameter b from Equation (14), where

$$\begin{array}{c}\hfill \underset{{\dot{r}}_i}{\underbrace{{\displaystyle \frac{\mathrm{d}{r}_i}{\mathrm{d}t}}}}=b^2\underset{r_i^{\prime }}{\underbrace{{\displaystyle \frac{\mathrm{d}{r}_i}{\mathrm{d}\tau }}}}.\end{array}$$

3.2. Definition of the Free Function $g_j\left(\tau \right)$

For this application, let the free function be expressed as a linear combination of basis functions multiplied with unknown coefficients according to

$$g_i\left(\tau \right)={\displaystyle \sum _{j=1}^m}{a}_{ij} T_j\left(\tau \right)={\mathbf{\xi }}_i^T \mathit{h}\left(\tau \right)$$

(17)

where $\mathit{h}\in {\mathbb{R}}^m$ is a vector of the basis functions $T_j\left(\tau \right)$ and ${\mathbf{\xi }}_i\in {\mathbb{R}}^m$ is the unknown coefficient vector associated with the ith vector component. Like the original work on the solution of ordinary differential equations with TFC [6], we utilize Chebyshev orthogonal polynomials in the paper. The orthogonal polynomial representations are generally preferred due to their approximating and convergence properties. For example, Chebyshev polynomials generate a function that minimizes the maximum error in its application, and, therefore, they are well suited to approximating other functions [13,14].

However, due to the structure of the constrained expression, the free functions can be defined by a variety of function approximation methods, including other orthogonal polynomial sets and techniques in the field of machine learning. In fact, the application to neural networks (NN) was studied in Ref. [11], and promising results have been obtained by using the Extreme Learning Machine (ELM) [15], which is a single-layer feed-forward NN. The interested reader is directed to Ref. [16] for an in-depth look at the ELM implementation with TFC and Ref. [17] for applications of this technique to problems in spaceflight.

3.3. Discretization of the Domain

In order to solve problems numerically, the problem domain (and therefore the basis function domain) must be discretized. Since this article uses Chebyshev orthogonal polynomials, an optimal discretization scheme is the Chebyshev–Gauss–Lobatto nodes [18,19] shown in Equation (18). For $N+1$ points, the discrete points are calculated as

$${\tau }_k=-cos\left({\displaystyle \frac{k\pi }{N}}\right)\hspace{1em}\mathrm{for}\hspace{1em}k=0,1,2,\cdots ,N.$$

(18)

Compared with the uniform distribution, the collocation point distribution results in a much slower increase in the condition number of the matrix to be inverted in the least-squares as the number of basis functions, $m_2$, increases. The collocation points can be realized in the problem domain through the relationship provided in Equation (13).

3.4. Solution of the Resulting Algebraic Equation

By defining $g_i\left(t\right)$ according to Equation (17), our definition of Equations (2) and (3) is converted into four algebraic equations:

$${\tilde{F}}_i(\tau ,\mathsf{\Xi })=0 \mathrm{for} i=x,y,z,c$$

(19)

defined from Equation (5) through Equation (8). Note that the tilde symbol (∼) is used to signify that the differential Equation (5) through Equation (8) have been transformed into a new set of equations with embedded constraints. Additionally, $\mathsf{\Xi }$ represents the unknown parameters of this transformed set of equations, and these are defined in the following equations.

Next, we evaluate the three differential equations and one algebraic equation given by Equation (19) at the discretization points to construct a loss vector of the residuals of these equations.

$${\mathbb{L}}_i\left(\mathsf{\Xi }\right)={\left\{\begin{array}{ccccc}{\tilde{F}}_i({\tau }_0,\mathsf{\Xi }),& \cdots ,& {\tilde{F}}_i({\tau }_k,\mathsf{\Xi }),& \cdots ,& {\tilde{F}}_i({\tau }_f,\mathsf{\Xi })\end{array}\right\}}^{\mathrm{T}}={\mathbf{0}}^{\mathrm{T}}$$

with the total loss vector of

$$\mathbb{L}\left(\mathsf{\Xi }\right)={\left\{\begin{array}{cccc}{\mathbb{L}}_x^{\mathrm{T}}\left(\mathsf{\Xi }\right),& {\mathbb{L}}_y^{\mathrm{T}}\left(\mathsf{\Xi }\right),& {\mathbb{L}}_z^{\mathrm{T}}\left(\mathsf{\Xi }\right),& {\mathbb{L}}_c^{\mathrm{T}}\left(\mathsf{\Xi }\right)\end{array}\right\}}^{\mathrm{T}}={\mathbf{0}}_{4N\times 1}^{\mathrm{T}}$$

(20)

where the unknown vector, of size composed of $3m+7$ terms, is defined as

$$\mathsf{\Xi }={\left\{\begin{array}{cccccc}{\mathbf{\xi }}_x^{\mathrm{T}},& {\mathbf{\xi }}_y^{\mathrm{T}},& {\mathbf{\xi }}_z^{\mathrm{T}},& {\mathbf{\alpha }}^{\mathrm{T}},& {\mathbf{\beta }}^{\mathrm{T}},& b\end{array}\right\}}^{\mathrm{T}}={\mathbf{0}}^{\mathrm{T}}$$

Since the governing dynamics of the circular-restricted three-body problem are nonlinear, the unknown parameters occur nonlinearly in the loss vector, Equation (20). Therefore, a nonlinear least-squares technique was utilized using Equation (21) to update the parameters:

$${\mathsf{\Xi }}_{k+1}={\mathsf{\Xi }}_k+\Delta {\mathsf{\Xi }}_k$$

(21)

where Equation (22) is the classical least-squares equation using the Moore–Penrose inverse of the Jacobian of Equation (20):

$$\Delta {\mathsf{\Xi }}_k=-{\left({\mathbb{J}}^{\mathrm{T}}\left({\mathsf{\Xi }}_k\right) \mathbb{J}\left({\mathsf{\Xi }}_k\right)\right)}^{-1}{\mathbb{J}}^{\mathrm{T}}\left({\mathsf{\Xi }}_k\right) \mathbb{L}\left({\mathsf{\Xi }}_k\right),$$

(22)

This nonlinear least-squares method, known as the Gauss–Newton algorithm, is a modification of Newton’s method for finding a minimum of a function—in our case, minimizing the sum of the squared residual values in Equation (20).

The iteration continues until one of the following conditions is met:

$$max|\mathbb{L}\left({\mathsf{\Xi }}_k\right)|<\epsilon \hspace{1em}max|\Delta {\mathsf{\Xi }}_k|<\epsilon ,$$

based on the user defined tolerance, $\epsilon$, or the iteration exceeding some user-defined maximum.

4. Numerical Test

In the numerical implementation, all Jacobians were calculated using the JAX package [20,21] utilizing automatic differentiation [22]. In the nonlinear least-squares steps, NumPy’s pinv() was used to calculate the Moore–Penrose pseudo-inverse through singular-value decomposition (SVD). Additionally, a just-in-time compiler was used to optimize the code. All timing was handled by the process_time() from the Python package time. For all problems, we considered the Earth–Moon system with the parameters given in

Table 1. Earth–Moon system parameters.

Variable	Value
Earth mass $m_1$ (kg)	$5.9724\times {10}^{24}$
Moon mass $m_2$ (kg)	$7.3460\times {10}^{22}$

. Additionally, for the TFC implementation, the parameters used are summarized in

Table 2. TFC algorithm parameters.

Variable	Value
N [number of points]	140 [Lyapunov]  or  200 [Halo]
m [basis terms]	130 [Lyapunov]  or  190 [Halo]
$\epsilon$ [tolerance]	$2.22\times {10}^{-16}$
Maximum iterations	20
Least-squares method	numpy.pinv()

.

4.1. Initialization

For all numerical tests, the unknown vector must be initialized. First, the terms ${\mathbf{\xi }}_x$, ${\mathbf{\xi }}_y$, and ${\mathbf{\xi }}_z$ were all initialized by a null vector, which ultimately represents the simplest interpolating expression for the state variables. This initialization represents the worst-case scenario when there is no estimation of the trajectory, but just fitting the constraints. Next, the other unknown values of $\mathbf{\alpha }$, $\mathbf{\beta }$, and b (which are associated with the position, velocity, and the period of the orbit) were initialized using Richardson’s third-order analytical method for Halo-type periodic motion [23]. This initialization was used to find the first orbit of the specified Jacobi constants. For the following orbits, the desired Jacobi constant was incrementally increased, and the converged values from the prior Jacobi constant level were used to initialize each following step.

The same process was utilized for the differential corrector method [5], which was implemented to compare with TFC. In the differential corrector inner-loop, the desired Jacobi constant was obtained by an iterative least-squares approach to update the initial guess.

4.2. Lyapunov Orbits

First, the TFC method was used to explore the computation of Lyapunov orbits, which lie in the x-y plane, the orbital plane of the two primaries. For our test, the Lyapunov orbits were computed over a range of Jacobi constants, starting close to the equilibrium point’s specific energy levels and terminating at a Jacobi constant of 2.92. The associated trajectories for the orbits around L1 and L2 are provided in Figure 2a,b.

For the solution of Lyapunov orbits, it was important to transform the time for orbits of lower Jacobi constants, i.e., orbits closer to the secondary body, $m_2$. This was done by the time transformation from Ref. [5], where $t\to \zeta$, such that

$$\frac{\mathrm{d}t}{\mathrm{d}\zeta }=R_2.$$

(23)

This produces the following chain of transformations $t\to \zeta \to \tau$, or rather $[0,T]\to [0,{\zeta }_f]\to [-1,+1]$, which produces the final transformation:

$${\dot{r}}_i\left(t\right)={\displaystyle \frac{\mathrm{d}{r}_i}{\mathrm{d}t}}=\frac{\mathrm{d}{r}_i}{\mathrm{d}\tau }·\frac{\mathrm{d}\tau }{\mathrm{d}\zeta }·\frac{\mathrm{d}\zeta }{\mathrm{d}t}=\frac{b^2}{R_2}{r}_i^{\prime }\left(\tau \right)$$

To stress the importance of this scaling, Figure 3 shows both the scaled and non-scaled TFC solution for Lyapunov orbits around L1 (Figure 3a) and L2 (Figure 3b).

It can be seen that this scaling has a drastic effect on the accuracy of the converged solution. For orbits around L1, without scaling, significant accuracy loss is observed starting at a Jacobi constant of around 3.15. For orbits around L2, this accuracy reduction starts at a Jacobi constant of around 3.06, but is just as drastic.

Additionally, a comparison with the differential corrector method is provided in terms of speed and accuracy. Figure 4 compares the residuals for both methods, where it can be seen that the TFC approach is around two orders of magnitude more accurate than the differential corrector at higher Jacobi constants. Furthermore, the computation of the TFC solution is slightly faster, a little over 0.25 s in the extreme case, as displayed in Figure 5.

To further understand the computation time associated with the TFC approach, consider the computational breakdown given in Figure 6. Additionally, the statistics associated with this breakdown, in terms of mean and standard deviation, are given in

Table 3. Statistics associated with the computational breakdown given in Figure 6.

Variable	Average (ms)	Standard Deviation (ms)
Loss vector	5.69	3.46
Jacobian	91.1	55.8
Least-squares	329.8	203.8

. In the TFC method, the three major computations in each nonlinear least-squares iteration are evaluating the loss vector, the Jacobian, and the least-squares. Figure 6 shows that around 77% of the total computation time is associated with the least-squares portion. This should come as no surprise since the Jacobian is a 4 N × 2 m matrix—$560\times 260$ for this problem.

While the least-squares, i.e., matrix inversion, is the main bottleneck with the computational efficiency, a marginal speed gain could be achieved by removing the reliance on automatic differential through JAX [20,21] for the computation of the loss vector and Jacobian. Admittedly, these two could be “hard-coded” since the partial derivatives populating the Jacobian can be analytically derived. However, the speed gain associated with this is unknown and should be the focus of future implementations and numerical tests.

Similar to the L1 Lagrange point test, Figure 2a displays the computed trajectories around the L2 Lagrange point. Additionally, as in Figure 4 and Figure 5, the accuracy and computation time for these tests are provided in Figure 7 and Figure 8. In Figure 7, the TFC method is more accurate, albeit only slightly. At a Jacobi constant level of around 3.00 and above, the differential corrector method does not converge, as shown by the jump in accuracy. At this Jacobi constant level, the TFC method’s accuracy starts to decrease before failing to converge at the Jacobi constant value of 2.92. This drop-off in the accuracy of the TFC method can easily be seen by analyzing Figure 8, where the red box highlights where the TFC method’s max iteration limit is reached. In other words, at the Jacobi constant levels lower than 3.00, the representation of the orbit through the free function composed of Chebyshev polynomials starts to fail.

In both Lyapunov orbits around L1 and L2, the reduction in accuracy can be attributed to the definition of the free function, $g\left(x\right)$, not fully capturing the solution space of the trajectory. As the Jacobi constant decreases, the complexity of the trajectory’s shape increases and therefore an increasing number of basis functions is needed. However, contrary to this, increasing the basis terms to more than 130 does not provide a more accurate solution in the numerical implementation. Therefore, when solving lower Jacobi constant level Lyapunov orbits, a different definition of the free function must be used; future effort should focus on exploring other basis types.

4.3. Halo Orbits

Next, the proposed technique was utilized to compute Halo orbits around L1 and L2. These orbits differ from the Lyapunov orbit because they are not restricted to the x-y plane and are therefore three-dimensional. In fact, this family of orbits is a bifurcation of the Lyapunov orbits computed in the previous section and is characterized by “northern” and “southern” bifurcations. However, using the TFC method to compute these Halo orbits, the only difference is in the initialization of the $\mathbf{\alpha }$, $\mathbf{\beta }$, and b parameters. Additionally, for Halo orbits, the time transformation introduced in Equation (23) was not used because it produced unfavorable convergence properties, i.e., the method would simply converge to the Lyapunov family of orbits. First, we look at the computation of the “northern” family of Halo orbits around L1 and L2, as plotted in Figure 9.

In these plots, we can see that, around the L1 equilibrium point, the method converged to Lyapunov orbits for higher Jacobi constants. However, as the Jacobi constant decreases below 3.025, the method does not converge to a periodic orbit, as shown by the increase in residuals around the Jacobi constant value of 3.025. At the other equilibrium point, L2, the method converges for all Jacobi constant values; however, at values below 3.025, the solution jumps to a circular orbit around the Moon. Therefore, these orbits were not plotted. This “jumping” from the Halo orbit bifurcation to the circle orbit bifurcation is the major drawback with the algorithm in its current state; however, it is noted that the differential corrector approach also suffers from this. The reason that the TFC method exhibits this behavior is that the numerical solution, i.e. the nonlinear least-squares optimization step, has no information on the desired bifurcation. Therefore, the solution becomes closely linked to (1) the initialization of the unknown parameters and (2) the definition of the free function $g\left(x\right)$. In other words, the initialized parameters will converge to the closest local minimum.

Like the Lyapunov orbit tests, the loss vector’s maximum residual was recorded and is plotted in Figure 10. For almost all converged solutions, the residuals were on the order of $\mathcal{O}\left({10}^{-14}\right)$; however, around a Jacobi constant level of 3.025, the accuracy decreases for orbits around L1. The convergence for Halo orbits took longer, with some cases taking 8 s, while, on average, the solution time was around 2 s, as shown in Figure 11. Similar to the “northern” Halo orbits plotted in Figure 9, the Halo orbits of the “southern” bifurcation were computed with similar findings and, therefore, omitted from this paper for brevity.

5. Conclusions

In this paper, we implemented a functional interpolation method, based on the Theory of Functional Connections (TFC), to computed periodic orbits, i.e., Lyapunov and Halo orbits, in the circular-restricted three-body problem. The TFC approach allowed for the analytical embedding of the periodic constraint, such that the orbit’s initial and final position and velocity repeated after some period T. This process reduces the search space of the numerical algorithm to only those trajectories satisfying the periodic constraint.

In general, the method produced comparable results to the differential corrector method; however, the ease of implementation of the TFC method provides one benefit. For example, when using the TFC method, the user only needs to code the loss vector, Equation (20), and the parameter update is handled through automatic differentiation. For readers interested in this implementation, the authors direct them to the TFC GitHub (https://github.com/leakec/tfc, accessed on 15 February 2021) [24], where both codes for the Lyapunov (https://github.com/leakec/tfc/tree/main/examples/Hunter_Johnston_Dissertation/Chapter_4/Example_4_8, accessed on 15 February 2021) and Halo (https://github.com/leakec/tfc/tree/main/examples/Hunter_Johnston_Dissertation/Chapter_4/Example_4_9, accessed on 15 February 2021) orbit cases are provided. In addition to the code used to produce the results of this paper, this toolbox provides examples to the solution of a wide variety of ODEs and PDEs.

In all, this work was the first time that the TFC method was utilized to enforce the constraints of periodic orbits. Future work in this area should focus on two major aspects in order to increase speed, accuracy, and robustness of the technique, which are summarized below:

• Definition of the free function. In this study, Chebyshev orthogonal polynomials were used to define the free function; however, a large number of basis terms were needed to accurately describe the trajectories for both Lyapunov and Halo orbits. One idea to remedy this is to use a hybrid basis composed of terms to capture the periodic and non-periodic portions separately.
• Faster least-squares implementation. As seen in Figure 6, the major bottleneck with computation time lies in the least-squares portion, which is an order of magnitude slower than the evaluation of the loss vector and Jacobian combined. The authors acknowledge that a more efficient algorithm than NumPy’s pinv() could be implemented.