1. Introduction
The transfer between two circular, coplanar Keplerian orbits of a spacecraft equipped with a continuous thrust propulsion system is a classic astrodynamics problem [
1] that is typically solved in an optimal framework by maximizing (or minimizing) a given scalar performance index as, for example, the propellant consumption or the total flight time. In this context, the optimization procedure described in the fundamental textbook by Bryson and Ho [
2] is a typical example.
In fact, when the optimization problem is approached using the well-known indirect method [
3,
4], i.e., a method based on the classical theory of the calculus of variations [
5], the solution requires the evaluation of the initial value of a set of costate variables, which is usually obtained by numerically solving a two-point boundary value problem (TPBVP). The computation time required by the numerical solution (and also by the convergence of the entire procedure) of the TPBVP associated with the optimization procedure is strongly dependent on the guess of unknown initial costates [
6,
7].
The purpose of this paper is to propose an analytical procedure to accurately approximate the initial costate variables in a coplanar, circle-to-circle minimum-time transfer, in order to initialize the numerical procedure that solves the TPBVP with an initial guess that is close to the actual solution. In particular, this paper considers a mission scenario in which the thrust vector is freely steerable during flight, and the magnitude of the propulsive acceleration varies over a finite range. The proposed procedure, which is consistent with the results from the classical literature [
8,
9,
10], is suitable for the case of a propulsive acceleration vector with a low (or very low) maximum magnitude. In this scenario, the minimum-time transfer trajectory can be approximated through a tight spiral [
11,
12], and some characteristics of the spacecraft osculating orbit can be described in an analytical (simple) form.
The proposed procedure is used here to analyze the optimal transfer in a set of both heliocentric and geocentric mission scenarios, as a function of the propulsive acceleration maximum magnitude. In particular, optimal circle-to-circle (rapid) spacecraft trajectories are studied in a classical LEO-GEO transfer and in a rendezvous mission towards Venus, Mars, Jupiter, and comet 29P/Schwassmann–Wachmann 1. Finally, the method discussed in this paper can also be used to rapidly obtain a set of accurate information about the optimal transfer trajectory in order to initialize more refined (trajectory) design procedures, such as those recently proposed by Huang et al. [
13].
A natural extension of the method proposed in this work is to consider the dynamics of a spacecraft, which include the change in mass due to propellant consumption. This important aspect of the thrust vector model introduces a further parameter in the (optimal) design of the mission, namely the specific impulse, the value of which influences the overall transfer performance. On the other hand, the proposed method can be adapted to a scenario based on the use of a propellantless propulsion system such as, for example, a scenario involving a photonic solar sail [
14,
15] or Janhunen’s Electric Solar Wind Sail [
16,
17].
2. Mathematical Preliminaries
Consider a two-dimensional scenario in which a spacecraft
S initially (time instant
) covers a circular, Keplerian orbit of assigned radius
around a primary body of gravitational parameter
μ. The spacecraft is equipped with a continuous-thrust propulsion system which gives a freely steerable thrust vector and a propulsive acceleration
whose magnitude
ranges in the interval
, where
is the maximum (finite) value of
. The continuous propulsive acceleration induced by the primary propulsion system is used to transfer the spacecraft to a circular, coplanar, target orbit of assigned radius
by minimizing the required flight time
, where
is the final time instant. In other terms, the spacecraft transfer mission coincides with a minimum-time, circle-to-circle orbit raising if
, or orbit lowering when
. Bearing in mind the symmetry of the transfer problem, introduce a classical polar reference frame
with the origin in the primary body center-of-mass
O, in which
r is the
O-
S distance and
is the polar angle measured from the
O-
S line at the initial time instant. The polar reference frame
is sketched in
Figure 1, where
(or
) is the radial (or transverse) unit vector.
In
Figure 1, the direction of the spacecraft propulsive acceleration vector is identified by the thrust angle
, which is defined as the angle (measured counterclockwise into the plane of the circular parking orbit) between the
O-
S line and the direction of
. Accordingly, the propulsive acceleration vector can be written as a function of the (control) angle
as
so that the polar form of the spacecraft equations of motion are the well-known
where
u (or
v) is the radial (or transverse) component of the spacecraft velocity. The system of differential Equation (
2) is completed by the 4 initial conditions that model the spacecraft traveling in the circular parking orbit, viz.
If the final value of the spacecraft polar angle
is unconstrained, the conditions for a circle-to-circle orbit transfer give the following 3 scalar constraints at the final time instant
The spacecraft equations of motion (
2) and the boundary constraints given by Equations (
3) and (
4) can be rewritten more conveniently in a dimensionless form by introducing the (dimensionless) variables
so that Equation (
2) becomes
where the prime symbol indicates the derivative with respect to the dimensionless time
. Accordingly, the boundary conditions written in Equations (
3) and (
4) become
and
in which
is the dimensionless target radius, and
is the dimensionless flight time to be minimized through a suitable selection of the spacecraft controls
. Note that Equations (
6)–(
8) are independent of both the actual value of the primary body gravitational parameter and the radius of the circular parking orbit. The optimization of the spacecraft transfer trajectory, i.e., the calculation of the time variation in both the thrust angle
and the propulsive acceleration (dimensionless) magnitude
that minimizes the flight time, is described in the next section.
3. Trajectory Optimization Using an Indirect Method
The circle-to-circle, minimum-time continuous-thrust orbit transfer optimization is studied using an indirect method, that is, a method based on the classical calculus of variations. In this context, and keeping in mind the mathematical model described in references [
2,
5], the dimensionless performance index
J to be maximized is written as
while the (dimensionless) Hamiltonian function
is given, according to Equation (
6), by
where
are the dimensionless costates, and
is the part of the Hamiltonian function that explicitly depends on the two control terms
, viz.
The optimal value of the controls
is obtained by using Pontryagin’s maximum principle [
18], i.e., by maximizing (at any time instant) the function
defined in the previous equation. More precisely, bearing in mind that
is linear in
and that the thrust vector direction (i.e., the value of
) is unconstrained, the optimal control law is simply given by
where
represents the dimensionless form of the maximum propulsive acceleration magnitude
. Note that, as expected, the control law reassumed in Equation (
12) is consistent with the result presented in reference [
2]. In particular, the first of Equation (
12) indicates that the rapid transfer is obtained by selecting the maximum value of
(i.e., no coasting arcs appear in the optimal transfer trajectory), while the last two equations indicate that the optimal thrust direction coincides with the classical Lawden’s primer vector direction [
19].
In this context, the time variation in the optimal thrust angle is obtained, through Equation (
12), by (numerically) solving the Euler–Lagrange equations [
2,
5], which give the
-derivative of the four costates
as
The expression of the Hamiltonian function is also used to evaluate the transversality condition [
2], which in this case gives the following two scalar constraints at the final time
In particular, taking Equation (14) into account, we obtain that
during the transfer, so Equations (
13) and (16) are simplified as
On the other hand, observing that the dimensionless time
does not compare explicitly in the expression of the Hamiltonian function given by Equation (
10), one has that the value of
is a constant of motion [
2]. Accordingly, the second of Equation (
17) can be rewritten as
Bearing in mind Equations (
10) and (
11), the optimal control law given by Equation (
12), the initial conditions reassumed in Equation (
7), and the expression of
given by Equation (
18), the previous equation gives the following constraint at the initial time instant
where
(or
) is the initial value of the dimensionless costate
(or
), which is one of the unknowns of the associated TPBVP. The previous equation states that the point
lies on a circle of radius equal to
, as sketched in
Figure 2, where the auxiliary angle
is introduced.
Accordingly, Equation (
22) and
Figure 2 allow the pair
to be calculated as a function of the single (geometric) dimensionless variable
as
Note that the auxiliary angle
coincides with the initial value of the thrust angle, that is,
; compare Equation (
12) with Equation (
23).
The TPBVP associated with the optimization problem is then formed by the 7 differential equations, Equations (
6), (
15), (
19) and (
20), the 4 initial conditions (
7), and the 3 final conditions (
8). In particular, the solution of the TPBVP requires the final time
, the auxiliary angle
, and the initial costate
to be numerically calculated by enforcing the (final) boundary constraints of Equation (
8). However, the convergence of the procedure used to solve the TPBVP, which is usually based on standard numerical methods as the single or multiple shooting procedure [
20], depends on the appropriate selection of a guessed value of the triplet
. In this respect, the next section describes an analytical method for selecting a triplet
sufficiently close to the effective value that resolves the TPBVP.
4. Analytical Method for the Approximation of Flight Time and Initial Costates
In this section, an analytical approximation of the three unknowns , that is, the flight time and the initial costates, is detailed. In particular, the proposed procedure can be applied when . This is the important case of a spacecraft equipped with a low-performance (or a very low-performance) propulsion system, which gives a propulsive acceleration whose maximum value is only a (small) fraction of the gravitational acceleration along the circular parking orbit. In this scenario, the total flight time is usually very high (a few dozen, or even hundreds, of parking orbit periods), and the knowledge of an accurate guess of the triplet is a crucial point in reducing the computation time required to solve the associated TPBVP. In this context, the last part of this section describes a method that can be used to rapidly estimate, as a function of the mission scenario characteristics in terms of and , the maximum value of compatible with the proposed, analytical approximate procedure.
In fact, the proposed method is based on the observation that, when
, the optimal transfer trajectory resembles a tight spiral, i.e., the value of the spacecraft (dimensionless) radial acceleration
is nearly zero during the entire transfer between the two (coplanar) circular orbits. In this case, the numerical simulations indicate that the direction of the propulsive acceleration vector is substantially transverse (i.e., the radial component of
is very low when compared with the transverse component) so that the thrust angle
is substantially a constant of motion with
Note that the latter equation can be rewritten in a compact form as
where
is the signum function. Accordingly, the value of the first of the three unknowns, i.e., the auxiliary angle
which coincides with the initial value of the thrust angle, is approximated through the simple equation
Recall that the value of
allows the two initial costates
and
to be calculated using Equation (
23) as a function of the maximum magnitude of the dimensionless propulsive acceleration
.
The second unknown, i.e., the value of the dimensionless flight time
, can be approximated using a simplified form of the spacecraft dynamics. In fact, bearing in mind that
during the entire transfer (see first of Equation (
12)), the last of Equation (
6) gives
where
is the dimensionless magnitude of the spacecraft-specific angular momentum vector
, i.e.,
. Furthermore, according to the third of Equation (
6), when
and
, the dimensionless magnitude
can be approximated as
so that, with the aid of the previous equation, one has
, and Equation (
27) can be easily solved to give an analytical approximation of the dimensionless (minimum) flight time
Note that the expression of the minimum flight time given by Equation (
29) is consistent with the classical result obtained by Alfano et al. [
8,
9,
10] using the elegant concept of the accumulated velocity change. The same concept has also been used by the author to evaluate the optimal performance of a solar-sail-based spacecraft in a typical heliocentric mission scenario [
21,
22].
The last of the three unknowns required to initialize the associated TPBVP, i.e., the value of the initial costate
, can be approximated by firstly enforcing in Equation (15) the initial conditions given by Equation (
7), viz.
Now, according to Equations (
23) and (
26), the value of
is given by
while
is substantially constant during the transfer (remember that the direction of the propulsive acceleration vector is essentially transverse during the transfer) so that
The value of the last unknown
is then obtained by substituting Equations (
31) and (
32) in Equation (
30), and the result is
or, equivalently,
.
Maximum Value of Compatible with the Proposed Method
The approximation of the triplet
given by Equations (
26), (
29) and (
33) is accurate when the optimal transfer trajectory resembles a tight spiral so that the osculating orbit is substantially circular and the radial component of the spacecraft acceleration
is nearly zero. The latter is, in fact, a scenario analyzed in detail by Alfano et al. [
8]. In particular, reference [
8] discusses that the accumulated velocity change
defined as
in this specific case (i.e., when the osculating orbit is nearly circular) can be approximated through a very simple equation that depends on the radius of both the parking and the target orbit, viz.
In particular, the latter equation indicates that the accumulated velocity change
coincides with the difference between the spacecraft velocity along the target and the parking orbit. Note also that Equation (
35) is consistent with the model proposed in this paper; see, for instance, Equation (
29).
Using the results detailed in reference [
8] in terms of plots of the accumulated velocity change
as a function of the characteristics of the propulsion system
and the radius of the target circular orbit, we obtain that Equation (
35) holds (i.e., the approximation of an osculating nearly circular orbit is valid) when the spacecraft completes, during the transfer, a number
n of revolutions around the primary body at least equal to 5. However, extensive numerical simulations indicate that the method proposed in this work still gives an accurate approximation of the initial costates when
. In fact, as underlined by Alfano et al. [
8], when
, the intensity of the propulsive acceleration becomes comparable to the gravitational acceleration of the primary body, and the shape of the transfer trajectory (i.e., the characteristics of the spacecraft osculating orbit) is substantially influenced by the thrust vector magnitude and direction.
The value of
n can be estimated using an analytical expression obtained through the combination of Equations (
6), (
25), and (
28). Indeed, assuming
during the transfer, with the aid of Equation (
28), the second and the last of Equation (
6) give
which can be easily solved to obtain the final value of the spacecraft polar angle
as a function of both the (dimensionless) maximum propulsive acceleration magnitude
and the target orbit radius
, viz.
Accordingly, the value of
n is obtained from the previous equation as
where
is the floor function [
20].
Therefore, the procedure can be summarized as follows. For a given (dimensionless) value of the target orbit radius
and the maximum propulsive acceleration magnitude
, the estimated number of complete revolutions during the transfer
n is calculated through Equation (
38). If the obtained value is
, the proposed analytical method is valid, and the expression of the triplet
given by Equations (
26), (
29), and (
33) can be used to obtain a reasonable approximation of the actual values that resolve the associated TPBVP. Otherwise, the actual value of the unknown initial costates can be very different from those determined analytically. On the other hand, Equations (
26), (
29), and (
33) can still be used to rapidly obtain a (very rough) guess to initialize the TPBVP numerical solution through standard methods [
20].
5. Model Validation and Numerical Simulations
The proposed procedure was tested in a set of circle-to-circle orbit, two-dimensional transfers that model some typical heliocentric (in which the primary body is the Sun and ) and geocentric (in which the primary body is Earth and ) mission scenarios.
In particular, regarding a potential heliocentric mission scenario, the initial parking orbit radius was set equal to . In fact, the case of models the situation in which the spacecraft begins the heliocentric phase of a typical interplanetary transfer after escaping from Earth using a parabolic (escape) orbit. In this context, three different (circular) target obits were considered; that is, three different values of were used in the numerical simulations, namely (i) , which models an Earth–Venus transfer; (ii) , which models an Earth–Mars transfer; and (iii) , which models an Earth–Jupiter transfer. Accordingly, in this selected set of heliocentric mission scenarios, the dimensionless value of the target orbit radius is . As regards the geocentric mission scenario, a classical (two-dimensional) LEO-GEO transfer is considered, in which the radius of the circular parking orbit is (i.e., a circular LEO of of altitude is assumed) and the target orbit has a radius . Therefore, in the geocentric case, the dimensionless value of the target orbit radius is .
The number
n of completed revolutions around the primary body, in the four proposed mission scenarios, can be calculated using Equation (
38) as a function of the value of
. The function
is shown in
Figure 3 when
and
. In particular, the stepped curves in
Figure 3 allow us to quickly estimate the maximum value of
for which the proposed analytical approximation can be used to initialize the associated TPBVP. In fact, according to
Figure 3, the condition
is met when (i)
in the Earth–Venus scenario; (ii)
in the Earth–Mars scenario; (iii)
in the Earth–Jupiter scenario; and (iv)
in the LEO-GEO scenario.
As a result, the validity of the proposed method was tested by considering, for each mission scenario, about 20 values of
ranging between
and
(with a step of
). For each value of
, the optimal (i.e., the rapid) transfer trajectory was calculated by solving the associated TPBVP using a single shooting method [
20], which was initialized according to the expressions given by Equations (
26), (
29), and (
33). In that context, the equations of motion (
6) and the Euler–Lagrange Equations (15), (
19), and (20) were integrated in double precision using the Adams–Bashforth method [
23] with an absolute and relative tolerance of
, while the tolerance on the convergence of the TPBVP was set to
. The solution of the numerical procedure gives the actual value of the triplet
, which was compared with the analytical approximation given by Equations (
26), (
29), and (
33) by introducing the ratio
defined as
where
indicates the value of □ calculated through the analytical approximation, while
is the actual value of □ (calculated numerically) that solves the TPBVP.
The results of the numerical simulations for the four exemplary mission scenarios are summarized in
Table 1,
Table 2,
Table 3 and
Table 4 in terms of
as a function of
. In particular, as also explicitly indicated in the table captions, the values of
are the results of the numerical solutions of the associated TPBVP, the value of
n is obtained analytically from Equation (
38), and the three ratios
are calculated according to Equation (
39). The last three columns in
Table 1,
Table 2,
Table 3 and
Table 4 confirm that the analytical approximations proposed in this work are consistent with the actual (numerical) values of the unknown costates. Furthermore, the expression of
n given by Equation (
38) gives a very good approximation of the actual number of completed revolutions around the primary body during the transfer. Finally, the value of
in all the transfers analyzed indicates that the expression of the flight time given by Equation (
29) is well suited to rapidly calculate the minimum transfer time without solving (numerically) the optimization problem, especially for low propulsive acceleration magnitude.
Case Study: Transfer towards Comet 29P/Schwassmann–Wachmann 1
The proposed method is now used to rapidly solve the TPBVP associated with the (heliocentric) minimum-time transfer towards comet 29P/Schwassmann–Wachmann 1 [
24,
25,
26,
27]. The latter was very recently considered by the author [
28] as a potential target for a solar-sail-based (rendezvous) interplanetary mission. The interested reader can refer to the Introduction section of reference [
28] to obtain a snapshot of the recent research concerning comet 29P/Schwassmann–Wachmann 1 (indicated as “comet SW1” in the rest of this paper). The current, heliocentric orbit of comet SW1 is nearly circular (the actual eccentricity is
), with an inclination of about
. The comet orbit is between Jupiter and Saturn [
29]. The numerical results discussed in reference [
28] indicate that a simplified circle-to-circle, two-dimensional transfer scenario can be used to approximate the real, three-dimensional Earth–comet SW1 trajectory in a classical interplanetary rendezvous mission. In this context, the comet orbit around the Sun is considered circular, with radius
, and coplanar to the heliocentric (assumed circular with radius
) orbit of Earth.
Therefore, in this heliocentric scenario (in which
), the dimensionless target radius is
, and Equation (
38) gives the function
, sketched in
Figure 4. According to
Figure 4, the condition
is approximately satisfied when
, so the optimal Earth–comet SW1 transfer is studied parametrically considering (again) a dimensionless maximum propulsive acceleration in the range
. Note that, in this case, one has
so that a value
(or
) corresponds to a dimensional propulsive acceleration of
(or
).
The TPBVP associated with the optimization problem was solved using the procedure previously described, and the results of the numerical simulations are summarized in
Table 5, which can be considered as an extension (to this scenario) of the results of
Table 1,
Table 2,
Table 3 and
Table 4. As expected, the values in
Table 5 confirm the accuracy of the proposed analytical model based on Equations (
26), (
29), and (
33).
The time variation in the thrust angle and the corresponding optimal transfer trajectory are two outputs of the optimization process. Assuming, for example,
, the polar form of the optimal transfer trajectory is sketched in
Figure 5, while the (dimensionless) time variation in the thrust angle
is shown in
Figure 6. Note that, according to
Figure 6, when the magnitude of the propulsive acceleration is sufficiently low (as in the case of
), the value of the thrust angle during the entire transfer is close to
, as modeled by Equation (
25). However, when the value of
increases (or the primary body–spacecraft distance is high, as at the end of the transfer), the thrust angle diverges from the value predicted by Equation (
25). In any case, the initial value of
, which coincides with the auxiliary angle
, is roughly equal to
, as expected.
6. Conclusions
In this paper, we studied the transfer between two circular, coplanar orbits of assigned characteristics, considering a spacecraft operated with a continuous thrust propulsion system that provides a magnitude of propulsive acceleration vector within an assigned range. The problem is addressed in an optimal framework using a classical indirect approach, minimizing the total flight time as a function of the target orbit radius and the characteristics of the propulsion system. The main contribution of this work is to use a simplified version of the spacecraft dynamics, not to approximate the (optimal) transfer trajectory but to accurately estimate the initial value of the unknown costates needed to solve the TPBVP associated with the optimization process.
The proposed procedure allows the initialization of the TPBVP with a guessed solution reasonably close to the actual one determined through a numerical approach based on the classic single shooting method. The method is designed to work well when the maximum magnitude of the propulsive acceleration vector is sufficiently small when compared with the gravitational acceleration of the primary body along the circular parking orbit. However, the procedure can still be used when the value of the propulsive acceleration is medium–high, in order to obtain a rough estimate of the initial costates. In this regard, the availability of a set of analytical equations, such as the one presented in this paper, is a useful (and important) tool in the preliminary phase of mission design.