OPTIMAL TRANSFER AND SCIENCE ORBIT DESIGN IN THE PROXIMITY OF SMALL BODIES Christopher J. Scott1 and Martin T. Ozimek2 This study details a robust and systematic approach for generating science orbits connected by fuel optimal transfers in the small body regime. The Clohessey-Wiltshire equations provide an analytical means to globally interpret the design space allowing a designer to quickly investigate families of science orbits and obtain closed-form optimal transfer. Solutions are refined with multiple shooting and nonlinear programming to transition the results into Hill’s equations, and finally, a high-fidelity ephemeris model. The final result is an accurate solution retaining the initial design characteristics with the potential for further fuel savings that are available from exploiting the high-fidelity perturbations.
INTRODUCTION This study demonstrates the calculation of optimal transfers and rendezvous about small bodies in orbit about a larger body such as the Sun or a planet. To date there have been many proposed and a few flown missions to small bodies including the Near-Earth Asteroid Rendezvous (NEAR) mission to Eros and Japan’s Hayabusa mission to Itokawa. The failed Phobos 2 spacecraft of the former Soviet Union was intended to pass within 50 m of Phobos and release a mobile lander. Numerous other proposed missions have intended to visit Martian system including the Mars-Moon Exploration, Reconnaissance and Landed Investigation (MERLIN), Phobos-Grunt, Phobos Reconnaissance and International Mars Exploration (PRIME), Phobos-Deimos Sample Return Mission (SRM), and Aladdin. This study was motivated by mission designs of these missions. Although the ∆V allocations for proximity operations are often small compared to the total ∆V budget, optimal use of the remaining fuel can significantly increase the mission lifetime and enable proper spacecraft disposal. There can be an order of magnitude difference in terms of ∆V among optimal and suboptimal transfers. A thorough understanding of the design space is required for the timely design of contingency maneuvers and the rapid redesign of nominal science orbits due to unforeseen targets of opportunity and more accurate estimates of the physical parameters of the small body or spacecraft state. Invariably the cadence of these activities will be accelerated compared to missions about well understood bodies. Both NEAR and Hayabusa required design updates due to failed maneuvers and other spacecraft anomalies. Moreover, the mathematics described herein relate to applications involving larger maneuvers such as transfers among distant retrograde orbits in the Earth-sun system. A transformed solution of the Clohessy-Wiltshire (CW) equations with phasing and energy constraints derived from Hill’s problem is used as an initial approximation in the search for optimal analytical solutions. It is advantageous to initially place a spacecraft on a relatively high-altitude orbit for initial shape and orbit determination, and regional characterization of the body. For example, OSIRIS-REx will perform reconnaissance of 1999 RQ36 for about year before landing. The characteristics of small bodies are generally not known to enough accuracy for precise science orbit planning. Because of the timescales in the Martian system such orbits are often called “pseudo-orbits” or “formation orbits” and are related to the distant retrograde orbits of the circular restricted three-body problem (CR3BP) and Hill’s problem. In 1
Mission Design Analyst, The Johns Hopkins University Applied Physics Laboratory, 11100 Johns Hopkins Road, Laurel MD, 20723. 2 Mission Design Analyst, The Johns Hopkins University Applied Physics Laboratory, 11100 Johns Hopkins Road, Laurel MD, 20723.
the orbital plane of the small body with respect to a larger central body these orbits appear as 2×1 relative ellipses. A number of key results are obtained from the analytical analysis including initial estimates of optimal transfer conditions. As the characteristic dimensions of these relative orbits decreases or the propagation times are increased the accuracy of the predictions degrades. To account for dynamically sensitive regions and the inherent inaccuracies of the approximation, the solutions are transitioned into Hill’s equations which are an approximation of the CR3BP valid for small mass ratios. With Hill’s equations many of the near symmetries of the CR3BP are transformed into exact symmetries simplifying numerical algorithms and analyses. Planar periodic retrograde orbits with an approximate 2:1 axis ratio exist for a large range of distances from the smaller primary. Additionally, a wide range of 3-D periodic retrograde orbits exist if certain phasing constraints are met. These periodic orbits are produced using minimum-norm correction and pseudo-arclength continuation1. The optimal solutions found using the CW equations are refined in Hill’s equations with SNOPT which utilizes sequential quadratic programming2. The algorithm proves to be robust even under conditions where the CW equations are inaccurate. The periodic orbits of Hill’s equations are transitioned into a full-ephemeris model of the solar system with multiple shooting. The transfers among these natural orbits are again refined using SNOPT for variable time transfers. MODELS Clohessy–Wiltshire (CW) equations The CW equations assume that the target satellite is “massless” and revolves in a circular orbit about the central body. 2 3 0 2 0
0
(1)
These equations comprise a set of linear coupled differential equations. Normalizing so that the characteristic length and the rate of rotation equal one yields the solution, Φ, . In order to express this solution in a geometrically descriptive form, Lovell et al.3 introduced the following formation parameters, 2 3 2 3 2
4 2 2
! " #
The equations of motion can then be written as $%& &' 2 3 &' $%& 2 2 where
3 2
# #
# &'
(2)
(3)
# $%&
Here β is the parametric angle and the subscript “o” indicates an initial condition. Because the xdirection points away from the center of the gravitating body, the spacecraft is at perigee when β equals zero. Likewise, it is at apogee when β equals (. Examining the x- and y-equations reveals that the projection of the relative trajectory into the y–x plane is a 2 ×1 ellipse, where a is the length of the minor
axis. Because motion in the z-direction is decoupled from motion in the x- and y-directions, the spacecraft follows the motion of a simple harmonic oscillator along this axis. If , equivalently , is equal to zero, then the formation has no secular changes. When equals 0 or ( the spacecraft crosses the orbital plane of the target. Hill’s Equations
Hill’s equations are an approximation of the CR3BP valid for small mass ratios. Using the rotating CW coordinate frame centered at the smaller primary and letting ) be the gravitational parameter and , the distance from the smaller primary, the following equations can be derived by setting the length unit to )⁄ and the time unit to 1⁄. 2 3 , (4) 2 ,
, This system possesses a constant of motion 1 1 1 . 3 2 , 2 and the following symmetries
(5)
/
, , , , , , 0 , , , , , , /
, , , , , , 0 , , , , , , /-
, , , , , , 0 , , , , , , /1
(6)
, , , , , , 0 , , , , , , /2
, , , , , , 0 , , , , , ,
Note that these symmetries are also present in the CW equations. Symmetries s1 and s2 are reflections in time while s3 is a reflection about the x-y plane. The s4 and s5 symmetries are a combination of s1 and s2, and, s1, s2 and s3, respectively. The last represents symmetry about the origin. Ephemeris Model A general n-body model is used for high-fidelity propagations given by the equations 34 5 6 7 8
9>4
34 39 :5; 56 < -
=5; 56 =
(7)
The position of each body is specified by JPL’s Navigation and Ancillary Information Facility (NAIF) standard DE421 ephemerides. PERIODIC ORBITS Hill’s Equations Under certain phasing constraints a subset of the formation orbits of the CW equations can be continued into the Hill’s model. It is descriptive to consider sets of constant J. For the CW equations to
remain a valid approximation the orbit cannot pass into regions that are gravitationally dominated by the small body. Thus, assuming that 1⁄, ? 1 and no drift, 1 .@ A 4 # 8
(8)
⁄
where the constant has been scaled by the mass ratio of Deimos and Mars, .@ .)CD . If this type of planar periodic orbit is iso-energetically continued into orbits traveling out-of-plane the energy can be expressed as .@ ⁄8 where is the initial y-amplitude. Thus with increasing # the family will evolve into an orbit that has diminishing amplitude in the x-y plane. We consider periodic orbits with two perpendicular crossings of the x-axis known as axi-ssymmetric periodic orbits. This sets the relative phasing between in-plane and out-of-plane motion to be ( or depending on the sign convention used for the z-amplitude.
It is useful to consider the altitude over a particular centric latitude and longitude on an energy surface. Considering a gravity locked body and locating this point by the unit vector E FGH GI GJ KL gives a centric latitude, M, and longitude, N, of M ,$&'GJ
N ( ,$:GI ⁄GH <
(9)
with the appropriate quadrant check. Here zero longitude corresponds to the hemisphere which faces the larger primary. If the dimension of the body facing the larger primary is D, then #OP 2Q and
#,#RH S #OP . This gives a theoretical upper limit for the centric latitude at a set energy,
Q 1 4 ! V " Z (10) M#RH ,$&' UW Y Q 12 ! " 1 T X Although in practice numerical sensitivity usually prohibits calculation of the higher latitude orbits, this equation shows that perfectly polar orbits are impossible without collision by energy considerations. Considering the orbit which passes directly above the point yields the amplitude ratio,
# GJ ! " [ \ GI
yeilding,
.@ A
with the appropriate quadrant checks,
(11)
GJ ]1 4 [ \ ^ 8 GI
GJ 1 GJ GI _ ` ` 2 # GH 2 GH GJ
(12)
(13)
where the subscript “f” indicates the flyover point on the orbit. The magnitude of the radius at the flyover can be reduced to a function of E and . ,_ aE
(14)
Similar relationships can be used to predict the projection of the family on the x-z, and y-z planes. The maximum distance along any direction is 3E · 5 where,
therefore,
5 $%&c &'6 sin ; 2 2 E · 5 g , , E$%&: Φ , , E<
and 3E · 5 g , , E . Projecting into the x-y, y-z, and x-z plane gives, 3E · 5
S4 3GI
,
1 3GJ
,
1 h 4 8GJ GJ1
(15)
(16)
(17)
This equation traces out a square in and is shown in Fig. 1. In the x-y and y-z planes the boundaries are ellipses with a dimension of in the y-direction and ⁄2 in the x- and z- directions. The x-z boundary is related to the area enclosed by well-studied Lissajous curves, with similar geometries appearing near the Lagrange points. Using pseudo-arclength continuation the entirety of formation orbits at a constant energy can be generated. With this method steps are taken tangent to the family, rather than with physical parameters such as initial position. This typically facilitates much larger step sizes, which allows for quick generation of entire families of periodic orbits. Hill’s problem possesses symmetry about the origin; therefore, a converged solution of an axi-symmetric periodic orbit, iO , will satisfy a set of constraints a quarter of a period into the orbit, jiO k l 0 . .
(18)
where . sets the energy. There is one more free variable, Fm nKo , than constraints, so that the direction tangent to the family is the null vector, ∆iO , of the Jacobian, QjiO . Setting a step size of, ∆&, the constraint vector can be appended,
qc r u 0 . . cst iO ∆iOs ∆&
(19)
Thus iterating, Os c Qq⁄q, until convergence tolerances are met yields the next solution, iO .
Figure 1 shows a set of formation orbits scaled to the Mars-Deimos system. Here the initial orbit is continued until the family returns to the original orbit. These orbits agree closely to the analytical results outlined above.
Figure 1. Formation orbits about Deimos using Hill’s equations. Analytically predicted boundaries are shown on the edges.
Using Eq. 9 the centric groundtracks are shown in Fig. 2.
Figure 2. Centric groundtracks for formation orbits shown in Fig. 1.
Full Ephemeris Multiple Shooting A fixed-time multiple shooting method was used to convert families of periodic orbits in Hill’s equations to full-ephemeris propagations. This method involved setting a number of patch points and integrating corresponding trajectory arcs throughout the entire segment. The procedure attenuates the inherent instability of a single shooting method over lengthy timespans. Combining these points into a single 6n vector v v (20) r w u 0 vP
In order to remove the discontinuities at the patch points the following 6(n-1) constraints must be satisfied, v{ v z { ~ j y v- v- } 0 w y } xv{P vP |
(21)
The resulting 6(n-1)×6n Jacobian is incorporated into a minimum norm correction routine to update the patch points, Os O jO o FjO jO o K jO . The integrations are performed using DDEABM, a variable order (1-12), Adams-Bashforth-Mouton predictor corrector method that has been shown to perform well relative to other methods for the desired tolerances.4 Applying this routine to orbits shown in Fig. 1 over four revolutions yields the orbits shown in Fig. 3. The n-body propagations were performed with the gravity of the Sun, Mars, Deimos, and Jupiter. Qualitatively the orbits are similar except for noticeable distortion in the y-direction and numerical instability for the nearly vertical orbits.
Figure 3. Formation orbits about Deimos using the n-body propagator with Sun, Mars, Deimos, and Jupiter.
OPTIMAL TRANSFERS Analytical Prediction Using this parameterization the initial and final states for a transfer over time _ are expressed as, R &'
R
$%&
R
r$%& - u
5 &'
# &'
# $%& R
$%&: Δ _ < _
5_ r _ &' : Δ _ < _ u where
#_ &':
_ Δ
Δ _ <
and
_ # Δ
.
&' : Δ _ < z ~ _ y_ $%&: Δ _ < - _ } y } x #_ $%& : Δ _ < | Letting
5
,
,
(22)
and
Φ Φ Ω Ω the velocities immediately after the first maneuver and before the second Φ Φ Ω Ω are s m Ω Ω and Ω Ω . Therefore,
t s m m
and
t Δ S s s s
Δ S:_ _ < :_ _ < : _ _ <
(23)
(24)
After significant simplification, the magnitudes of the maneuvers take the following form ΔV a :, Δ, Δ, _ , , , _ , _ <
: , Δ, Δ, _ , , , _ , _ <$%& ! N : , Δ, Δ, _ <"
(25)
: , Δ, Δ, _ <$%& !2 N :, Δ, Δ, _ <" J : # , Δ , Δ , _ <$%& !
ΔV a :, Δ, Δ, _ , , , _ , _ <
N - : # , Δ , _ , Δ <"
: Δ, Δ, _ , , , _ , _ <$%& ! N : , Δ, Δ, _ <"
(26)
: , Δ, Δ, _ <$%& !2 N :, Δ, Δ, _ <" J : # , Δ , Δ , _ <$%&:
Using the harmonic addition theorem,
N- # , Δ , Δ <
ΔV ΔV F @ :, Δ, Δ, _ , , , _ , _ <
7 @ :, Δ, Δ, _ , , , _ , _ <$%& ! Φ :, Δ, Δ, _ <"
(27)
@ :, Δ, Δ, _ <$%& !2 Φ :, Δ, Δ, _ <" J@ : # , Δ , Δ , _ <$%& !
Φ- : # , Δ , Δ , _ <"
Substituting in the phasing and amplitude constraints and considering no drift yields, ΔV ΔV F:, Δ, Δ, O , _ <
:, Δ, Δ, _ , O , _ <$%& !2 Φ:, Δ, Δ, _ , O , _ <"
(28)
Thus, setting the initial and final energies, z-amplitude, and phase shift 3-D transfers can be reduced to an optimization this function reduces to a single variable, , of periodicity π.
ΔV ΔV ,where q signifies an arbitrary variable, corresponds to either a subset or all of the zeros of ΔV or equivalently V. To show this first consider, 2ΔV . Thus, 0 implies that ΔV ΔV 0. Expanding the expression ¡ ¡
Equation 28 has utility because the zeros of
¢ ¢ΔV ¢ΔV ΔV ΔV 2ΔV 2ΔV ¢£ ¢£ ¢£
and letting ΔV ΔV gives,
2ΔV !
¤
Therefore, ΔV ΔV 0 when
"
or
or 2ΔV !
¡
¤
"
(30)
0 and ΔV ΔV . Expanding the expression,
¢ ¢ΔV ¢ΔV 2ΔV ΔV 2ΔV 2ΔV ¢£ ¢£ ¢£
and letting ΔV ΔV gives,
(29)
2ΔV !
Similarly 2ΔV ΔV 0 when
¤
or
"
¡
or 2ΔV !
¤
"
0 and ΔV ΔV . Writing out
(31)
(32) , ¡
∂ΔV ¢ ¢ ΔV ΔV 2ΔV ΔV ∂q ¢£ ¢£
(33)
indicates that if the entire expression equals zero then each term independently goes to zero when ΔV ΔV or 2ΔV ΔV correspond to a subset of the zeros of ΔV ΔV . Thus the zeros of
¡
or
.
Suppose that we analytically solve for a zero £ which satisfies
If while continuing this solution along the variable £,
§
§
¨ 0 and
ΔV §
ΔV ¨ 0.
changes sign then the solution has bifurcated
to a maximum. Detection of this event indicates that two minimums exist in its neighborhood with ΔV © ΔV .
For a function taking the form ΔV ΔV an analytical expression can be found for the stationary values of . ∂ ΔV ΔV 2&'2 Φ 0 ∂
(34)
( Φ ª J %«« 2 where n is an integer. The minimum becomes ¬ΔV ΔV =
§®
¯
(35)
Generally Φ is a complicated expression, but takes on simplified forms when certain conditions are met. For no phasing change, ¬Φ|±² _ (36)
From the viewpoint of two-body mechanics about the larger primary, this corresponds to mean anomaly value of (⁄2 _ ⁄2 where n is odd. In the case of rendezvous near the origin, _ _ A 0 and Φ³ _
(37)
with no phasing constraints. For transfers among orbits with the same z-amplitude, including transfers within the x-y plane ¬ ∂ ΔV ΔV ` ∂Δ J
® ´§ ±J´ ,§ ,±²
∂VT ` ∂Δ J
¬
® ´§ ±J´ ,§ ,±²
0
(38)
Without phasing constraints, it may be necessary to find global minimums of V with Δ and as free
variables.
§
Thus, necessary conditions for a local minimum are " §
!
§
0,
=0, §
¨ 0, and
¨ 0. In practice, we find Δ 0 for these conditions to be met when #
#_ . The equations reveal that a finite phasing change is optimal if # © #_ . Maps such as the one for planar transfers, shown in Fig. 4a can be quickly produced using the above equations and continuation with _ .
(a)
(b)
Figure 4. (a) Location of minimums versus ¶· for all planar transfers. (b) Representative contours of ¸¹ versus ºm and º for planar transfers.
A representative plot of ΔV versus and Δ for a set of planar transfers for a set Δ and _ is shown in Fig. 4b. Inspection reveals that the proper choice of departure time and phasing change can result in an order of magnitude difference in terms of ΔV. Focusing on only planar transfers ΔV can be expressed as,
ΔV a $%&2 N
ΔV a $%&2 N
where Φ
»¤ s»
ª ΔV ΔV 2a $%&2 Φ
§
_
. By construction a ¨ 0, ¨ 0 , and a ¨ . It can be shown that the ratio ¬ ½ ¼
function of only _ . Combining the terms ∂ ¬ ΔV ¾ ¢ ®,±²
(39)
±²
8 N N a N N Â $%& 2 $%& 1 ÃÄÃ 2 2 À a $%& !N N " 2 ¬ 8 N N a N N Á $%& 2 $%& 1 %«« À N N 2 2 a $%& ! " ¿ 2
is a
(40)
Thus, bifurcation is depends only on _ . Letting
N N a N N 2 $%& 1 ÃÄÃ 2 2 ¬ ÅQn:_ < Æ N N a N N $%& 2 $%& 1 %«« 2 2 $%&
(41)
Plotting the first function over a large range of _ shows that maximums will remain as maximums for all practical values of transfer time.
Figure 5. Ço:¶· < versus ¶· for even n.
The roots of ÅQn:_ < for n odd are tabulated below. ¶·
0 1.93481 5.80764 6.63778 8.54956 9.1349 12.3487 12.7679 15.5367 18.7065 18.9877
Table 1: Roots of ÅQn:_ < for n odd.
Analytically Predicted Minimums 2 4 2 4 2 4 2 4 2 4 2
Length of interval 1.93481 3.87283 0.83014 1.91178 0.58534 3.2138 0.4192 2.7688 3.1698 0.2812
Putting the square of the total velocity in the following form, ΔV
a È2 $%&2 N $%&2 N
a a h :$%&2 N $%&2 N < $%&2 N $%&2 N É
shows that all possible zeros of ¬ ½ §
±²
(42)
are dependent only on _ . Thus, the bifurcation diagram,
shown in Fig. 4, applies to all transfers among relative ellipses for all and Δ.
FULL EPHEMERIS OPTIMIZATION The analytical predictions of the CW equations are used as an initial guess for transfers propagated using Hill’s equations. The solution is optimized using an SQP algorithm implementing the Broyden– Fletcher–Goldfarb–Shanno method (Found in MATLAB’s fmincon function). We seek to minimize ΔV subject to the state equations given in Eq. 4 and the nonlinear constraints
Where
j Ê
v0 Ë1 Ìc :% < Í Î v:a < Ë2 Ì· :% Δ a <
% Î Î Δ¸t , Ë1 Ë1 , Δ¸t Δ¸ Δ¸
(43)
(44)
and Î F0 0 0Ko . v0 and v:_ < are the initial and final states along the transfer, Ìc and Ì· are the states on the first and final ellipses at a normalized time, t, where 0 corresponds to 0. These constraints implicitly ensure the proper change in phasing is enforced. Figures 6 and show the CW optimal predicted ΔV and those found using Hill’s equations for no phasing change. The transfers start from a single orbit shown in Fig. 1, a =1810.44 km and end at each of the other orbits. With the exception of few points the data is generally indistinguishable at the scales for shorter transfer times, _ 0.16 0.79 days.
Figure 6. Optimal ¸¹ predicted by the CW equations, solid dots, and Hill’s equations, open dots, for º Î at the time shown.
Figure 7. Optimal ¸¹ predicted by the CW equations, solid dots, and Hill’s equations, open dots, for º Î at the time shown.
ΔV at the optimal location of and Δ 0 versus tf is plotted in Fig. 8. If tf is left as free variable then a simple gradient descent algorithm will generally make minor adjustments to the initial value. However, an algorithm which allows wider variation in the search space will generally result much larger variations in transfer times. By the CW equations ΔV generally oscillates with transfer time. The results below show that that these regular oscillations are broken in higher fidelity models yielding a single advantageous transfer time.
Figure 8. Optimal ¸¹ predicted by the CW equations versus ¶· for º Î.
If no phasing constraints are necessary then a family of transfers can be seeded from, # #_ since the optimal phasing change can be predicted analytically. A simple continuation algorithm works in this
case by setting the guesses for Δ and _ of the next step to be equal to converged values from the previous step.
Depending on the sign of
§
the corresponding value of can be predicted
analytically or calculated by a variety of simple search algorithms (i.e., tracking the bifurcations or calculating the roots of ). §
The solution obtained from Hill’s model is used as an initial guess for the full-ephemeris model for transfers among relative ellipses resulting from multiple shooting refinement. States from the initial and final orbits are sampled via a cubic spline. Note that a phasing constraint is not included in this optimization. ΔV is minimized subject to the state equations derived from Eq. 7 and . ΔV L 0 C Ó L 0 C Ó L 0 C Ó
SÔ :_ < L Ó Ô :_ < C Ó Ô :_ < C Ó L 0 z ~ 0 y L }
0 y L } y Ó } y Ó } x _ |
Ô Ó L _ j Ô Ó L _ Î
Ô Ó L _
(45)
(46)
(47)
Where the subscript “D” refers to a condition on the departure science orbit, the subscript “A” refers to a condition on the arrival science orbit, and the subscript “T” refers to a condition on the transfer trajectory. All numerical integrations of the transfer trajectory utilize DDEABM and the n-body equations of motion using ephemeris locations of the bodies via the SPICE toolkit. These numerical integrations are compiled from C into MATLAB MEX-functions to combine the speed of C-code with the useful functionality of existing MATLAB routines. The parameters Ó and Ó are time-like free variables that completely parameterize a departure and arrival state on the science orbits, respectively, without need for any numerical integration. (Here, Ó and Ó specifically chosen as Julian dates that are consistent with ephemeris positions of the gravitating bodies in the n-body problem.) The variable Ó thus serves as a reference epoch for the transfer trajectory so that it is consistent with the ephemeris positions at each location on the departure orbit. The total time of flight, tf, appears in parentheses because the optimization approach considers either time-fixed or time-free solutions. In the latter case, tf serves as an independent design variable, but in the former case tf is “frozen” to a constant user-selected value. Both cases may serve a mission benefit. For example, if ΔV-savings are of paramount concern without regard to time-of-flight, then the time-free optimal solution should be sought. If time-of-flight constraints are crucial to the overall design, then a user may fix the value, thus removing it as a free parameter in the optimization. Moreover, the time-fixed formulation is useful for parametric studies involving ΔV and flight time. The formulation follows a preliminary assumption of time-invariance of the arrival science orbit, which is not generally true of ephemeris trajectories. Instead, the arrival science orbit, like the departure orbit, depends on epoch in the n-body problem to locate the gravitating bodies. For this study, the proximity orbits selected typically exhibit only small variations from their Hill’s model counterparts, so an initial solution is still usually “close” to an actual insertion point on the arrival orbit. Nevertheless, the above approach is only a preliminary design that allows for robust optimization since each science orbit depends only on a single independent spline variable. Once the problem is solved once, the time invariance assumption may be iteratively removed by the following algorithm.
Essentially, this pseudo-algorithm algorithm gradually removes the discrepancy in the transfer transfe trajectory and ephemeris time.. To avoid “chatter” that may prevent the achievement of time continuity, the jmax flag allows the guaranteed nteed solution of the problem via a fixed fixed-time time feasible solution that meets all constraints and closely matches the time-invariant invariant optimal results. Typically jmax = 3 for this study. The typical cause for the chatter include the small variations in the cconverged quasi-periodic periodic arrival orbit states due to new locations of the ephemeris bodies upon each iteration. The algorithm represents a robust strategy for computing near-optimal optimal transfer trajectories while preserving the dynamical structure of the science scie orbits.
Figure 9. Ephemeris Transfer Optimization Diagram
Both phasing change and transfer time are left as free variables to optimize ΔV using both MATLAB’s implementation of its SQP algorithm (MSQP, fmincon function) and the Sparse Nonlinear OPTimizer (SNOPT). The results in Fig. 9 show that fmincon converges on a large span of values for either model depending on the initial guess. For the n-body propagation SNOPT consistently converges to a limited number of minuma that are superior to the majority of MSQP solutions.
(a)
(b) Figure 9. Optimal ΔV predicted by the Hill’s equations and full-ephemeris propagation versus tf . (a) Õ Öt× km (b) Õ -1401 km
A variety of locally optimal transfers to orbits of the same energy, in the sense of Eq. 8, produced by the methodology outlined above are shown in Fig. 10. If the transfer time is long enough such that the CW equations lose significant accuracy SNOPT will still generally converge on a solution.
(a)
(c)
(b)
(d)
Figure 10. Full-ephemeris transfers found using SNOPT. (a) Õ Öt× km, (b) Õ -1401 km, (c) Õ Öt× km, (d) Õ ttÖ ØÙ
In the plane, long duration transfers can be found that take advantage of weakly chaotic orbits found near stable distant retrograde orbits.5 The plots in Fig. 11 demonstrate that these dynamics can be exploited in full-ephemeris models. The transfers represent incremental changes of the relative orbit that would be necessary while the spacecraft performs reconnaissance of the body.
(a)
(b)
Figure 10. Full-ephemeris planar transfers. (a) Õ tÎÎ km, ¸¹ =3.3 m/s, ¶· =3 day (b) Õ ÖÎÎ km, ¸¹ =9.0 m/s, ¶· =6.5 day
CONCLUSION
The analytical analysis provided a qualitative description of the design space for higher fidelity models and a prediction of the evolution of ΔV with parameters such as transfer time. Phasing and energy constraints related to the properties of Hill’s equations reduced the dimension of the problem by two greatly simplifying the analysis. For realizable transfer times the analysis provided precise guesses to seed the Hill’s and full-ephemeris models, and an accurate estimation of ΔV. A simple algorithm provides a solution to the assumption of time invariance of the science orbits. This algorithm in conjunction with multiple shooting and robust SQP algorithms such as SNOPT provides a reliable means to connect theoretical developments with highly accurate propagations suitable for detailed mission design. This reliability facilitates further automation of the design process and the evaluation of numerous mission designs. The introduction of additional constraints, such as tangential departure or insertion, or, if operationally feasible, a third maneuver could yield improvements. REFERENCES 1
“SNOPT: An SQP Algorithm for Large-Scale Constrained Optimization”, SIAM Review, Vol. 47, No. 1, pp. 99-131, 2005.
2
Pavlak, T., “Mission Design Applications in the Earth-Moon System: Stationkeeping”, M.S. Thesis, Purdue University, 2010.
Transfer Trajectories and
3
Lovell, T. A., Tragesser, S. G., and Tollefson, M. V., “A Practical Guidance Methodology for Relative Motion of LEO Spacecraft Based on the Clohessy-Wiltshire Equations,” AAS Paper 04-252, Feb. 2004. 4
Montenbruck,O., “Numerical Integration Methods for Orbital Motion,” Celestial Mechanics and Dynamical Astronomy, Vol. 53, pp. 59-69, 1992.
5
Scott, C.J. and Spencer, D.B., “Transfers to Sticky Distant Retrograde Orbits”, Journal of Guidance, Control, and Dynamics, Vol. 33, No. 6, pp. 1940-1946, 2010.
