Research Article  Open Access
Geoffrey G. Wawrzyniak, Kathleen C. Howell, "Generating Solar Sail Trajectories in the EarthMoon System Using Augmented FiniteDifference Methods", International Journal of Aerospace Engineering, vol. 2011, Article ID 476197, 13 pages, 2011. https://doi.org/10.1155/2011/476197
Generating Solar Sail Trajectories in the EarthMoon System Using Augmented FiniteDifference 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 singlespacecraft relay to support communications with an outpost at the lunar south pole. Although trajectory design within the context of the EarthMoon 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 finitedifference schemes are simple to understand and implement. Two augmented finitedifference methods (FDMs) are developed and compared to a HermiteSimpson 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 higherfidelity models. The primary purpose of using a simple, lowerfidelity, augmented finitedifference 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 spacebased 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 EarthMoon link [1].
In contrast to a constellation of multiple spacecraft, alternative communications strategies that rely on only one satellite do exist, using current or nearfuture technology. Advanced propulsion concepts, such as lowthrust 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 EarthMoon 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. Differentialcorrections schemes [10] and shooting methods (both adaptations of the initial value problem), collocation approaches, and finitedifference 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 finitedifference 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 highprecision scheme; solutions from those methods may not persevere without additional control in a higherfidelity model (e.g., a model that includes planetary ephemerides) or in actual flight. If a solution based on a highprecision approach (e.g., shooting or highdegree collocation) is migrated from a model with limited accuracy to one of higher accuracy (e.g., circular restricted threebody model versus ephemeris model or actual flight), it may be no more accurate in that highfidelity model than a solution from a lowprecision approach (e.g., finitedifferencing). 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, higherfidelity 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 EarthMoon circular restricted threebody (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 twopoint 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 EarthMoon 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 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 angularvelocity 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: where represents the mass fraction of the smaller body, or , and and are the distances from the larger and smaller bodies, respectively, that is, Solar gravity is neglected in this model. At a distance of 1 AU, an appropriate distance to assume for a sailcraft in the EarthMoon system, the applied acceleration from a solar sail is modeled as where is the sailface normal, is a unit vector in the suntospacecraft 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 , where is the sail pitch angle, or the angle between the solar incidence direction and the sail normal.
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 EarthMoon distance), to the square of the characteristic time, ( days), that is, A recent sailcraft design for NASA's Space Technology competition (ST9) that was built by L'Garde possesses overall characteristic acceleration, , of , while the characteristic acceleration of the sail and its support structure alone is closer to (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, 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 where is .
The sail modeled here is a perfectly reflecting, flat solar sail. Billowing is not incorporated in this force model; however, could be less than , 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 finiteelement 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 sailface normal but, instead, is increasingly offset from the sailface 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 and , 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 finitedifference methods for generating solar sail trajectories.
3. Augmented FiniteDifference Methods
The derivation for an augmented finitedifference method begins with generic, nonlinear, secondorder, twopoint BVP that can be represented as 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. Centraldifference 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 where is defined as and and the symbol โโ indicates a numerical approximation. Equation (10) arises from three centraldifferences: one for the velocity at , which is midway between and , another at , midway between and , 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 , , and are fixed at a value of , then the midpoint of [ is and (9) and (10) simplify to Note that is used for when calculating for a periodic solution. Similarly, substitutes for in , so the problem does not include redundant constraints. The calculations for and are similar. The differences between the actual velocity and acceleration and their respective CDAs at are where and and are continuous on . 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 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 (FDMR). The second FDM is formulated to directly solve for position and velocity along the discretized trajectory (FDMRV).
The FDMR formulation is based on approximating (1) via the CDAs for velocity and acceleration given by (11) and (12), respectively, that is, 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: where โโ 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 FDMR are collected into the subvector, where the length of is . The subvectors corresponding to each node are collected into a column vector, which represents the discretized trajectory as well as the associated control history and slack variables. The velocities associated with the trajectories generated by the FDMR can be reconstructed with CDAs of during postprocessing.
The FDMRV process is similarly formulated, but there are two distinct sets of ODEs to solve, that is, Note the difference in the velocity terms in (16) and (21). With this FDMRV alternative, the subvector is since is now explicitly part of the solution. This formulation results in an vector with length .
Both the FDMR and the FDMRV 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, controldirection 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 , that is, where the is the Jacobian, . Superscripts on vectors refer to iteration number. The initial guess is denoted ; 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 (i.e., ), (23) is rearranged into a leastnorm problem, that is, Equations (23) and (24) reflect the NewtonRaphson 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, When an initial guess is in the neighborhood of a solution, convergence is quadratic. Specifying is sufficient, since any further iteration approaches double precision. If the initial guess is not in the neighborhood of a solution, bears little resemblance to ; the step from to may even be chaotic. However, the new 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 leastnorm update in (24) minimizes the norm of the correction to the state vector; the state vectors for the two formulations are different, and the leastnorm 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 FDMR approach, (8) is written as The pseudogravity gradient associated with the true solution, , is expanded as a function of an approximate solution, , that is, where and is the Hessian of . Using the solution from the FDMR algorithm, (15) is subtracted from (26), resulting in where . It is assumed that is continuously differentiable to fourth order. Rearranging and multiplying by yields where is the identity matrix and is the skewsymmetric gyroscopic matrix. Some terms in (29) can be simplified to a form that will be useful later, that is Now, let correspond to the difference element that possesses the largest magnitude. Equation (29) reduces to an upper bound on , Using (30), (31) simplifies to The inequality in (34) represents three scalar equations. Thus, the position error corresponding to the FDMR formulation is as . To calculate the velocity from the FDMR algorithm, the CDA from (11) is used, also with an error . The error analysis for the alternate FDMRV 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 FDMR formulation beginning with (26).
The augmented finitedifference 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, , or , which translates to an error of of the EarthMoon 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 finitedifference 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 , The presence of depends on whether the formulation is the FDMR or FDMRV, and the length of the vector is either or , 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, bandeddiagonal 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: 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 FDMRV algorithm) in : 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 NewtonRaphson 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, 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 FDMR or an FDMRV 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, . 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 Formulating the constraint in terms of as opposed to 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, , maintains the visibility of the spacecraft from an outpost located near the south pole of the Moon, and the spacecraft altitude constraint, , is imposed for radio power restrictions. Altitude is defined as the distance from the lunar south pole, that is, The third path constraint requires that the sailface normal, or control , is always directed away from the Sun (the sunlight vector is ), or in (7). Of the inequality constraints, only this attitude requirement is mandated. Together, these three inequality constraints are written as 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: Combining the path constraints at all epochs results in total elements in . When using the path constraints in (42), . The first two path constraints in (42) as well as the periodicity constraint in (38) are illustrated in Figure 2. The 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.
Only the first 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 FDMRV 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 FDMR formulation is similar to that of the FDMRV, except that the diagonal corresponding to does not exist in the FDMR option. The size of the Jacobian for the FDMR option associated with 101 nodes is 710 rows by 909 columns; approximately 0.62% of the entries are nonzero.
The advantage of a finitedifference 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]. Thirdparty 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 . 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 lowerlevel 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.

In summary, a finitedifference method yields a solution for the path by replacing the path derivatives with their finitedifference 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 finitedifference 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 is selected. A generous maximumaltitude constraint of one EarthMoon 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 mm/s^{2} [34]. JAXA's IKAROS is the first mission to have flown using a sail as its only propulsive device, and NASA recently deployed the NanoSailD2. Both of these sailcraft are relatively small (200โm^{2} and 10โm^{2}, resp.) and are designed to demonstrate inspace 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 nearfuture technological improvements. An assessment of the technology level (i.e., characteristic acceleration) is required for a sailonly LSP mission. The focus of this analysis is not to demonstrate that a sailonly LSP mission is feasible based on current technology. Instead, the results of this investigation demonstrate the use of finitedifference methods for the sample application of an LSP mission. Therefore, a high characteristic acceleration, that is, mm/s^{2}, is used to generate example orbits.
To demonstrate the FDMR and FDMRV 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 outofplane force contributed by the sail along the initial trajectory. Derived analytically by McInnes [33, pages 115โ118, 223224], the sailpitch angle that maximizes the outofplane thrust below the Moon is . The initial guess for the thrust vector is then The initial guess for each trajectory is a circle offset below the Moon. Three sample initial guesses appear in Figure 4. The darkblue circle with a radius of 59,000โkm is offset below the Moon by 23,000โkm. The radius of the lightblue 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 EarthMoon circular restricted threebody 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 sailface 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
When the estimates from Figure 4 initialize the FDMR algorithm, the three orbits in Figure 5 result. While the darkblue orbit is not far from its initial approximation, the red orbit and the lightblue orbit are obviously dissimilar from their respective initial guesses. This indicates that a good initial guess for the finitedifference 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 FDMRV process, the trajectories in Figure 6 result. For both the FDMR and FDMRV algorithms, the darkblue 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
(a) 3D view
(b) Side view
For all trials of the FDMR and FDMRV 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 seventhdegree GaussLobatto 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 NewtonRaphson 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
Although the darkblue and red orbits appear to be similar in both the FDMR and FDMRV results, they differ by up to 8000โkm and 1800โkm, respectively, at certain epochs along the respective paths. The difference in the lightblue orbits in Figures 5 and 6 strongly suggests that the different methods, FDMR and FDMRV, 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 HermiteSimpson 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 HermiteSimpson 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 HermiteSimpson approach are clearly similar to those produced by the FDMRV process and only differ by, at most, 3000โkm for the red orbit, 450โkm for the lightblue orbit, and 70โkm for the darkblue orbit. The darkblue and red paths from the HermiteSimpson process are also similar to the darkblue and red paths from the FDMR approach (differing by, at most, 8000โkm and 3500โkm, resp.), but the lightblue paths are significantly different between the two strategies. Since the theoretical local error in the HermiteSimpson 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 HermiteSimpson 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 imaginarystep method is not. Thus, the error for the CDA employed in the HermiteSimpson collocation method is minimized at a step size greater than zero (on the order of ), but the error for the imaginarystep 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 HermiteSimpson 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 FDMR and FDMRV algorithms are used to initialize the HermiteSimpson scheme. All orbits reconverge quadratically in three or four steps when the HermiteSimpson strategy is incorporated. The differences in position and velocity, as well as the direction of the sailface 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 finitedifference methods, and, as such, the augmented finitedifference methods produce accurate and appropriate reference orbits for preliminary mission analysis.
(a) FDMR solutions as initial guesses
(b) FDMRV solutions as initial guesses
7. Conclusions
Augmented finitedifference 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 finitedifference methods may be larger than other similar methods that examine the trajectory as a whole, such as collocation, finitedifference 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 sailangle 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 finitedifference 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 boundaryvalue 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.
References
 A. Ely and E. Lieb, โConstellations of elliptical inclined lunar orbits providing polar and global coverage,โ in Advances in the Astronautical Sciences, vol. 123, pp. 1447โ1462, Lake Tahoe, Calif, USA, 2005, paper no. AAS 05343. View at: Google Scholar
 M. T. Ozimek, D. J. Grebow, and K. C. Howell, โDesign of solar sail trajectories with applications to lunar south pole coverage,โ Journal of Guidance, Control, and Dynamics, vol. 32, no. 6, pp. 1884โ1897, 2009. View at: Publisher Site  Google Scholar
 G. G. Wawrzyniak and K. C. Howell, โThe solar sail lunar relay system: an application of solar sails in the Earth–Moon system,โ in Proceedings of the 59th International Astronautical Congress, Glasgow, Scotland, September 2008, paper no. IAC08.C1.3.14. View at: Google Scholar
 D. J. Grebow, M. T. Ozimek, and K. C. Howell, โDesign of optimal lowthrust lunar polesitter missions,โ in Advances in the Astronautical Sciences, vol. 134, pp. 741โ760, Savannah, Ga, USA, 2009, paper no. AAS 09148. View at: Google Scholar
 R. W. Farquhar, โThe control and use of libration point satellites,โ Tech. Rep. TRR346, National Aeronautics and Space Administration, February 1970. View at: Google Scholar
 C. R. McInnes, โSolar sail trajectories at the lunar ${L}_{2}$ Lagrange point,โ Journal of Spacecraft and Rockets, vol. 30, no. 6, pp. 782โ784, 1993. View at: Google Scholar
 C. R. McInnes, A. J. McDonald, J. F. L. Simmons, and E. W. MacDonald, โSolar sail parking in restricted threebody systems,โ Journal of Guidance, Control, and Dynamics, vol. 17, no. 2, pp. 399โ406, 1994. View at: Google Scholar
 J. Simo and C. R. McInnes, โSolar sail trajectories at the Earth–Moon Lagrange points,โ in Proceedings of the 59th International Astronautical Congress, Glasgow, Scotland, 2008, paper no. IAC08.C1.3.13. View at: Google Scholar
 J. Simo and C. R. McInnes, โSolar sail orbits at the EarthMoon libration points,โ Communications in Nonlinear Science and Numerical Simulation, vol. 14, no. 12, pp. 4191โ4196, 2009. View at: Publisher Site  Google Scholar
 K. C. Howell and H. J. Pernicka, โNumerical determination of Lissajous trajectories in the restricted threebody problem,โ Celestial Mechanics, vol. 41, no. 1–4, pp. 107โ124, 1987. View at: Publisher Site  Google Scholar
 H. B. Keller, Numerical Methods for TwoPoint BoundaryValue Problems, Blaisdell, Waltham, Mass, USA, 1968.
 S. D. Conte and C. de Boor, Elementary Numerical Analysis, McGrawHill, New York, NY, USA, 1980.
 D. R. Kincaid and E. W. Cheney, Numerical Analysis: Mathematics of Scientific Computing, American Mathematical Society, Providence, RI, USA, 2002.
 J. Stoer and R. Bulirsch, Introduction to Numerical Analysis, Springer, New York, NY, USA, 3rd edition, 2002.
 E. W. Brown, An Introductory Treatise on Lunar Theory, Cambridge University Press, London, UK, 1896.
 D. L. Richardson, โAnalytic construction of periodic orbits about the collinear points,โ Celestial Mechanics, vol. 22, no. 3, pp. 241โ253, 1980. View at: Publisher Site  Google Scholar
 S. H. Strogatz, Nonlinear Dynamics and Chaos: with Applications to Physics, Biology, Chemistry, and Engineering, Westview Press, Cambridge, Mass, USA, 2001.
 D. Izzo, โGlobal optimization and space pruning for spacecraft trajectory design,โ in Spacecraft Trajectory Optimization, B. A. Conway, Ed., Cambridge University Press, New York, NY, USA, 2010. View at: Google Scholar
 M. Pontani and B. A. Conway, โSwarming theory applied to space trajectory optimization,โ in Spacecraft Trajectory Optimization, B. A. Conway, Ed., Cambridge University Press, New York, NY, USA, 2010. View at: Google Scholar
 R. S. Wilson and K. C. Howell, โTrajectory design in the SunEarthMoon system using lunar gravity assists,โ Journal of Spacecraft and Rockets, vol. 35, no. 2, pp. 191โ198, 1998. View at: Google Scholar
 B. G. Marchand, K. C. Howell, and R. S. Wilson, โImproved corrections process for constrained trajectory design in the nbody problem,โ Journal of Spacecraft and Rockets, vol. 44, no. 4, pp. 884โ897, 2007. View at: Publisher Site  Google Scholar
 O. von Stryk and R. Bulirsch, โDirect and indirect methods for trajectory optimization,โ Annals of Operations Research, vol. 37, no. 1, pp. 357โ373, 1992. View at: Publisher Site  Google Scholar  MathSciNet
 J. T. Betts, โSurvey of numerical methods for trajectory optimization,โ Journal of Guidance, Control, and Dynamics, vol. 22, no. 2, pp. 193โ207, 1998. View at: Google Scholar
 R. D. Russell and L. F. Shampine, โA collocation method for boundary value problems,โ Numerische Mathematik, vol. 19, no. 1, pp. 1โ28, 1972. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 E. J. Doedel, โFinite difference collocation methods for nonlinear two point boundary value problems,โ SIAM Journal on Numerical Analysis, vol. 16, pp. 173โ185, 1979. View at: Google Scholar
 C. R. Hargraves and S. W. Paris, โDirect trajectory optimization using nonlinear programming and collocation,โ Journal of Guidance, Control, and Dynamics, vol. 10, no. 4, pp. 338โ342, 1987. View at: Google Scholar
 P. J. Enright and B. A. Conway, โOptimal finitethrust spacecraft trajectories using collocation and nonlinear programming,โ Journal of Guidance, Control, and Dynamics, vol. 14, no. 5, pp. 981โ985, 1991. View at: Google Scholar
 P. J. Enright and B. A. Conway, โDiscrete approximations to optimal trajectories using direct transcription and nonlinear programming,โ Journal of Guidance, Control, and Dynamics, vol. 15, no. 4, pp. 994โ1002, 1992. View at: Google Scholar
 S. Tang and B. A. Conway, โOptimization of lowthrust interplanetary trajectories using collocation and nonlinear programming,โ Journal of Guidance, Control, and Dynamics, vol. 18, pp. 599โ604, 1995. View at: Google Scholar
 A. L. Herman and B. A. Conway, โDirect optimization using collocation based on highorder GaussLobatto quadrature rules,โ Journal of Guidance, Control, and Dynamics, vol. 19, no. 3, pp. 592โ599, 1996. View at: Google Scholar
 R. G. Melton, โComparison of direct optimization methods applied to solar sail problems,โ in Proceedings of the AIAA/AAS Astrodynamics Specialists Conference and Exhibit, Monterey, Calif, USA, August 2002, paper no. AIAA 20024728. View at: Google Scholar
 G. G. Wawrzyniak and K. C. Howell, โInvestigating the design space for solar sail trajectories in the Earth—Moon system using augmented finitedifference methods,โ Purdue University Document, School of Aeronautics and Astronautics, February 2011. View at: Google Scholar
 C. R. McInnes, Solar Sailing: Technology, Dynamics and Mission Applications. Space Science and Technology, Springer, New York, NY, USA, 1999.
 D. Lichodziejewski, B. Derbes, D. Sleight, and T. Mann, โVacuum deployment and testing of a 20m solar sail system,โ in Proceedings of the 47th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Newport, RI, USA, May 2006, paper no. AIAA20061705. View at: Google Scholar
 B. Derbes, D. Lichodziejewski, and G. Veal, โA “yank and yaw” control system for solar sails,โ in Advances in the Astronautical Sciences, vol. 19, pp. 2893โ2907, 2004. View at: Google Scholar
 L. RiosReyes and D. J. Scheeres, โSolarsail navigation: estimation of force, moments, and optical parameters,โ Journal of Guidance, Control, and Dynamics, vol. 30, no. 3, pp. 660โ668, 2007. View at: Publisher Site  Google Scholar
 B. Campbell, An analysis of thrust of a realistic solar sail with focus on a flight validation missionin a geocentric orbit, Ph.D. Thesis, The George Washington University, Washington, DC, USA, 2010.
 The MathWorksTM, Inc, โMATLAB R, Version 7.11.0 (R2010b),โ August 2010, http://www.mathworks.com/. View at: Google Scholar
 A. L. Herman, Improved collocation methods with application to direct trajectory optimization, Ph.D. Thesis, University of Illinois at UrbanaChampaign, 1995.
 J. T. Betts, Practical Methods for Optimal Control Using Nonlinear Programming, SIAM, Philadelphia, Pa, USA, 2001.
 C. R. Ortiz Longo and S. L. Rickman, โMethod for the calculation of spacecraft umbra and penumbrashadow terminator points,โ Tech. Rep. 3547, Lyndon B. Johnson Space Center, Houston,Tex, USA, 1995. View at: Google Scholar
 D. E. Salane, โAdaptive routines for forming Jacobians numerically,โ Tech. Rep. SAND861319, Sandia National Laboratories, Albuquerque, NM, USA, 1986. View at: Google Scholar
 S. A. Forth and M. Edvall, User Guide for MAD—A Matlab Automatic Differentiation Package, TOMLAB/MAD, Version 1.4 The Forward Mode, Engineering Systems Department, Defence College of Management & Technology, Cranfield University, Shrivenham, UK, 2007.
 L. F. Shampine, โAccurate numerical derivatives in MATLAB,โ ACM Transactions on Mathematical Software, vol. 33, no. 4, Article ID 1268781, 2007. View at: Publisher Site  Google Scholar
 M. Macdonald and C. McInnes, โSolar sail mission applications and future advancement,โ in Proceedings of the 2nd International Symposium on Solar Sailing, pp. 1โ26, New York City College of Technology, City University of New York, New York, NY, USA, July 2010. View at: Google Scholar
 J. Martins, I. Kroo, and J. Alonso, โAn automated method for sensitivity analysis using complex variables,โ in Proceedings of the 38th Aerospace Sciences Meeting and Exhibit, Reno, Nev, USA, January 2000, paper no. AIAA 20000689. View at: Google Scholar
Copyright
Copyright © 2011 Geoffrey G. Wawrzyniak and Kathleen C. Howell. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.