Generating Solar Sail Trajectories in the Earth-Moon System Using Augmented Finite-Difference Methods

Abstract

Using a solar sail, a spacecraft orbit can be offset from a central body such that the orbital plane is displaced from the gravitational center. Such a trajectory might be desirable for a single-spacecraft relay to support communications with an outpost at the lunar south pole. Although trajectory design within the context of the Earth-Moon restricted problem is advantageous for this problem, it is difficult to envision the design space for offset orbits. Numerical techniques to solve boundary value problems can be employed to understand this challenging dynamical regime. Numerical finite-difference schemes are simple to understand and implement. Two augmented finite-difference methods (FDMs) are developed and compared to a Hermite-Simpson collocation scheme. With 101 evenly spaced nodes, solutions from the FDM are locally accurate to within 1740 km. Other methods, such as collocation, offer more accurate solutions, but these gains are mitigated when solutions resulting from simple models are migrated to higher-fidelity models. The primary purpose of using a simple, lower-fidelity, augmented finite-difference method is to quickly and easily generate accurate trajectories.

1. Introduction

When a permanent outpost on the Moon to support extended human expeditions is eventually established, the astronauts at the facility will require a continual communications link with the Earth. A leading candidate location for the lunar base is at the south pole, which is not always in view of either antennas on the Earth or space-based transceivers, such as the TDRSS satellites, in geosynchronous orbit. Therefore, at least two spacecraft in traditional polar orbits about the Moon are required for an Earth-Moon link [1].

In contrast to a constellation of multiple spacecraft, alternative communications strategies that rely on only one satellite do exist, using current or near-future technology. Advanced propulsion concepts, such as low-thrust ion engines, as well as solar sails, supply a force in addition to gravity and can actually offset an orbit from the Moon [2–4], that is, the orbit plane is displaced from the gravitational center and the spacecraft appears to hover above the surface. Such an orbit allows one spacecraft to be in view of a location on the lunar surface at all times. This dynamical system is complex due to the gravitational effects of the two primaries and, in the case of solar sails, the fact that the Sun, which influences the direction and magnitude of the resulting sail thrust vector, continually moves relative to a fixed Earth-Moon system. A previous approach, based on an understanding of the dynamical structure and using that knowledge to design an orbit, has been successful in developing some solar sail orbits in the vicinity of the lunar Lagrange points [5–9]. However, motion below one of the lunar poles requires an alternative strategy.

Trajectories can be represented as solutions to boundary value problems (BVPs) in terms of ordinary differential equations (ODEs). For a BVP, conditions are specified at the beginning and, also, at the end of the time domain. The specific values of the states at the extremes may or may not be fixed in a BVP with a periodicity constraint; periodicity only requires that the values at the extremes are equal. Differential-corrections schemes [10] and shooting methods (both adaptations of the initial value problem), collocation approaches, and finite-difference methods [11–13] are common numerical processes to solve BVPs and are similar in that they all employ linear corrections to a nonlinear system (these three techniques are discussed in greater detail below). Variational methods are also used to solve BVPs when the problem possesses certain minimality properties [14]. Classical perturbation techniques have long supplied analytical approximations for trajectories [5, 15–17]. Finally, recent efforts using evolutionary algorithms to solving optimization problems also provide an approach to solving boundary value problems [18, 19].

Shooting methods [10, 12, 20–23] and collocation schemes [2, 24–31] possess high degrees of accuracy and are frequently employed for trajectory generation and optimization. Shooting methods typically require some knowledge of the design space to initiate the procedure. Collocation schemes do not require a priori knowledge of the solution corresponding to the BVP, but the implementation is often neither intuitive nor straightforward. Furthermore, the partial derivatives of the constraints with respect to the states are nontrivial in collocation schemes and are often approximated numerically [2].

A more rudimentary method for solving BVPs is the finite-difference method (FDM), in which the derivatives that appear in the differential equations are replaced with their respective finite differences and evaluated at node points along the trajectory [12]. The solution process is iterative. The trajectory is discretized, and the equations that represent the relationships at the nodes are solved simultaneously. Using a central difference approximation, the smallest local errors associated with an FDM approach are proportional to Δ𝑡². However, model accuracy and uncertainty may obviate any precision gains from a high-precision scheme; solutions from those methods may not persevere without additional control in a higher-fidelity model (e.g., a model that includes planetary ephemerides) or in actual flight. If a solution based on a high-precision approach (e.g., shooting or high-degree collocation) is migrated from a model with limited accuracy to one of higher accuracy (e.g., circular restricted three-body model versus ephemeris model or actual flight), it may be no more accurate in that high-fidelity model than a solution from a low-precision approach (e.g., finite-differencing). Moreover, the FDM is simple to formulate and offers an improved understanding of the design space. Partial derivatives are easily obtained via standard analytical techniques, without having to resort to numerical or automatic differentiation (each is often computationally expensive). Once the analyst is familiar with the feasible options, higher-fidelity methods for solving the BVP can be used to refine trajectories or combined with optimization techniques to produce optimal trajectories. Another option is to employ the FDM technique to generate trajectories that fit path constraints to quickly explore a broad design space. Using the MATLAB computing environment, Wawrzyniak and Howell [32] survey over 10 million combinations of initial guesses for the path and control profile, as well as a range of characteristic accelerations and path constraints appropriate to the lunar south pole coverage problem. The time to generate the millions of trajectories for the survey is approximately one week when the algorithms run on eight cores over five platforms.

In the following sections, an FDM is described and applied to solar sail orbits in the Earth-Moon circular restricted three-body (CR3B) system. Also included is a modified FDM for the same application, where the velocity is incorporated as part of the solution at each node along the discretized trajectory. Neither method is strictly a straightforward, “textbook” FDM [11–13], since both are augmented with trajectory and control constraints. Additionally, in contrast to a simple two-point BVP, where the states are fixed at the extremes, the trajectory is required to be periodic; nowhere along the arc is the solution prespecified.

A description of the dynamical model in the Earth-Moon CR3B system is followed by an algorithm using the augmented FDMs for generating trajectories. The strategies are based on minimizing the difference between accelerations from the equations of motion (as well as the velocities as another alternative) and the corresponding values numerically derived from positions along a discretized trajectory. An error analysis is also included. A separate study uses these FDMs for surveying the solution space and assessing the required solar sail and spacecraft characteristics necessary for the lunar south pole (LSP) coverage problem [32].

2. System Dynamical Model

The problem that represents this application is defined within the context of the CR3B system, that is, the problem is formulated in a frame, 𝑅, that is rotating with respect to an inertial system, 𝐼. A CR3B model that incorporates the gravity contributions of two primary bodies is geometrically advantageous for understanding the problem. Consistent with McInnes [33], the nondimensional vector equation of motion for a spacecraft at a location 𝐫 relative to the barycenter (center of mass of the primaries) is𝐚+2^𝑅𝐼𝝎^𝑅×^𝑅𝐯+∇𝑈(𝐫)=𝐚_𝑠(𝑡),(1) where the first term is the acceleration relative to the rotating frame (more precisely expressed as ^𝑅𝑑²𝐫/𝑑𝑡², where the left superscript 𝑅 indicates a derivative in the rotating frame) and the second term is the corresponding Coriolis acceleration, which requires the velocity relative to the rotating frame, ^𝑅𝐯 (more precisely ^𝑅𝑑𝐫/𝑑𝑡). Vectors are denoted with boldface. Derivatives of the position vector, ^𝑅𝐯 and ^𝑅𝐚, are assumed to be relative to the rotating frame and, consequently, 𝑅 is dropped. The angular-velocity vector, ^𝐼𝝎^𝑅 (or 𝝎), relates the rate of change the rotating frame with respect to the inertial frame. The applied acceleration, from a solar sail in this case, is indicated on the right side by 𝐚_𝑠(𝑡). The pseudogravity gradient, ∇𝑈(𝐫), combines centripetal and gravitational accelerations:∇𝑈(𝐫)=(𝝎×(𝝎×𝐫))+(1−𝜇)𝑟³¹𝐫_𝟏+𝜇𝑟³²𝐫_𝟐,(2) where 𝜇 represents the mass fraction of the smaller body, or 𝑚₂/(𝑚₁+𝑚₂), and 𝑟₁ and 𝑟₂ are the distances from the larger and smaller bodies, respectively, that is, 𝑟₂𝑟₁=(𝜇+𝑥)²+𝑦²+𝑧²,=(𝜇+𝑥−1)²+𝑦²+𝑧².(3) Solar gravity is neglected in this model. At a distance of 1 AU, an appropriate distance to assume for a sailcraft in the Earth-Moon system, the applied acceleration from a solar sail is modeled as𝐚_𝑠(𝑡)=𝛽𝓘(𝑡)⋅𝐮̂²𝐮̂,(4) where 𝐮̂ is the sail-face normal, 𝓘(𝑡) is a unit vector in the sun-to-spacecraft direction, and 𝛽 is the sail characteristic acceleration in nondimensional units. These vectors appear in Figure 1. Observed from the rotating frame, 𝑅, the Sun moves in a clockwise direction about the fixed primaries. The sail mass, 𝑚₃, is negligible compared to the masses of the Earth and Moon, which are 𝑚₁ and 𝑚₂, respectively. The term (𝓘(𝑡)⋅𝐮̂) is also expressed as cos𝛼, where 𝛼 is the sail pitch angle, or the angle between the solar incidence direction and the sail normal.

Figure 1 

Earth-Moon system model.

To generate the magnitude of the sail acceleration in dimensional units, 𝑎₀, 𝛽 is multiplied by the system characteristic acceleration, 𝑎^∗, which is the relationship between the dimensional and nondimensional acceleration in (1). In fact, 𝑎^∗ is the ratio of the characteristic length, 𝐿^∗ (384,400 km for the Earth-Moon distance), to the square of the characteristic time, 𝑡^∗ (2𝜋𝑡^∗=27.321 days), that is, 𝑎^∗=𝐿^∗𝑡^∗2=2.7307mm/s².(5) A recent sailcraft design for NASA's Space Technology competition (ST9) that was built by L'Garde possesses overall characteristic acceleration, 𝑎₀, of 0.58mm/s², while the characteristic acceleration of the sail and its support structure alone is closer to 1.70mm/s² (0.212 to 0.623 in units of nondimensional acceleration, resp.) [34].

The sunlight direction is expressed relative to the rotating frame and is a function of time, that is,𝓘(𝑡)=cos(Ω𝑡)𝐱̂−sin(Ω𝑡)𝐲̂+0𝐳̂,(6) where Ω is the ratio of the synodic rate of the Sun as it moves along its path to the system rate, approximately 0.9192. One physical constraint is imposed on the attitude of the spacecraft: the sail normal, 𝐮̂, which is coincident with the direction of the resultant force in an ideal model, is always directed away from the Sun. This is written mathematically as𝓘(𝑡)⋅𝐮̂≥cos𝛼_max,(7) where 𝛼_max is 90^∘.

The sail modeled here is a perfectly reflecting, flat solar sail. Billowing is not incorporated in this force model; however, 𝛼_max could be less than 90^∘, as sail luffing (i.e., flapping) is assumed to occur at high pitch angles [35]. Higher fidelity models include optical models [33], parametric models that incorporate billowing in addition to optical effects [33, 36], and realistic models based on finite-element analysis that incorporates optical properties and manufacturing flaws [37]. Optical effects can represent a nonperfectly reflecting solar sail; some energy is absorbed, and some is reflected diffusely as well as specularly. An ideal sail reflects only specularly. In all of these models, the resulting acceleration from a solar sail is not perfectly parallel to the sail-face normal but, instead, is increasingly offset from the sail-face normal as the sail is pitched further from the sunlight direction [33]. Fully accounting for realistic solar sail properties attenuates the sail characteristic acceleration by nearly 25% and places an upper limit on the pitch angle between 50^∘ and 60^∘, depending on the properties of the sail [37]. Nevertheless, this analysis will employ an ideal sail to lend insight into the technology level that is required to solve the LSP coverage problem and into the use of augmented finite-difference methods for generating solar sail trajectories.

3. Augmented Finite-Difference Methods

The derivation for an augmented finite-difference method begins with generic, nonlinear, second-order, two-point BVP that can be represented as𝐫̈=𝐟(𝑡,𝐫,𝐫̇),(8) where 𝐫(𝑡₁) and 𝐫(𝑡_𝑛) are fixed at the extremes of the time span, [𝑡₁,𝑡_𝑛] in interval notation [11–14]. The solution of (8) is a trajectory, 𝐫(𝑡). To solve the BVP, the classic FDM discretizes the domain into nodes, or epochs, 𝑡₁,𝑡₂,…,𝑡_𝑛 and replaces the derivatives in (8) by their respective finite differences. Central-difference approximations (CDAs) are commonly used because they are more accurate than forward or backward differences. The first and second time derivatives of 𝐫(𝑡) are, of course, the velocity and acceleration vectors and are approximated by CDAs as𝐯_𝑖̃=𝐫^𝑖+1−𝐫_𝑖−1𝑡_𝑖+1−𝑡_𝑖−1𝐚̃,(9)_𝑖=2𝐫^𝑖+1−𝐫_𝑖𝑡_𝑖−𝑡_𝑖−1−𝐫_𝑖−𝐫_𝑖−1𝑡_𝑖+1−𝑡_𝑖𝑡_𝑖−𝑡_𝑖−1𝑡_𝑖+1−𝑡_𝑖−1𝑡_𝑖+1−𝑡_𝑖,(10) where 𝐫(𝑡_𝑖) is defined as 𝐫_𝑖 and 𝑡_𝑖∈(𝑡_𝑖−1,𝑡_𝑖+1) and the symbol “⋅” indicates a numerical approximation. Equation (10) arises from three central-differences: one for the velocity at 𝑡_𝑖−1/2, which is midway between 𝑡_𝑖−1 and 𝑡_𝑖, another at 𝑡_𝑖+1/2, midway between 𝑡_𝑖 and 𝑡_𝑖+1, and a third using the first two intermediate velocities resulting in the acceleration at 𝑡_𝑖. The FDM does not require uniform node placement (mesh refinement and nonuniform node spacing may improve the accuracy of the FDM), but when the time steps between 𝑡_𝑖−1, 𝑡_𝑖, and 𝑡_𝑖+1 are fixed at a value of Δ𝑡, then the midpoint of [𝑡_𝑖−1,𝑡_𝑖+1] is 𝑡_𝑖 and (9) and (10) simplify to𝐯_𝑖̃=𝐫^𝑖+1−𝐫_𝑖−1𝐚̃2Δ𝑡,(11)_𝑖=𝐫^𝑖+1−2𝐫_𝑖+𝐫_𝑖−1Δ𝑡².(12) Note that 𝐫_𝑛−1 is used for 𝐫₀ when calculating 𝐯₁̃ for a periodic solution. Similarly, 𝐫₁ substitutes for 𝐫_𝑛 in 𝐯̃_𝑛−1, so the problem does not include redundant constraints. The calculations for 𝐚₁ and 𝐚_𝑛−1 are similar. The differences between the actual velocity and acceleration and their respective CDAs at 𝑡_𝑖 are𝐯_𝑖−𝐯_𝑖̃=−Δ𝑡²𝐚6𝐫^𝜌_𝑣,(13)_𝑖−𝐚̃_𝑖=−Δ𝑡²12𝐫^iv𝜌_𝑎,(14) where 𝜌_𝑣,𝑎∈(𝑡_𝑖−1,𝑡_𝑖+1) and 𝐫^(𝑡) and 𝐫^iv(𝑡) are continuous on [𝑡_𝑖−1,𝑡_𝑖+1]. Derivations of these truncation errors are available in the literature [11–13]. The approximations in (11) and (12) replace the first and second time derivatives in (8), and the result is an equation at each epoch of the form𝐫_𝑖−1+2𝐫_𝑖+𝐫_𝑖+1Δ𝑡²=𝐟𝑡,𝐫,𝐫^𝑖+1−𝐫_𝑖−12Δ𝑡.(15) Equation (15) may still be nonlinear, but it can be linearized into a system of equations, one for each 𝑡_𝑖, excluding the extremes, and solved iteratively. Due to the approximations in (11) and (12), the classic FDM results in a solution that will approximate each component of the true solution to the BVP by an error proportional to Δ𝑡² at each node [11, 13].

For the current problem, the FDM is augmented to incorporate path and system constraints, as well as a control history. Two variations of this FDM are developed. In the first, the position is the only explicitly discretized trajectory state (FDM-R). The second FDM is formulated to directly solve for position and velocity along the discretized trajectory (FDM-RV).

The FDM-R formulation is based on approximating (1) via the CDAs for velocity and acceleration given by (11) and (12), respectively, that is,𝐚̃_𝑖+2𝝎×𝐯_𝑖̃+∇𝑈𝐫_𝑖−𝐚_𝑠𝑡_𝑖=0.(16) Note that the velocity term in (16) is replaced with its CDA, 𝐯_𝑖̃. To solve the BVP, the trajectory is discretized at 𝑡₁,𝑡₂…,𝑡_𝑛. Position is expressed in terms of Cartesian coordinates at each node: 𝐫_𝑖=𝑥_𝑖𝑦_𝑖𝑧_𝑖^T,(17) where “T” indicates a vector transpose. The control, 𝐮_𝑖, that is, the sail pointing vector in this application, as well as a set of 𝑚 slack variables, 𝜂_𝑖, are also included at each node. Note that the control, 𝐮_𝑖, is not necessarily a unit vector until the FDM converges on a solution. The positions, control, and slack variables at node 𝑖 for the FDM-R are collected into the subvector, 𝐪_𝑖=𝐫^T_𝑖𝐮^T_𝑖𝜼^T_𝑖^T_(6+𝑚),(18) where the length of 𝐪_𝑖 is (6+𝑚). The subvectors corresponding to each node are collected into a column vector,𝐗=𝐪^T₁𝐪^T₂⋯𝐪^T_𝑛−1𝐪^T_𝑛^T_(6+𝑚)𝑛,(19) which represents the discretized trajectory as well as the associated control history and slack variables. The velocities associated with the trajectories generated by the FDM-R can be reconstructed with CDAs of 𝐫(𝑡) during postprocessing.

The FDM-RV process is similarly formulated, but there are two distinct sets of ODEs to solve, that is,𝐚̃𝐯_𝑖̃−𝐯_𝑖=0,(20)_𝑖+2𝝎×𝐯_𝑖+∇𝑈𝐫_𝑖−𝐚_𝑠𝑡_𝑖=0.(21) Note the difference in the velocity terms in (16) and (21). With this FDM-RV alternative, the subvector is𝐪_𝑖=𝐫^T_𝑖𝐯^T_𝑖𝐮^T_𝑖𝜼^T_𝑖^T_(9+𝑚),(22) since 𝐯_𝑖 is now explicitly part of the solution. This formulation results in an 𝐗 vector with length (9+𝑚)𝑛.

Both the FDM-R and the FDM-RV algorithms constrain the difference between the evaluated EOMs at node 𝑖 and the associated numerically derived approximations to be zero. Additional constraints are required for periodicity, control-direction magnitude, and path characteristics. All of these constraints are dependent on 𝐗 and are collected in the column vector 𝐅(𝐗). An expansion of 𝐅(𝐗) about 𝐗^𝑗 yields a linear approximation for 𝐅(𝐗^𝑗+1), that is, 𝐅𝐗^𝑗+1=𝐅𝐗^𝑗+𝐷𝐅𝐗^𝑗𝐗^𝑗+1−𝐗^𝑗+𝒪Δ𝐗²,(23) where the 𝐷𝐅(𝐗^𝑗) is the Jacobian, 𝜕𝐅(𝐗^𝑗)/𝜕𝐗^𝑗. Superscripts on vectors refer to iteration number. The initial guess is denoted 𝑗=0; a converged solution by 𝑗=𝑓. Partial derivatives used in the Jacobian are derived analytically, but are not specifically presented here. By assuming that 𝐗^𝑗 is in the neighborhood of 𝐅(𝐗)=0 (i.e., ‖𝐗^𝑗+1−𝐗^𝑗‖≪1), (23) is rearranged into a least-norm problem, that is,𝐗^𝑗+1=𝐗^𝑗−𝐷𝐅𝐗^𝑗^T𝐷𝐅𝐗^𝑗⋅𝐷𝐅𝐗^𝑗^T⁻¹𝐅𝐗^𝑗.(24) Equations (23) and (24) reflect the Newton-Raphson method of root finding for several variables. For efficient computation, 𝐷𝐅(𝐗) is preallocated and stored in memory as a sparse matrix [38]. Iterations based on (24) continue until some convergence tolerance is met, that is,‖‖𝐗^𝑗+1−𝐗^𝑗‖‖‖‖𝐗^𝑗‖‖≤tol.(25) When an initial guess is in the neighborhood of a solution, convergence is quadratic. Specifying tol≈1×10⁻⁷ is sufficient, since any further iteration approaches double precision. If the initial guess is not in the neighborhood of a solution, 𝐗^𝑗+1 bears little resemblance to 𝐗^𝑗; the step from 𝐗^𝑗 to 𝐗^𝑗+1 may even be chaotic. However, the new 𝐗^𝑗+1 might be in the neighborhood of an alternate solution and subsequently lead to convergence. The converged discretized solution, 𝐗^𝑓, that satisfies these constraints is the trajectory that solves the EOMs to within a theoretical error.

Algebraically, these FDMs are equivalent. In practice, the two formulations yield slightly different results and, sometimes, dramatically different results. This problem is sensitive to the initial guess; because the initial condition may or may not explicitly include velocity, the vector of initial guesses is not the same for the two FDMs. Additionally, the least-norm update in (24) minimizes the norm of the correction to the state vector; the state vectors for the two formulations are different, and the least-norm update for an initial guess that includes velocity states may move the solution to a different region when compared to an update that does not contain velocity states between iterations. The convergence criterion (25) also depends on the norm of the update to the state vector. Because the size of the state vectors differs between the methods, the vector norm of the update will also differ for converged solutions. The current formulations do not control specified states for convergence (e.g., position and not velocity), and other variations on the implementation of these methods will yield different results. Given that the objective is a simple and quick method to approximate the design space, neither strategy should be considered superior to the other; both achieve the stated goal and results from both are insightful.

At each epoch, the converged subvector 𝐪_𝑖 differs from the true solution, 𝐪_𝑖, due to the truncation error associated with the FDMs and the limited machine precision. The proceeding error analysis is described by Kincaid and Cheney [13, pages 589–592]. For a step size Δ𝑡 that is greater than some value, the error due to truncation in (9)–(14) will dominate the machine error. Assuming that there exists a true set of position states, 𝐫, that solve (1) and (13) as formulated for the FDM-R approach, (8) is written as𝐫_𝑖−1+2𝐫_𝑖+𝐫_𝑖+1Δ𝑡²−Δ𝑡²12=−𝝎×𝐫^iv𝜌_𝑎𝐫_𝑖+1−𝐫_𝑖−1Δ𝑡+Δ6𝑡²𝝎×−∇𝑈𝐫^𝜌_𝑣𝐫_𝑖+𝐚_𝑠𝑡_𝑖.(26) The pseudogravity gradient associated with the true solution, 𝐫_𝑖, is expanded as a function of an approximate solution, 𝐫_𝑖, that is,∇𝑈𝐫_𝑖=∇𝑈𝐫_𝑖+𝐌𝐫_𝑖𝐞_𝑖+𝒪𝐞^T_𝑖𝐞_𝑖,(27) where 𝐞_𝑖=𝐫_𝑖−𝐫_𝑖 and 𝐌(𝐫_𝑖) is the Hessian of 𝑈. Using the solution from the FDM-R algorithm, (15) is subtracted from (26), resulting inΔ𝑡⁻²=𝐞−Δ_𝑖−1𝑡+₋₁2𝐞_𝑖+𝐞_𝑖+1𝝎×𝐞_𝑖+1−𝐞_𝑖−1−𝐌_𝑖𝐞_𝑖+Δ𝑡²𝐡_𝑖,(28) where 𝐡_𝑖=(1/12)𝐫^iv(𝜌_𝑎)+(1/6)𝐫^(𝜌_𝑣). It is assumed that 𝐫(𝑡) is continuously differentiable to fourth order. Rearranging and multiplying by Δ𝑡² yields(−𝐈+Δ𝑡𝝎×)𝐞_𝑖−1+2𝐈−Δ𝑡²𝐌_𝑖𝐞_𝑖+(−𝐈−Δ𝑡𝝎×)𝐞_𝑖+1=−Δ𝑡⁴𝐡_𝑖,(29) where 𝐈 is the identity matrix and 𝝎× is the skew-symmetric gyroscopic matrix. Some terms in (29) can be simplified to a form that will be useful later, that is(−𝐈+Δ𝑡𝝎×)+2𝐈−Δ𝑡²𝐌_𝑖+(−𝐈−Δ𝑡𝝎×)=−Δ𝑡²𝐌_𝑖.(30) Now, let 𝝀 correspond to the difference element 𝐞_𝑖 that possesses the largest magnitude. Equation (29) reduces to an upper bound on 𝝀,(−𝐈+Δ𝑡𝝎×)𝝀+2𝐈−Δ𝑡²𝐌_𝑖𝝀+(−𝐈−Δ𝑡𝝎×)𝝀≥−Δ𝑡⁴𝐡_𝑖.(31) Using (30), (31) simplifies to−Δ𝑡²𝐌_𝑖𝝀≥−Δ𝑡⁴𝐡_𝑖,(32)𝐌_𝑖𝝀≤Δ𝑡²𝐡_𝑖,(33)𝝀≤Δ𝑡²𝐌^𝑖−1𝐡_𝑖.(34) The inequality in (34) represents three scalar equations. Thus, the position error corresponding to the FDM-R formulation is 𝒪(Δ𝑡²) as Δ𝑡→0. To calculate the velocity from the FDM-R algorithm, the CDA from (11) is used, also with an error 𝒪(Δ𝑡²). The error analysis for the alternate FDM-RV algorithm incorporates the fact that the velocity, 𝐯=𝐟(𝑡,𝐫), possesses errors that can also be expressed by (13), and, therefore, the derivation continues with the same steps as those for the FDM-R formulation beginning with (26).

The augmented finite-difference methods sacrifice precision for simplicity. At each node, the error in the trajectory is proportional to Δ𝑡², significantly less precise than common collocation methods [26, 39]. In this analysis, Δ𝑡=0.067, or 𝑛=101, which translates to an error of 0.452% of the Earth-Moon distance, approximately 1740 km, in each direction. To improve the precision of the FDM, the step size between nodes, Δ𝑡, can be decreased. The Jacobian 𝐷𝐅(𝐗) becomes quite large, imposing a significant computational cost. Furthermore, the smaller the step size, the larger the machine error [12]. Alternatively, Richardson extrapolation, or extrapolation to the limit, could be employed to further minimize the error of the solutions [12]. Incorporating a process to determine the specific step size that minimizes some combination of truncation and machine errors or incorporating an additional process, like Richardson extrapolation, adds complexity to a method that is formulated to be simple to understand and implement, as well as quick to yield results.

Prior to any analysis, the design space for this problem is not well known, so a trajectory that is precise to within Δ𝑡² relative to an actual solution is meaningful. If the goal is a general understanding of the design space, accomplished in a relatively quick analysis, then these results can be very insightful.

4. Algebraic Constraint Vector: F(X)

At the core of the augmented finite-difference methods is the algebraic constraint vector, 𝐅(𝐗), which contains the algebraic approximations to the equations of motion and other constraints necessary for the desired periodic solar sail trajectory. In this formulation, each element in 𝐅(𝐗) should be zero to simultaneously satisfy the ODEs and path constraints. For a periodic solar sail trajectory, subject to path constraints 𝐠(𝐗),𝐅(𝐗)=Δ𝐚(𝐗)^TΔ𝐯(𝐗)^T𝐓(𝐗)^T𝑦₁𝐍(𝐗)^T𝐠(𝐗)^T^T.(35) The presence of Δ𝐯(𝐗) depends on whether the formulation is the FDM-R or FDM-RV, and the length of the 𝐅(𝐗) vector is either (2+(4+𝑚)𝑛) or (2+(7+𝑚)𝑛), accordingly. The order of these elements within the vector is arbitrary and no difference in performance is apparent when 𝐅(𝐗) and 𝐗 are rearranged such that 𝐷𝐅(𝐗) is a sparse, banded-diagonal matrix. Using the configuration in (35), each element set appears as block diagonal in the corresponding Jacobian, 𝐷𝐅(𝐗). Each set is subsequently discussed, and the corresponding size of each set is indicated in a subscript outside of the braces.

The first element set in (35) is the difference in acceleration. The acceleration at a given epoch, 𝑖, as evaluated from the equations of motion, depends only on the position, velocity, and control, 𝐚_𝑖=𝐟(𝐫_𝑖,𝐯_𝑖,𝐮_𝑖). The numerically derived acceleration, 𝐚̃_𝑖 from (10) or (12), depends on the state at multiple epochs. A valid trajectory possesses the same acceleration whether computed from the EOMs or from the CDAs, within the truncation error. Therefore, the difference between accelerations resulting from the EOMs and those from numerical calculation should nominally be zero for a converged solution; thus, the difference forms the first set of equality constraints in the composite constraint vector:Δ𝐚(𝐗)=⎧⎪⎨⎪⎩𝐚¹𝐫₁,𝐯₁,𝐮₁−𝐚̃₁𝐚_𝑛−1⋮(𝐫,𝐯,𝐮)||_𝑛−1−𝐚̃_𝑛−1⎫⎪⎬⎪⎭_3(𝑛−1).(36) Likewise, the difference between the velocity from the EOMs and the numerically derived velocity is zero for a valid trajectory. Velocity from the EOMs at each epoch, 𝐯_𝑖, is available from 𝐗 (see (19) and (22)). This difference forms a second set of equality constraints (only in the FDM-RV algorithm) in 𝐅(𝐗):Δ𝐯(𝐗)=⎧⎪⎨⎪⎩𝐯¹−𝐯₁̃𝐯_𝑛−1⋮−𝐯̃_𝑛−1⎫⎪⎬⎪⎭_3(𝑛−1),(37) where 𝐯_𝑖 may be approximated as 𝐯_𝑖̃. These differences in acceleration and velocity may also be denoted as defects and are located at the node corresponding to each 𝑡_𝑖. In an initial guess for the trajectory state vector, the defects Δ𝐚_𝑖 and Δ𝐯_𝑖 at each node will most likely not be zero. The Newton-Raphson iteration process adjusts the path of the trajectory to resolve these differences.

The next constraints enforce periodicity and unit length of the control vector. To enforce periodicity, the goal is an originating and final state vector that are equal, which is represented as a set of constraints in the algebraic constraint vector, 𝐅(𝐗), that is,𝐓(𝐗)=⎧⎪⎪⎨⎪⎪⎩𝐫^𝑛−𝐫₁𝐮𝐯_𝑛𝑛−𝐯₁−𝐮₁𝜼_𝑛−𝜼₁⎫⎪⎪⎬⎪⎪⎭_(6∣9+𝑚),(38) where 𝑚 is the number of path constraints at the first and last node, expressed as slack variables. The length of 𝐓(𝐗) depends on the application of an FDM-R or an FDM-RV approach. Although the only boundary condition is periodicity, for consistency the position at the boundary is constrained to be in the 𝑥𝑧-plane. Therefore, the 𝑦-coordinate at the first point along the trajectory is constrained, 𝑦₁=0. In practice, without this constraint, 𝑦₁ remained very close to zero, approximately 40 km from the 𝑥𝑧 plane once the solution is converged. With this constraint, 𝑦₁ is reduced to 0.004 mm from the plane. Next, the magnitude of the control vector must remain of unit length in the iteration process𝐍(𝐗)=⎧⎪⎨⎪⎩𝐮^T1𝐮𝐮⋮₁−1_T𝑛−1𝐮_𝑛−1−1⎫⎪⎬⎪⎭_(𝑛−1).(39) Formulating the constraint in terms of 𝐮^T_𝑖𝐮_𝑖−1 as opposed to 𝐮^T_𝑖𝐮_𝑖−1 results in a simpler partial derivative while retaining the intended effect; this constraint is incorporated into the implementation of the FDM since the control history in this analysis is described in terms of a unit vector.

Path constraints are included as inequality constraints, which are converted to equality constraints by use of slack variables, a successful numerical adaptation from nonlinear programming [2, 40]. The slack variables are incorporated into 𝐗 and are associated with the other state elements at their particular epoch, 𝑖. In this analysis, the elevation angle constraint, 𝐸_min, maintains the visibility of the spacecraft from an outpost located near the south pole of the Moon, and the spacecraft altitude constraint, 𝐴_max, is imposed for radio power restrictions. Altitude is defined as the distance from the lunar south pole, that is,𝐴_𝑖=𝑥_𝑖−1+𝜇²+𝑦^2𝑖+𝑧_𝑖+𝑅_𝑚².(40) The third path constraint requires that the sail-face normal, or control 𝐮_𝑖, is always directed away from the Sun (the sunlight vector is 𝓘_𝑖), or 𝛼_max=90^∘ in (7). Of the inequality constraints, only this attitude requirement is mandated. Together, these three inequality constraints are written as 𝐸_min≤𝐸_𝑖≡arcsin−𝑧^𝑖+𝑅_𝑚𝐴𝐴_𝑖,(41a)_maxcos𝛼_max≥𝐴_𝑖,(41b)≤𝓘_𝑖⋅𝐮_𝑖.(41c)For the given problem and model, adding a path constraint to avoid the penumbra and umbra of the Earth or the Moon's shadows is unnecessary because of the elevation constraint. A shadow constraint could be added for another application or shadowing effects could be incorporated into the dynamical model directly [41]. Additional inequality path constraints could include limits on the body turn rates and accelerations governed by the attitude control system of the spacecraft. The associated slack variables are squared and added to the inequality constraints, resulting in the following equality constraints: 𝐠_𝑖𝐫_𝑖,𝜼_𝑖=⎧⎪⎪⎨⎪⎪⎩sin𝐸^min+𝑧^𝑖+𝑅_𝑚𝐴_𝑖+𝜂^2𝐸,𝑖𝐴_𝑖−𝐴_max+𝜂^2𝐴,𝑖cos𝛼_max−𝓘_𝑖⋅𝐮_𝑖+𝜂^2ℓ,𝑖⎫⎪⎪⎬⎪⎪⎭_𝑚.(42) Combining the path constraints at all epochs results in 𝑚(𝑛−1) total elements in 𝐠(𝐗). When using the path constraints in (42), 𝑚=3. The first two path constraints in (42) as well as the periodicity constraint in (38) are illustrated in Figure 2. The cos𝛼_max constraint appears in Figure 1. The sailcraft in Figure 2 orbits below the Moon. Feasible solutions could exist that do not orbit below the Moon, but do orbit below either the 𝐿₁ or the 𝐿₂ point.

Figure 2 

Path constraints for an orbit below the Moon (Moon image from nasa.gov).

Only the first (𝑛−1) epochs are required in formulating Δ𝐚, Δ𝐯, 𝐍, and 𝐠(𝐗) in the constraint vector 𝐅(𝐗). Due to periodicity, the first and 𝑛th epochs are identical. Activating the 𝑛th epoch yields a problem that is overconstrained and the Jacobian, 𝐷𝐅(𝐗), is not full rank.

As mentioned, the Jacobian is a sparse matrix. A diagram of the sparsity pattern based on partial derivatives of 𝐅(𝐗) with respect to 𝐗 for the FDM-RV formulation appears in Figure 3. Diagonals corresponding to the specific element sets in 𝐅(𝐗) are labeled in the figure. For a scenario employing 101 nodes, the size of the Jacobian is 1013 rows by 1212 columns where approximately 0.38% of the elements are nonzero. The sparsity pattern for the FDM-R formulation is similar to that of the FDM-RV, except that the diagonal corresponding to 𝐯(𝐗) does not exist in the FDM-R option. The size of the Jacobian for the FDM-R option associated with 101 nodes is 710 rows by 909 columns; approximately 0.62% of the entries are nonzero.

Figure 3 

Sparsity pattern for the Jacobian 𝐷𝐅(𝐗) in the FDM-RV formulation.

The advantage of a finite-difference approach is its ease of implementation and the speed improvements enabled by its simplicity. Partial derivatives are easily accessible via analytical derivation, especially for an idealized force model such as (1). As complexity of the formulation increases, which is the case with collocation, for example, analysts rely on numerical or automatic differentiation to generate partial derivatives. MATLAB supplies the function NUMJAC.M, which numerically approximates a derivative using a forward difference approximation [42]. Third-party software for MATLAB exists for automatic (a.k.a. algorithmic) differentiation. TOMLAB Optimization's MAD (MATLAB Automatic Differentiation) suite employs automatic differentiation [43], and Shampine's PMAD (Poor Man's Automatic Differentiation) exploits complex step differentiation [44]. Both MAD and PMAD result in highly accurate derivatives, on the order of machine precision when compared to the partial derivatives determined analytically, while derivatives calculated by NUMJAC.M are accurate to 10⁻⁸. Because analytical derivatives are directly encoded into the algorithm, computation times are significantly faster when compared to the other options listed. Table 1 lists the computations times for generating 𝜕∇𝑈/𝜕𝐫, the most computationally intensive step in constructing an entry in the Jacobian, for the four differentiation strategies. Clearly, if analytical derivatives are easily available, they should be employed. Furthermore, if the subroutine to generate the derivatives analytically is programmed in a lower-level language (e.g., C or FORTRAN) and then complied as a MATLAB Executable (MEX), further time savings are realized (with compiled code the computation of 𝜕∇𝑈/𝜕𝐫 takes less than 0.001 seconds). The formulation of the FDMs for this analysis exploits the computational advantage of analytical derivatives and compiled MEX code.

Table 1 

Computation times for various differentiation strategies.

In summary, a finite-difference method yields a solution for the path by replacing the path derivatives with their finite-difference approximations. The differences between the approximate derivatives and the derivatives determined by evaluating the equations of motion at a specific point along a path are minimized by iteratively solving a linearized system of equations. Inequality path constraints are added to this set of equality constraints by way of slack variables, and subsequently augmenting the constraint vector and linearized system of equations, resulting in a feasible path.

5. Results

The objective when employing a finite-difference method is to quickly generate a feasible trajectory where system constraints and sail characteristics are given. Because model simplification into the CR3B system ignores lunar librations, the inclined orbits of the primaries, additional perturbations to the orbit, and lunar surface features, a conservative minimum elevation angle of 𝐸_min=15^∘ is selected. A generous maximum-altitude constraint of one Earth-Moon distance (384,400 km) is also selected (in practice most results are an order of magnitude closer to the lunar south pole and this constraint is not active) to establish the constraints in (42).

As mentioned, a recent sail design possesses a characteristic acceleration of 0.58mm/s² [34]. JAXA's IKAROS is the first mission to have flown using a sail as its only propulsive device, and NASA recently deployed the NanoSail-D2. Both of these sailcraft are relatively small (200 m² and 10 m², resp.) and are designed to demonstrate in-space deployment of the sail. Designs for larger sailcraft exist, and future solar sail spacecraft are likely to be hybridized with other propulsion devices [45]. However, it is instructive to demonstrate the FDMs with characteristic accelerations for sail spacecraft that may be available with near-future technological improvements. An assessment of the technology level (i.e., characteristic acceleration) is required for a sail-only LSP mission. The focus of this analysis is not to demonstrate that a sail-only LSP mission is feasible based on current technology. Instead, the results of this investigation demonstrate the use of finite-difference methods for the sample application of an LSP mission. Therefore, a high characteristic acceleration, that is, 1.70mm/s², is used to generate example orbits.

To demonstrate the FDM-R and FDM-RV approaches, three sets of initial guesses for the path and control history are selected to initialize the process. The initial guesses corresponding to the slack variables are always established such that 𝐠⁰(𝐗)=𝟎 (the superscript “0” indicates an initial guess). For each orbit, the initial guess of the control history is one that maximizes the out-of-plane force contributed by the sail along the initial trajectory. Derived analytically by McInnes [33, pages 115–118, 223-224], the sail-pitch angle that maximizes the out-of-plane thrust below the Moon is 𝛼^∗=−35.26^∘. The initial guess for the thrust vector is then𝐮̂^0𝑖=cosΩ𝑡_𝑖cos𝛼^∗𝐱̂+sinΩ𝑡_𝑖cos𝛼^∗𝐲̂−sin𝛼^∗𝐳̂.(43) The initial guess for each trajectory is a circle offset below the Moon. Three sample initial guesses appear in Figure 4. The dark-blue circle with a radius of 59,000 km is offset below the Moon by 23,000 km. The radius of the light-blue circle is 14,000 km and the 𝑧-offset is 54,000 km. The dimensions of the red circle are 67,000 km by 75,000 km. As previously described, the system is formulated in the Earth-Moon circular restricted three-body system, where the reference frame is fixed to the Earth and the Moon. Therefore, the Sun is initially along the negative 𝑥-axis and moves in this reference frame clockwise around the Earth and the Moon, per (6). In this frame, the spacecraft is initially on the opposite side of the Moon from the Sun and moves clockwise around the Moon. The arrows on the orbits represent the initial guess for the direction of the sail-face normal, 𝐮̂, at the specified point along the orbit. Note that their direction is generally away from the Sun. The arrows are spaced in time by a little more than one day.

(a) 3D view

(b) Side view

Figure 4 

Three example initial guesses.

When the estimates from Figure 4 initialize the FDM-R algorithm, the three orbits in Figure 5 result. While the dark-blue orbit is not far from its initial approximation, the red orbit and the light-blue orbit are obviously dissimilar from their respective initial guesses. This indicates that a good initial guess for the finite-difference method is not essential for converging on a trajectory that solves the equations of motion. When the same initial guesses are used as inputs for the FDM-RV process, the trajectories in Figure 6 result. For both the FDM-R and FDM-RV algorithms, the dark-blue orbits converged in less than 10 iterations. The other two orbits required between 15 and 20 iterations for both methods. When individual iterations are examined for these two orbits, it is observed that the general shape of the orbit emerges quickly from the initial guess and the large number of iterations results from small attitude and path changes at the extremities in the ±𝑦 directions to accommodate the elevation constraint. All solutions converge quadratically over the final four or five iterations.

(a) 3D view

(b) Side view

Figure 5 

Three converged orbits from the FDM-R process.

(a) 3D view

(b) Side view

Figure 6 

Three converged orbits from the FDM-RV process.

For all trials of the FDM-R and FDM-RV formulations (and in the collocation formulation described in the next section) it is observed that the corrections to the components of the initial guess of the control history are not smooth at the first and 𝑛th point; this phenomenon is apparent in the control history from the seventh-degree Gauss-Lobatto collocation formulation presented in Figure 7 of Ozimek et al. [2]. All corrections to the discretized path and attitude profiles are done simultaneously, and the first point and 𝑛th points are the only two points in this simulation that are constrained to be equal to each other to ensure periodicity. Because the remaining points in the attitude profile do appear to be continuous and smooth, a postfit modification is applied such that the first and 𝑛th points are shifted by interpolating with neighboring points in the periodic profile. This postfit modification is only applied when the Newton-Raphson solver is not within the radius of convergence (i.e., when the convergence criterion in (25) exceeds 0.1). Control profiles from subsequent iterations appear smooth thereafter.

(a) 3D view

(b) Side view

Figure 7 

Three converged orbits from the Hermite-Simpson collocation process.

Although the dark-blue and red orbits appear to be similar in both the FDM-R and FDM-RV results, they differ by up to 8000 km and 1800 km, respectively, at certain epochs along the respective paths. The difference in the light-blue orbits in Figures 5 and 6 strongly suggests that the different methods, FDM-R and FDM-RV, can produce different solutions to the equations of motion using the same initial guess. While the theoretical errors in the solutions from the respective algorithms are equivalent, the differences in the formulation of the FDM (either -R or -RV) drive the evolution of the trajectory from iteration to iteration. It is not appropriate to suggest that the solution from one method is more “correct” than the other; both methods theoretically solve the equations of motion to within their errors bounds, and both methods produce orbits that are not obvious when the design space is not well understood.

6. Comparison to a Collocation Method

As stated, prior to any analysis, the design space for this problem is unknown. Entry into the design space via analytic techniques is difficult because of coupling of the sail attitude and thrust vector. However, based on the FDM algorithms, solutions precise to 𝒪(Δ𝑡²), or 1740 km for this step size, in each direction emerge. In work by Ozimek et al. [2], a Hermite-Simpson collocation method for the LSP coverage problem that has a theoretical local error on the order of 𝒪(Δ𝑡⁵), or 0.5 km for the same step size as that from the FDM algorithms. The Hermite-Simpson collocation method is adapted here as a basis to evaluate the accuracy of the FDM algorithms.

Using the same initial guesses as those represented in Figure 4 to seed the collocation method, the three orbits in Figure 7 result. The three orbits from the Hermite-Simpson approach are clearly similar to those produced by the FDM-RV process and only differ by, at most, 3000 km for the red orbit, 450 km for the light-blue orbit, and 70 km for the dark-blue orbit. The dark-blue and red paths from the Hermite-Simpson process are also similar to the dark-blue and red paths from the FDM-R approach (differing by, at most, 8000 km and 3500 km, resp.), but the light-blue paths are significantly different between the two strategies. Since the theoretical local error in the Hermite-Simpson method is on the order of 𝒪(Δ𝑡⁵), or 0.5 km for this step size, these differences between the similar orbits are not surprising [30]. Given that the initial guesses for the three sample orbits do not resemble their associated converged trajectories, it is remarkable that the converged trajectories from all three numerical processes are similar.

There is one noteworthy difference between the Hermite-Simpson collocation method employed by Ozimek et al. and the adaptation used for this analysis. Ozimek et al. employ a small step along the imaginary axis to calculate the partial derivatives related to the defects. For the current analysis, a simple central difference approximation (CDA) represents the same partial derivatives. Both are numerical schemes for computing difficult partial derivatives and subject to roundoff errors, but the CDA is also subject to truncation error, while the imaginary-step method is not. Thus, the error for the CDA employed in the Hermite-Simpson collocation method is minimized at a step size greater than zero (on the order of 10⁻⁸), but the error for the imaginary-step method approaches the roundoff error as the step size approaches zero (on the order of machine precision) [46]. This difference in precision should not have large ramifications for the process to converge, but must be noted. Furthermore, because the partial derivatives are numerically formulated for the Hermite-Simpson collocation scheme in this analysis, the three orbits in Figure 7 took approximately 2 seconds to generate. The orbits in Figures 5 and 6 took approximately 1 second to generate.

The results from any procedure can be used to seed another, usually more precise, process. The solutions produced by the FDM-R and FDM-RV algorithms are used to initialize the Hermite-Simpson scheme. All orbits reconverge quadratically in three or four steps when the Hermite-Simpson strategy is incorporated. The differences in position and velocity, as well as the direction of the sail-face normal, along the 29.5 day trajectories appear in Figure 8. Each color in the figure corresponds to the respective “initial” orbits appearing in Figures 5 and 6. The differences are well within the theoretical errors for the finite-difference methods, and, as such, the augmented finite-difference methods produce accurate and appropriate reference orbits for preliminary mission analysis.

(a) FDM-R solutions as initial guesses

(b) FDM-RV solutions as initial guesses

Figure 8 

Differences between the solutions from the Hermite-Simpson method and their initial guesses supplied by the FDM-R (a) and FDM-RV (b) solutions.

7. Conclusions

Augmented finite-difference methods are useful tools for generating trajectories in poorly understood dynamical regimes, such as the mechanics of flying a solar sail in a multibody system. These schemes are simple to understand and implement, notably in the presence of path constraints. While the theoretical errors of finite-difference methods may be larger than other similar methods that examine the trajectory as a whole, such as collocation, finite-difference methods can provide a reasonable entrance to the design space of a complicated nonlinear problem. They also quickly return results in a regime where little intuition concerning the trajectory and sail-angle history exist. Because of its speed and simplicity, this method may also serve as the basis for generating a large survey of orbit options. Once a viable trajectory is uncovered by the finite-difference method, other techniques may be employed to refine it further. The FDM approach developed for this analysis need not be applied to just solar sail problems, but should be applicable to a large variety of mission design problems.

Acknowledgments

The authors would like to thank Purdue Professor William Crossley for his explanation of optimization techniques in the literature used to solve similar problems and Professor Ahmed Sameh, also from Purdue University, for his insights on numerical methods used to solve boundary-value problems. Former Purdue graduate students Daniel Grebow and Zubin Olikara were invaluable to the initial development and understanding of the numerical methods examined in this paper. The authors also are grateful for the comments from the Editor and reviewers.