Fast 2-impulse non-Keplerian orbit-transfer using the Theory of Functional Connections
Allan K. de Almeida Junior, Hunter Johnston, Carl Leake, Daniele Mortari
FFast 2-impulse non-Keplerian orbit-transfer using the Theory of FunctionalConnections (cid:63)
Allan K. de Almeida Junior · Hunter Johnston · CarlLeake · Daniele MortariAbstract
This study applies a new approach, the Theory of Functional Connections (TFC), to solve the two-point boundary-value problem (TPBVP) in non-Keplerian orbit transfer. The perturbations considered aredrag, solar radiation pressure, higher-order gravitational potential harmonic terms, and multiple bodies. Theproposed approach is applied to Earth-to-Moon transfers, and obtains exact boundary condition satisfactionand with very fast convergence. Thanks to this highly efficient approach, perturbed pork-chop plots of Earth-to-Moon transfers are generated, and individual analyses on the transfers’ parameters are easily done at lowcomputational costs. The minimum fuel analysis is provided in terms of the time of flight, thrust applicationpoints, and relative geometry of the Moon and Sun. The transfer costs obtained are in agreement with theliterature’s best solutions, and in some cases are even slightly better.
Keywords
Two point boundary value problem · Transfers in space · Bi-circular four body problem · Theoryof Functional Connections
The main problem in evaluating perturbed orbit transfer costs is solving a two-point boundary-value problem(TPBVP) on a nonlinear dynamics. This problem is traditionally solved using shooting methods. This paperfocuses its attention on Earth-to-Moon co-planar orbit transfer as an example of TPBVP orbit transfer byalso including Sun gravitational perturbation.The three-body problem (3BP) [38,7,13] has traditionally been used as an initial guess to solve multi-body astrodynamics problems. Earth/Moon transfers have recently attracted the interest of the scientificcommunity following the work done by [42]. This transfer can be done through the use of low thrust maneuvers[35,12,11] or two-impulse thrust, a Hohmann maneuver, solved via the patched restricted three-body problem[6,34], full ephemeris n -body problem [2,21], restricted three body-problem [36,48,44,23,6], restricted four- body problem [1], and bi-circular restricted four-body problem [50,49,24,5,26,25,45,33,31,40,32].The Earth/Moon transfer TPBVP has also been solved using a gradient method based on the Lambertproblem as an initial guess [39], a direct transcription and multiple shooting method [45], or a more recentlydeveloped method based on the gradient method and Jacobi integral variational equation [8].In general, most of the techniques used to solve this problem in the mentioned references, such as theshooting method, are computationally intensive, because they need to integrate the differential equations (cid:63) A previous version was presented as: De Almeida Junior, A.K., Johnston, H., Leake, C., and Mortari, D. “Evaluationof transfer costs in the Earth/Moon system using the Theory of Functional Connections,” AAS 20-596, AstrodynamicsSpecialist Conference, August 9-12, Lake Tahoe, CA.Allan K. de Almeida JuniorINPE - National Institute for Space Research, S˜ao Jos´e dos Campos, SP, BrazilAerospace Engineering, Texas A&M University, College Station, TX.E-mail: [email protected] JohnstonAerospace Engineering, Texas A&M University, College Station, TX. E-mail: [email protected] LeakeAerospace Engineering, Texas A&M University, College Station, TX. E-mail: [email protected] MortariAerospace Engineering, Texas A&M University, College Station, TX. E-mail: [email protected] a r X i v : . [ a s t r o - ph . E P ] F e b Allan K. de Almeida Junior et al. (dynamics) multiple times, because the shooting method iterates the initial conditions to match the bound-ary conditions. In this study, the authors adopt the recently developed
Theory of Functional Connections (TFC) [28,27,29,30,18] to solve the TPBVP. This method provides functionals that analytically satisfy theboundary conditions. This way, the initial constraint optimization problem becomes unconstrained, and asingle unconstrained nonlinear least-squares approach can be applied to solve the optimization problem. Incontrast, previous methods need to solve the problem multiple times.The efficiency also involves the accuracy in terms of position and velocity errors. In particular, the codeuses automatic differentiation and a just-in-time (JIT) compiler [3,9] to implement the nonlinear least-squares. As a result, the run time to find a solution for a given set of parameters is on the order of seconds.This fast method for solving the TPBVP with nonlinear differential equations has facilitated the generationof multi-body, perturbed, orbit transfer pork-chop plots and individual analyses on the transfer parameters.The two-impulse transfer is done using the circular co-planar restricted 3BP and bi-circular, bi-planarfour-body problem by including the Sun. A set of four parameters of this transfer are investigated to identifythe optimal combination that minimizes fuel consumption. These parameters are 1) first impulse applicationpoint on the departure LEO, 2) second impulse application point on the arrival orbit around the Moon,3) time of flight, and 4) Sun relative direction with respect to the Earth/Moon system. In particular, thetransfer time considered is short, up to seven days.This article is organized as follows. The mathematical background, the definition of the problem, andthe models adopted are shown (next section). In the results section, the parameters’ influence over theequivalent fuel cost are analyzed for a transfer between a LEO and another circular orbit close to the Moon.Such results are compared with those available in the literature, and the major findings are summarized inthe conclusions section.
In this section, the problem is defined, and the mathematical models and tools are explained. First, arecent technique to solve the two-point boundary value problems for a system of differential equations issummarized. Then, the problem is defined, and the mathematical model used to evaluate the transfer costsis presented.2.1 Summary of the Theory of Functional ConnectionsThe numerical technique to solve differential equations based on the Theory of Functional Connections(TFC) [28,30] has been shown to produce solutions with machine-level error within milliseconds for bothlinear [27] and nonlinear [29] differential equations and for a wide array of constraint cases. Several workshave utilized the TFC framework including, hybrid systems [15], optimal control problems [10,16], quadraticand nonlinear programming [22], and other applications [14,20,19].
This approach’s foundation 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 con-strained expressions, the differential equation can be transformed into an algebraic expression that canultimately be solved via a variety of optimization techniques. In general, the TFC methodology provides functional interpolation by embedding a set of n linear constraints. The general expression from [18] can beadopted for vector equations using the following functional, x ( t, g ( t )) = g ( t ) + k (cid:88) j =1 φ j ( t ) ρ j ( t, g ( t )) , (1)where g ( t ) is a free function and ρ j ( t, g ( t )) are called projection functionals since they project the free functiononto the space of functions that satisfy the constraints. They are derived by setting the constraint functionequal to zero and replacing x ( t ) with g ( t ). In addition, φ j ( t ) are called switching functions; by definition,they are equal to 1 when evaluated at the constraint they are referencing and equal to 0 when evaluatedat all other constraints. The switching functions are composed of a set of linearly independent functionscalled support functions, s k ( t ), with unknown coefficients α ij , such that φ j ( t ) = s i ( t ) α ij . The following step-by-step procedure can be used to derive the switching functions:1. Choose k support functions, s k ( x );2. Write each switching function as a linear combination of the support functions with unknown coefficients; ast 2-impulse non-Keplerian orbit-transfer using the Theory of Functional Connections (cid:63)
3. Based on the switching function definition, write a system of equations to solve the unknown coefficients.Once the switching functions and projection functionals have been constructed, they can be substituted intoEq. (1) to form the constrained expression. For the interested reader, this process is described in far moredetail with examples in [18].2.2 Application of the Theory of Functional Connections to solve two-point boundary value problemsFor convenience, a step-by-step derivation of the constrained expression used in the paper is presentedbelow. Additionally, a summary of the numerical technique is also provided. Since we deal with a bi-planar,bi-circular four-body problem, a two-dimensional vector is defined as X = { x, y } T , where X = x and X = y .This problem’s dynamics are governed by a nonlinear second-order vectorial differential equation subject toboundary value and derivative constraints. In general, a differential equation of this kind can be expressedas F i ( t, X j , ˙ X j , ¨ X j ) = 0 subject to: (cid:40) X j ( t ) = X j X j ( t f ) = X f j where i, j = 1 , . (2)Since the constraints for each component of the vector X i are of the same form, Eq. (1) can easily be adaptedto the vectorial form X i ( t, g i ( t )) = g i ( t ) + k (cid:88) j =1 φ j ( t ) ρ j ( t, g i ( t )) . (3)Following the process detailed in the previous section, one must first choose the support functions s k ( t ). Forthese specific constraints, the support functions are defined in the same manner as [16], where s k ( t ) = t ( k − .Now, the definition of the switching functions is used to produce a set of equations. For example, the firstswitching function has the two equations, φ ( t ) = 1 and φ ( t f ) = 0 . These equations are expanded in terms of the support functions, φ ( t ) = (1) · α + ( t ) · α = 1 φ ( t f ) = (1) · α + ( t f ) · α = 0which can be compactly written as (cid:20) t t f (cid:21) (cid:26) α α (cid:27) = (cid:26) (cid:27) . The same is done for the other switching function to produce a set of equations that can be solved by matrixinversion, (cid:20) t t f (cid:21) (cid:20) α α α α (cid:21) = (cid:20) (cid:21)(cid:20) α α α α (cid:21) = (cid:20) t t f (cid:21) − = t f t f − t − t t f − t − t f − t t f − t , leading to the switching functions, φ ( t ) = t f − tt f − t and φ ( t ) = t − t t f − t . From their definition, the projection functionals are ρ ( t, g i ( t )) = X i − g i ( t ) and ρ ( t, g i ( t )) = X i f − g i ( t f ) . Collecting these terms in the form of Eq. (3), the constrained expression can be expressed as X i ( t, g i ( t )) = g i ( t ) + t f − tt f − t (cid:16) X i − g i ( t ) (cid:17) + t − t t f − t (cid:16) X i f − g i ( t f ) (cid:17) . Allan K. de Almeida Junior et al.
Now, substituting the constrained expression into the differential equation, Eq. (2), the differential equationis transformed into one subject to no constraints,˜ F i ( t, g j ( t ) , ˙ g j ( t ) , ¨ g j ( t )) = 0 i, j = 1 , . (4)In order to solve this new equation, three major steps must be taken: 1) define the free-function g i ( t ), 2)discretize the domain, and 3) solve the resulting algebraic equation. The following text summarizes thisprocess. g j ( t ) : For this application, let the free-function be expressed as a linear combination of basis functions multipliedwith unknown coefficients according to g i ( t ) = h T ξ i , (5)where h ∈ R m is a vector of the basis functions and ξ i ∈ R m is the unknown coefficient vector associated withthe i th vector component. In general, the basis functions that make up h will not coincide with the problemdomain, e.g., Chebyshev and Legendre polynomials are defined in [ − , z ∈ [ z , z f ] and the problem be defined for t ∈ [ t , t f ]. In order to use the basis functions, amap between the basis function domain, z , and problem domain, t , must be created. The simplest map islinear, z = z + z f − z t f − t ( t − t ) ←→ t = t + t f − t z f − z ( z − z ) . (6)Using this map, the subsequent derivatives of the free function can be computed asd n g i d t n = (cid:18) d z d t (cid:19) n d n h T d z n ξ i , where, by defining c := d z d t = z f − z t f − t , the derivative computations can be simply written as d n g i d t n = c n d n h T d z n ξ i , where c is referred to as the mapping coefficient. In order to solve problems numerically, the problem domain, and therefore the basis function domain, mustbe discretized. Since this article uses Chebyshev orthogonal polynomials, an optimal discretization schemeis the Chebyshev-Gauss-Lobatto nodes [17,47]. For, N + 1 points, the discrete points are calculated as, z k = − cos (cid:18) kπN (cid:19) for k = 0 , , , · · · , N. (7)When compared with the uniform distribution, the above collocation point distribution results in a much easier condition of the matrix to be inverted in the least-squares as the number of basis functions, m ,increases. The collocation points can be realized in the problem domain through the relationship providedin Eq. (6). ast 2-impulse non-Keplerian orbit-transfer using the Theory of Functional Connections (cid:63) By defining g i ( t ) according to Eq. (5), our definition of Eq. (4) is converted into˜ F i ( t, ξ j ) = 0 i, j = 1 , . (8)Now, if the original differential equation (2) is linear, then Eq. (8) is linear in the unknown ξ j coefficientsand can be solved with any least-squares technique. However, if the original differential equation is nonlinear,as is the governing differential equation of this paper, a nonlinear least-squares approach must be used. First,using Eq. (8) one must construct the loss function of the system, which takes the form, L ( Ξ ) = ˜ F ( t , Ξ )...˜ F ( t f , Ξ )˜ F ( t , Ξ )...˜ F ( t f , Ξ ) = where Ξ = { ξ T ξ T } T and the function ˜ F i ( t, ξ j ) is evaluated at the Chebyshev-Gauss-Lobatto nodes throughEqs. (6) and (7). In order to solve for the unknown coefficients, the parameter update equation is providedby, Ξ k +1 = Ξ k + ∆Ξ k (9)where the ∆Ξ k can be defined using classic least-squares, ∆Ξ k = − (cid:16) J T ( Ξ k ) J ( Ξ k ) (cid:17) − J T ( Ξ k ) L ( Ξ k )or in this paper through a QR decomposition method. It follows that the Jacobian matrix is, J ( Ξ ) = ∂ ˜ F ( t ,Ξ ) ∂ ξ ∂ ˜ F ( t ,Ξ ) ∂ ξ ... ... ∂ ˜ F ( t f ,Ξ ) ∂ ξ ∂ ˜ F ( t f ,Ξ ) ∂ ξ ∂ ˜ F ( t ,Ξ ) ∂ ξ ∂ ˜ F ( t ,Ξ ) ∂ ξ ... ... ∂ ˜ F ( t f ,Ξ ) ∂ ξ ∂ ˜ F ( t f ,Ξ ) ∂ ξ . q is defined such that it locates the origin of a rotating frame of reference with respect to an inertialone. The equation of motion of a particle in this rotating frame of reference is [43] (chapter 7)d r d t + 2 ω × d r d t + ω × ( ω × r ) + d ∗ ω d t × r + d ∗ q d t = a , (10)where r locates the spacecraft in the rotating frame, ω is the angular velocity of the rotating frame, and a isthe specific force acting on the particle. In this equation, the derivatives of the fourth and fifth members ofthe left side must be taken with respect to the time in the inertial frame, while the others are taken in therotating frame - not taking into consideration the motion of their bases with respect to the inertial frame.In the bi-planar, bi-circular four-body problem, the Earth and the Moon rotate in a circular orbit aroundtheir barycenter, which in turn rotates in a circular orbit around the Sun. Moreover, the two orbital planescoincide with each other. The rotating frame of reference can be seen in Fig. 1. In our problem, the centerof the inertial frame of reference coincides with the Sun’s position. The vector q locates the barycenter ofthe Earth and the Moon, which are located along the x axis of the rotating frame of reference. The position of the origin of the rotating frame (the barycenter of the Earth/Moon system) written in the inertial frameis q = { R s cos θ , R s sin θ } T , where R s is the distance between the Sun and the barycenter of the pairEarth/Moon and θ = ( ω r t + γ − π ) is the angle that locates the center of the rotating frame with respectto the positive side of the x -axis of the inertial frame, where ( γ − π ) is an initial phase and ω r is the Allan K. de Almeida Junior et al. r s MoonEarth L xy r r Spacecraft ω r MoonEarth L xy Spacecraft ω Sun q γ Fig. 1
The Earth/Moon rotating frame of reference. constant angular velocity of the barycenter of the Earth/Moon system in the inertial frame of reference. Theacceleration of this barycenter with respect to the inertial frame written in the rotating frame of referenceis d ∗ q d t = µ s R s (cid:26) cos ( ω s t + γ )sin ( ω s t + γ ) (cid:27) , (11)where µ s is the gravitational parameter of the Sun and ω s = ω r − ω . The distance of the Earth and theMoon from their center of mass is given by d = Rµ / ( µ + µ ) and d = Rµ / ( µ + µ ), respectively, where µ and µ are the gravitational parameters of the Earth and the Moon, respectively, and R = d + d is thetotal distance between them. The position vectors of the spacecraft relative to the Earth ( r ) and the Moon( r ) can be written as r = r + d ˆ i and r = r − d ˆ i , where [ˆ i , ˆ j ] are unit vectors along the [ x, y ] axes,respectively. Their magnitudes are r = (cid:107) r (cid:107) and r = (cid:107) r (cid:107) . Similarly, the position vector of the spacecraft relative to the Sun is given by r s = r − p s , where p s = R s { cos ( ω s t + γ ) , sin ( ω s t + γ ) } T locates the Sun inthe rotating frame of reference.Hence, the equation of motion, Eq. (10), of a spacecraft under the gravitational influence of the Earth,the Moon, and the Sun written in the rotating frame of reference isd r d t + 2 ω × d r d t + ω × ( ω × r ) = − µ r r − µ r r (12) − µ s r s r s − µ s R s { cos ( ω s t + γ ) , sin ( ω s t + γ ) } T . In the case where µ s = 0, Eq. (12) is reduced to the planar circular restricted three-body problem. Thissimpler case will be investigated first in section 3. Afterwards, the bi-circular, bi-planar restricted four-bodyproblem will be used to model the spacecraft’s motion.2.4 The equivalent fuel cost ∆V The initial and final positions of the transfer are denoted with A and B , respectively. A first impulse isapplied to the spacecraft at the point A such that it travels to the point B during a time of flight T underthe gravitational influence of the main bodies. The initial and final velocities ( V A and V B , respectively) of ast 2-impulse non-Keplerian orbit-transfer using the Theory of Functional Connections (cid:63) the boundary value problem in the rotating frame are evaluated through the use of the TFC. The equivalentfuel cost ∆V is the sum of the magnitude of the initial and the final impulses, and hence it is given by ∆V = (cid:107) δ V A (cid:107) + (cid:107) δ V B (cid:107) , where (cid:107) δ V A (cid:107) , the magnitude of the difference between the initial velocity and the one that the spacecrafthad just before the impulse, is given by (cid:107) δ V A (cid:107) = (cid:107) V A − V Ai (cid:107) , and (cid:107) δ V B (cid:107) , the magnitude of the difference between the final velocity and the velocity required for thespacecraft to achieve its final orbit, is similarly given by (cid:107) δ V B (cid:107) = (cid:107) V Bf − V B (cid:107) . r around the Earth with an angular velocity ω = (cid:112) µ /r in thisinertial frame of reference. A second inertial frame of reference is also defined such that its center coincideswith the center of the Moon. The spacecraft will be inserted in a second orbit around the Moon of radius ρ and angular velocity Ω = ± (cid:112) µ /ρ in this inertial frame. The first impulse is applied when the spacecraft’sinitial position vector relative to the Earth in the circular orbit makes an angle α with the x -axis. The secondimpulse is applied when the spacecraft crosses the final circular orbit around the Moon at the point B , atwhich point the orbit makes an angle β with the x -axis in the rotating frame of reference. Both points, A and B , are shown in Fig. 2. initial orbit final orbitA BEarth Moon α β xy r ρ Fig. 2
The initial (A) and final (B) points of the transfer.
The position vector at the initial time of motion in the rotating frame of reference is r A = (cid:26) − d + r cos αr sin α (cid:27) . The initial velocity before the impulse in the rotating frame of reference is V Ai = (cid:26) − ( ω − ω ) r sin α ( ω − ω ) r cos α (cid:27) . The position vector at the final time of motion in the rotating frame is r B = (cid:26) d + ρ cos βρ sin β (cid:27) Allan K. de Almeida Junior et al. and the final velocity of the circular orbit after the second impulse (at B ) in the rotating frame is V Bf = (cid:26) − ( Ω − ω ) ρ sin β ( Ω − ω ) ρ cos β (cid:27) . In the case of a counter-clockwise final orbit around the Moon, the angular velocity is Ω = (cid:112) µ /ρ . In thecase the orbit is clockwise, the angular velocity is Ω = − (cid:112) µ /ρ .2.6 Error estimationThe initial and final positions are known, and the initial and final velocities are obtained through the TFCmethod. The initial position and initial velocity form a set of initial conditions. This set of initial conditionsis numerically integrated for the same system via a fourth-order Runge-Kutta integrator. A position error P E is defined as the magnitude of the difference between the final position (the point B ) and the position givenby the Runge-Kutta integrator ran for the corresponding time of travel T . Analogously, a velocity error V E is defined as the magnitude of the difference between the velocity at the final time of motion t f obtained bythe solution via TFC and the Runge-Kutta method using the initial conditions given by the TFC solution.Note that P E and V E are not true errors since they also accumulate numerical errors generated by the RKintegrator. On the other hand, there is no true solution available for this system, and such a combination oferrors may give us a good estimation of the upper limits of the true errors for this model.2.7 Values of the parametersThe values of the parameters for the Earth/Moon system used in this research are given in Table 1. Tofacilitate comparisons, these values are the same as those used by [41,48,49,45]. Table 1
Values of the parameters for the Earth/Moon system [41]. R . × m R s . × m µ . × m / s µ . × m / s µ s . × m / s ω . × − ω s − . × − . × mRadius of the Moon 1 . × m The results for the transfer are presented in this section. The spacecraft departs from an initial circular orbitaround the Earth with an altitude of 167 km and will be inserted into a final orbit around the Moon withan altitude of 100 km, according to Fig. 2.A spiral is used as an initial guess to find the possible solutions of the system. The spiral is defined asa linear variation on time of the radius and the angle in polar coordinates, from the initial position to thefinal one.
Each iteration of Eq. (9) takes less than a tenth of a second to be evaluated (usually much less) for aresidual on the order of 10 − . For most cases, a solution is found between 10 to 20 iterations. The caseswhere no convergence is obtained are exceptional. The number of iterations depends on the 3BP/4BP caseconsidered, on the boundary values, on the time of flight, and, obviously, on the initial guess adopted. ast 2-impulse non-Keplerian orbit-transfer using the Theory of Functional Connections (cid:63) µ s = 0 inEq. 12. This model is considered in this subsection. In the CR3BP, this system has three main parametersrelated to the trajectory of the transfer: the position to apply the first impulse in the initial orbit, describedby the parameter α , the position at which the spacecraft applies the second impulse near the Moon, β , andfinally, the transfer time T . In order to investigate the relations of all of these parameters with the best ∆V , the following procedure is used. First, the parameter T is fixed; thus, the dependency of α and β onthe minimum ∆V is to be found for a specific time of flight T . The minimum ∆V is numerically found bysweeping the parameters α and β in their neighborhood one at a time and maintaining the lowest valueamong the ones numerically found. After that, the neighborhood in which each parameter is evaluated isdecreased, and the process is repeated several times for all the parameters. The whole process is repeateduntil the change in the ∆V is less than 1 × − m/s, and the change in the angles are less than 1 × − rad.Although the position error of the solution obtained through TFC is always smaller than one meter (aboutfive orders of magnitude smaller than one meter in most cases), the obtained solution is not necessarily theone with the lowest ∆V . Thus, the main source of error in the procedure to find the parameters’ relationsis the convergence. We deal with a complex non-linear problem with multiple solutions, and the solutionmay converge to ones that are not of interest. Thus, although the process stops when the above changesin ∆V and angles are reached, it does not necessarily mean that this range is the error for the solution ofinterest, i.e., the minimum ∆V for that time of transfer T . On the other hand, the accuracy is good enoughto investigate the parameters’ variations as functions of the time of flight for the best ∆V , as can be seen inthe following results.Using the process described above, the trajectories with the lowest costs are found, and the parametersrelated to these trajectories are disclosed. The results can be seen in Fig. 3, where the dependencies of theparameters α (in red) and β (in blue) on the time of flight T are shown for a final counter-clockwise orbit(above) and a final clockwise orbit (below). The values of the transfer costs are shown in gray in km/s. Thered and blue curves are curve fits for the dots using tenth order Hermite polynomials. The parameter α decreases as T increases for short-time transfers up to two days and increases with T after this value in arange from 222 ◦ to 272 ◦ . This behavior for α is identical for both the counter-clockwise and the clockwisecases. In fact, there is no considerable difference in the values of this parameter for both cases. In thecounter-clockwise case, the parameter β increases as T increases for T up to almost three days and decreaseswith T after this value in a range from 270 ◦ to 208 ◦ . Considering the clockwise case, the parameter β alwaysdecreases when T increases in a different range from about 70 ◦ to -70 ◦ (or 290 ◦ ). The best ranges for α and β can be seen in Fig. 4. Using the procedure described above, all parameters are varied to search for theminimum ∆V , including the time of flight T . The values of all parameters associated with the minimum ∆V are shown in Table 2 for the counter-clockwise and clockwise cases. Variables associated Counter-clockwise Clockwise ∆V (m / s) 3946 .
93 3952 . (cid:107) δ V A (cid:107) (m / s) 3134 .
60 3137 . (cid:107) δ V B (cid:107) (m / s) 812 .
33 814 . T (days) 4 . . α (rad) 4 . . α . ◦ . ◦ β (rad) 4 . . β . ◦ . ◦ P E (m) 4 . × − . × − V E (m / s) 2 . × − . × − ˙ x | t = t (m / s) 9745 .
19 10007 . y | t = t (m / s) − . − . Table 2
Values associated to the minimum ∆V found for the transfer in the CR3BP. In order to analyze the influence of the Sun, the complete system given by Eq. (12) is considered with theparameters given in Table 1. To analyze the influence of the Sun over the trajectory related to the minimum ∆V for the previous CR3BP, the parameters α , β , and T shown in Table 2 are fixed, and the angle related αβ ° ° ° ° ° ° ° T ( days ) A ng l e s α and β Δ V ( km / s ) Δ V m i n ( k m / s ) αβ ° ° ° ° T ( days ) A ng l e s α and β Δ V ( km / s ) Δ V m i n ( k m / s ) Fig. 3
Relations between all of the parameters in CR3BP transfer: the angles α (in red) and β (in blue) associated withthe trajectory with the minimum ∆V found (in gray) as functions of the time of flight T . The dots are the values foundin this research, while the curves are fits to these dots. The upper figure is the counter-clockwise case and lower one is theclockwise. to the position of the Sun is varied. The cost ∆V as a function of γ is shown in Fig. 5. Note that the ∆V changes by about 2 m/s due to the Sun’s influence over the transfer compared with the CR3BP. There aretwo ∆V local minima for γ around 95 ◦ or 275 ◦ . The dependencies of all the parameters α , β , and γ on thetime of flight T are shown in Fig. 6 for the local minimum of γ around 95 ◦ (upper) and for the local minimumof γ around 275 ◦ (lower): both for the counter-clockwise case. The clockwise case is shown in Fig. 7 for thesame two best positions of the Sun. The dots correspond to the found minima, and the curves represent fitting to the points using tenth order Hermite polynomials. The values related to the lowest ∆V found inthis research for the bi-planar, bi-circular 4BP are shown in Table 3. The values of α and β that producethe lowest ∆V differ by up to 0.7 ◦ from the CR3BP: this can be seen by comparing the values reported inTables 2 and 3. The Sun’s influence changes the cost of the transfers by about 2 m/s, depending on the Sun’s ast 2-impulse non-Keplerian orbit-transfer using the Theory of Functional Connections (cid:63) Earth Satellite in theinitial orbitRangeof α To theMoon
MoonSatellite in thefinal orbit
Range of β counter - clockwise Range of β clockwise To theEarth
Fig. 4
Range of the angles of departure α and arrival β that minimizes the costs ∆V for the counter-clockwise ( β c − c ) andclockwise ( β clockwise ) cases for transfers from 1 to 7 days. position. There are two minima and two maxima separated by a 180 ◦ phase with respect to the positionof the Sun, γ : the position of the Sun for the lowest ∆V is about 96 ◦ or 276 ◦ (for T ≈ . ◦ or 300 ◦ for longer transfers times of about seven days. Thevariation of γ for the counter-clockwise case shown in Fig. 6 is also shown in an isolated form in Fig. 8 for abetter visualization of the range of γ . The range of the values of γ that minimizes the fuel cost is shown inFig. 9.Since the values of α , β , γ , and T are disclosed in Figs. 6 and 7 for the lowest ∆V , the effects of thevariations of pairs of these parameters can now be investigated. The first analysis is on the pair α and T . Inthis case, the values of β and γ are constrained using the curve given by the Hermite polynomials fits of thedots in Fig. 6 (upper). The cost ∆V as a function of α and T is shown in Fig. 10. The red curve in this figureis the same red curve of Fig. 6 (upper), and it shows the values of α as a function of T for the minimum ∆V . Note that the data shown in Fig. 10 is obtained for the counter-clockwise case, but the pattern of α isnot changed for the clockwise case. The second analysis is on the pair β and T . In this case, the values of α and γ as functions of T are given by the Hermite polynomial fits shown in Fig. 6 (upper). In this secondcase, the cost ∆V as a function of the pair β and T is shown in Fig. 11. The blue curve represents the valuesof β for the minimum ∆V , and it is the same curve shown in Fig. 6 (upper). Note that this analysis is alsodone for the counter-clockwise case. The third analysis shows the cost ∆V as a function of the pair β and T ,but this time for the clockwise case. In this case, the parameters α and γ are given by the curve fits of Fig.7 (upper). The result is shown in Fig. 12. The blue curve represents the same blue curve of Fig. 7 (upper), and it represents the values of β that should be used to minimize the cost ∆V . Note in Figs. 10 - 12 that asmall variation in the parameters α and β leads to a large variation in the cost ∆V . ° ° ° ° ° ° ° ° γ Δ V ( k m / s ) Fig. 5
The ∆V cost as a function of the position of the Sun γ for the parameters α , β and T fixed in the values given inTable 2.2 Allan K. de Almeida Junior et al. αβγ ° ° ° ° ° T ( days ) A ng l e s α , β , and γ Δ V ( km / s ) Δ V m i n ( k m / s ) αβγ ° ° ° ° ° T ( days ) A ng l e s α , β , and γ Δ V ( km / s ) Δ V m i n ( k m / s ) Fig. 6
The counter-clockwise case: dependencies among all of the parameters in 4BP transfer: the angles α (in red), β (inblue), and γ (in green) associated with the trajectory with the minimum ∆V found (in gray) as functions of the time offlight T . T - ∆V available in the literature. The data for the best ∆V and their respective time of flight T available in the literature for a time of flight up to seven days are shownin Table 4. In the Hohmann-like transfer, the spacecraft is under the influence of only the Earth or only theMoon gravitational attraction in the trajectory’s initial and final parts, respectively. This type of transfer islargely used to compare fuel cost efficiency ∆V and time of flight [36,2,4,37,45,32]. The PCR3BP is used with hybrid genetic algorithms and deterministic simplex methods to solve the boundary-value problem andobtain the data in [23]. The PCR3BP is also used with a computer software called AUTO to obtain the datain [48]. The same software is used to obtain optimal values for the same transfer based on the co-circular,co-planar restricted four-body problem in [49]. Topputo‘s data are obtained from Pareto optimal solutions ast 2-impulse non-Keplerian orbit-transfer using the Theory of Functional Connections (cid:63) αβγ ° ° ° ° ° ° ° T ( days ) A ng l e s α , β , and γ Δ V ( km / s ) Δ V m i n ( k m / s ) αβγ ° ° ° ° T ( days ) A ng l e s α , β , and γ Δ V ( km / s ) Δ V m i n ( k m / s ) Fig. 7
The clockwise case: dependencies among all of the parameters in 4BP transfer: the angles α (in red), β (in blue),and γ (in green) associated with the trajectory with the minimum ∆V found (in gray) as functions of the time of flight T . using the direct transcription and shooting strategy through the bi-circular restricted four-body problem[45]. The gray dots shown in Figs. 6 and 7 form the same Pareto optimal solutions given near the point (i)in [45].Note from Table 4 that the TFC counter-clockwise for the PCR3BP is 4 .
64 m/s better than the bestcase available in the literature. In this case, the velocity of the spacecraft at the point B before the impulseis { . , − . } T m/s while the velocity after the impulse - according to section 2.5 - is V Bf = { . , − . } T m/s. The angle between these two velocities vectors is less than 10 − rad. Note from section 2.5 that the velocity V Bf is tangential to the orbit, because the final orbit is circular in the frameof reference fixed in the Moon. Hence, the velocity at point B given by TFC (before the impulse) is alsotangential to the orbit. In fact, the obtained solutions for the minimum ∆V tend to converge to initial andfinal tangential orbits, and this result agrees with [46]. Note that the restriction that the initial and final ° ° ° ° T ( days ) γ Fig. 8
Values of γ that minimizes the costs ∆V for the counter-clockwise transfer. Earth MoonSun Rangeof γ Rangeof γ Fig. 9
Range of the relative position of the Sun that minimizes the costs ∆V with respect to the Earth and Moon fortransfers from 1 to 7 days. Variables associated Counter-clockwise Clockwise ∆V (m / s) 3944 .
83 3949 . (cid:107) δ V A (cid:107) (m / s) 3134 .
41 3137 . (cid:107) δ V B (cid:107) (m / s) 810 .
421 812 . T (days) 4 .
625 4 . α (rad) 4 . . α . ◦ . ◦ β (rad) 4 . . β . ◦ . ◦ γ (rad) 1 . . γ . ◦ . ◦ P E (m) 3 . × − . × − V E (m / s) 2 . × − . × − ˙ x | t = t (m / s) 9799 . . y | t = t (m / s) − . − . Table 3
Values associated to the minimum ∆V found for the transfer using the bi-circular bi-planar R4BP. velocities must be tangential to the initial and final orbits, respectively, is present in the methods to evaluatethe ∆V in all the literature references shown in Table 4. In contrast, the TFC method does not require thisassumption. ast 2-impulse non-Keplerian orbit-transfer using the Theory of Functional Connections (cid:63) ° ° ° ° ° ° ° ° T ( days ) α Δ V ( km / s ) Fig. 10
The cost ∆V as function of α and the total time of flight T . The yellow dot represents the pair α and T for theminimum ∆V . ° ° ° ° °
01 2 3 4 5 6 7135 ° ° ° ° ° T ( days ) β Δ V ( km / s ) Fig. 11
The counter-clockwise case: the cost ∆V as function of β and the total time of flight T . The yellow dot representsthe pair α and T for the minimum ∆V .6 Allan K. de Almeida Junior et al. ° ° ° ° ° ° ° ° T ( days ) β Δ V ( km / s ) Fig. 12
The clockwise case: the cost ∆V as function of β and the total time of flight T . The yellow dot represents the pair α and T for the minimum ∆V . Model Literature reference T (days) ∆V (m/s)PCR3BP [23] 3 . . . . . . . . . . . . . . .
625 3944 . Table 4
The best ∆V found in the literature for a short transfer time (up to seven days) from the Earth to the Mooncompared with the results for the best counter-clockwise and clockwise cases found in this paper. The altitudes of the orbitsaround the Earth and the Moon are 167 km and 100 km, respectively. This study takes advantage of TFC’s ability to embed constraints into functionals and uses them to solveTwo-Point Boundary Value Problems (TPBVP) appearing in perturbed orbit transfer problems. The TPBVPexample selected is the Earth-to-Moon, two-impulse orbit transfer, where the third body influence movesfrom being a small perturbation to the dominant force. The TFC functionals analytically satisfy the boundaryconditions and turn the problem into an unconstrained one. This allows us to solve the problem by nonlinearleast-squares with automatic differentiation and a just-in-time (JIT) compiler, rather than shooting methodsadopted by competing approaches. The initial guess is a spiral trajectory, and it allows us to obtain fastconvergence with residual error on the order of 10 − . This fast procedure can be used to quickly buildnon-Keplerian pork-chop plots that can be used for analysis and/or preliminary mission design.The dynamics of this TPBVP has, in general, multiple solutions that are associated with local minima. This means that, for assigned boundary conditions, there might be several different solutions satisfying theequations of motion. Therefore, the convergence to the desired minimum ∆V tot solution depends on the initialguess. On the other hand, since the convergence is almost always obtained, the method can be used to findall sub-optimal solutions. In general, no more than five distinct solutions for assigned boundary conditions ast 2-impulse non-Keplerian orbit-transfer using the Theory of Functional Connections (cid:63) are found. The longer the transfer time, the more solutions there are. However, when all the solutions havebeen identified, the selection of the optimal is trivial.This fast approach allowed us to analyze the behavior of the ∆V tot as a function of all four mainparameters: the time of flight, the relative Sun direction, and the locations where the two impulses areapplied. The analysis is performed to investigate the total costs around the optimal value of the parametersto quantify the effective contribution. This way, by discretization, pork-chops plots can be built for anyparameter pair. Summarizing, the TFC method to solve this astrodynamics TPBVP allowed us: – to quickly find all distinct solutions; – to build perturbed trajectory pork-chops to understand the cost variations as a function of the parameters’variations.Specifically, two Earth-to-Moon trajectories, with transfer times up to seven days were analyzed for twotarget orbits around the Moon, one clockwise and one counter clockwise, respectively. The transfer minimum ∆V tot cost obtained in the bi-planar, bi-circular 4BP case coincides with the best values found in literature.In the PCR3BP case, the proposed TFC-based approach provides a slightly better solution than those foundin literature. Acknowledgements
This work was supported by FAPESP - S˜ao Paulo Research Foundation through grants 2019/18480-5,2018/07377-6 and 2016/24561-0, and by the NASA Space Technology Research Fellowship, Leake [NSTRF 2019] Grant
References
1. Nima Assadian and Seid H. Pourtakdoust. Multiobjective genetic optimization of earth–moon trajectories in the re-stricted four-body problem.
Advances in Space Research , 45(3):398 – 409, 2010.2. Edward A. Belbruno and James K. Miller. Sun-perturbed earth-to-moon transfers with ballistic capture.
Journal ofGuidance, Control, and Dynamics , 16(4):770–775, 1993.3. James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, and SkyeWanderman-Milne. JAX: composable transformations of Python+NumPy programs, 2018.4. C. Circi and P. Teofilatto. On the dynamics of weak stability boundary lunar transfers.
Celestial Mechanics andDynamical Astronomy , 79:41–72, 2001.5. S. Da Silva Fernandes and C.M.P. Marinho. Sun influence on two-impulsive Earth-to-Moon transfers. In
Proceedings ofthe 22nd International Symposium on Space Flight Dynamics , feb 2011.6. Sandro da Silva Fernandes and Cleverson Maranh˜ao Porto Marinho. Optimal two-impulse trajectories with moderateflight time for earth-moon missions.
Mathematical problems in engineering , 2012(971983), 2012.7. Florin Diacu. Poincar´e and the three-body problem. by june barrow-green.
Historia Mathematica , 26(2):175 – 178,1999.8. Luiz Arthur Gagg Filho and Sandro da Silva Fernandes. A method based on jacobi integral variational equation forcomputing earth-moon trajectories in the four-body problem.
Acta Astronautica , 165:312 – 330, 2019.9. Roy Frostig, Matthew Johnson, and Chris Leary. Compiling machine learning programs via high-level tracing. In
SysMLConference , 2018.10. Roberto Furfaro and Daniele Mortari. Least-squares Solution of a Class of Optimal Space Guidance Problems via Theoryof Connections.
ACTA Astronautica , 2019.11. M. Guelman. Earth-to-moon transfer with a limited power engine.
Journal of Guidance, Control, and Dynamics ,18(5):1133–1138, 1995.12. Albert L. Herman and Bruce A. Conway. Optimal, low-thrust, earth-moon orbit transfer.
Journal of Guidance, Control,and Dynamics , 21(1):141–147, 1998.13. Philip Holmes. Poincar´e, celestial mechanics, dynamical-systems theory and “chaos”.
Physics Reports , 193(3):137 –163, 1990.14. Hunter Johnston, Carl Leake, Yalchin Efendiev, and Daniele Mortari. Selected applications of the theory of connections:A technique for analytical constraint embedding.
Mathematics , 7(6), 2019.15. Hunter Johnston and Daniele Mortari. Least-squares solutions of boundary-value problems in hybrid systems, 2019.16. Hunter Johnston, Enrico Schiassi, Roberto Furfaro, and Daniele Mortari. Fuel-efficient powered descent guidance onlarge planetary bodies via theory of functional connections, 2020.17. C. Lanczos.
Applied Analysis , page 504. Dover Publications, Inc., New York, 1957.18. Carl Leake, Hunter Johnston, and Daniele Mortari. The multivariate theory of functional connections: Theory, proofs,and application in partial differential equations, 2020.19. Carl Leake, Hunter Johnston, Lidia Smith, and Daniele Mortari. Analytically embedding differential equation constraintsinto least squares support vector machines using the theory of functional connections.
Machine Learning and KnowledgeExtraction , 1(4):1058–1083, 2019.20. Carl Leake and Daniele Mortari. Deep Theory of Functional Connections: A New Method for Estimating the Solutionsof Partial Differential Equations.
Machine Learning and Knowledge Extraction , 2(1):37–55, 2020.21. Hanlun Lei, Bo Xu, and Yisui Sun. Earth–moon low energy trajectory optimization in the real system.
Advances inSpace Research , 51(5):917 – 929, 2013.22. Tina Mai and Daniele Mortari. Theory of Functional Connections Applied to Nonlinear Programming under EqualityConstraints, 2019. Paper AAS 19-675 of the 2019 AAS/AIAA Astrodynamics Specialist Conference, Portland, ME,August 11-15, 2019.8 Allan K. de Almeida Junior et al.23. Giovanni Mengali and Alessandro A. Quarta. Optimization of biimpulsive trajectories in the earth-moon restrictedthree-body system.
Journal of Guidance, Control, and Dynamics , 28(2):209–216, 2005.24. G. Mingotti and F. Topputo. Ways to the Moon: a survey. In
Paper AAS 11–283, 21th AAS/AIAA Space FlightMechanics Meeting , feb 2011.25. Giorgio Mingotti, Francesco Topputo, and Franco Bernelli-Zazzera. Efficient invariant-manifold, low-thrust planar tra-jectories to the moon.
Communications in Nonlinear Science and Numerical Simulation , 17(2):817–831, 2012.26. Ashley Moore, Sina Ober-Bl¨obaum, and Jerrold E. Marsden. Trajectory design combining invariant manifolds withdiscrete mechanics and optimal control.
Journal of Guidance, Control, and Dynamics , 35(5):1507–1525, 2012.27. Daniele Mortari. Least-Squares Solution of Linear Differential Equations.
MDPI Mathematics , 5(4), 2017.28. Daniele Mortari. The Theory of Connections: Connecting Points.
MDPI Mathematics , 5(4), 2017.29. Daniele Mortari, Hunter Johnston, and Lidia Smith. High Accuracy Least-squares Solutions of Nonlinear DifferentialEquations.
Journal of Computational and Applied Mathematics , 352:293 – 307, 2019.30. Daniele Mortari and Carl Leake. The Multivariate Theory of Connections.
MDPI Mathematics , 7(3), 2019.31. Kaori Onozaki, Hiroaki Yoshimura, and Shane D. Ross. Tube dynamics and low energy Earth–Moon transfers in the4-body system.
Advances in Space Research , 2017.32. K. Oshima, F. Topputo, and T. Yanao. Low energy transfers to the moon with long transfer time.
Celestial Mechanicsand Dynamical Astronomy , 131(4):0, 2019.33. Kenta Oshima, Francesco Topputo, Stefano Campagnola, and Tomohiro Yanao. Analysis of medium-energy transfers tothe Moon.
Celestial Mechanics and Dynamical Astronomy , 2017.34. L. Peng, Y. Wang, G. Dai, Y. Chang, and F. Chen. Optimization of the earth-moon low energy transfer with differentialevolution based on uniform design. In
IEEE Congress on Evolutionary Computation , pages 1–8, 2010.35. Daniel P´erez-Palau and Richard Epenoy. Fuel optimization for low-thrust Earth-Moon transfer via indirect optimalcontrol.
Celestial Mechanics and Dynamical Astronomy , 130(2):21, February 2018.36. H.J. Pernicka, D.P. Scarberry, S.M. Marsh, and T.H. Sweetser. A search for low delta-v earth-to-moon trajectories.
TheJournal of the Astronautical Sciences , 42:77, 1995.37. E. Perozzi and A.D. Salvo. Novel spaceways for reaching the moon: an assessment for exploration.
Celestial Mechanicsand Dynamical Astronomy , 102:207–218, 2008.38. Henri Poincar´e. Sur les ´equations de la dynamique et le probleme des trois corps.
Acta Math , 13(1):270, 1890.39. AFBA Prado and R. Broucke. Transfer orbits in restricted problem.
Journal of Guidance, Control, and Dynamics ,18(3):593, 1995.40. Yi Qi and Shijie Xu. Optimal earth-moon transfers using lunar gravity assist in the restricted four-body problem.
ActaAstronautica , 134:106 – 120, 2017.41. Carles Sim´o, Gerard G´omez, `Angel Jorba, and Josep Masdemont. The bicircular model near the triangular librationpoints of the rtbp. In
From Newton to chaos , pages 343–370. Springer, 1995.42. Theodore H Sweetser. An estimate of the global minimum dv needed for earth-moon transfer. In sfm , pages 111–120,1991.43. Keith R. Symon.
Mechanics . Addison-Wesley Inc., 1 edition, 1953.44. F. Topputo, M. Vasile, and F. Bernelli-Zazzera. Earth-to-moon low energy transfers targeting l1 hyperbolic transitorbits.
Annals of the New York Academy of Sciences , 1065:55–76, 2005.45. Francesco Topputo. On optimal two-impulse earth–moon transfers in a four-body model.
Celestial Mechanics andDynamical Astronomy , 117(3):279–313, 2013.46. B. F. Villac and D. J. Scheeres. Escaping trajectories in the hill three-body problem and applications.
Journal ofGuidance, Control, and Dynamics , 26(2):224–232, 2003.47. K Wright. Chebyshev Collocation Methods for Ordinary Differential Equations.
The Computer Journal , 6(1):358–365,1964. Issue 4.48. Kazuyuki Yagasaki. Computation of low energy earth-to-moon transfers with moderate flight time.
Physica D: NonlinearPhenomena , 197(3):313 – 331, 2004.49. Kazuyuki Yagasaki. Sun-perturbed earth-to-moon transfers with low energy and moderate flight time.
Celestial Me-chanics and Dynamical Astronomy , 90:197, 2004.50. Hiroshi Yamakawa, Jun’ichiro Kawaguchi, Nobuaki Ishii, and Hiroki Matsuo. A numerical study of gravitational captureorbit in the earth-moon system. In