Abstract
A complete first-order analytical solution, which includes the short periodic terms, for the problem of optimal low-thrust limited-power transfers between arbitrary elliptic coplanar orbits in a Newtonian central gravity field is obtained through canonical transformation theory. The optimization problem is formulated as a Mayer problem of optimal control theory with Cartesian elements—position and velocity vectors—as state variables. After applying the Pontryagin maximum principle and determining the maximum Hamiltonian, classical orbital elements are introduced through a Mathieu transformation. The short periodic terms are then eliminated from the maximum Hamiltonian through an infinitesimal canonical transformation built through Hori method. Closed-form analytical solutions are obtained for the average canonical system by solving the Hamilton-Jacobi equation through separation of variables technique. For transfers between close orbits a simplified solution is straightforwardly derived by linearizing the new Hamiltonian and the generating function obtained through Hori method.
1. Introduction
In the general analysis of optimal space trajectories, two idealized propulsion models are of most frequent use [1]: LP and CEV systems. In the power-limited variable ejection velocity systems—LP systems—the only constraint concerns the power, that is, there exists an upper constant limit for the power. In the constant ejection velocity limited thrust systems—CEV systems—the magnitude of the thrust acceleration is bounded. In both cases, it is usually assumed that the thrust direction is unconstrained. The utility of these idealized models is that the results obtained from them provide good insight about more realistic problems.
The purpose of this work is to present a complete first-order analytical solution, which includes the short periodic terms, for the problem of optimal low-thrust limited-power transfers (no rendezvous) between arbitrary elliptic coplanar orbits in a Newtonian central gravity field. This analysis has been motivated by the renewed interest in the use of low-thrust propulsion systems in space missions verified in the last two decades.
The optimization problem is formulated as a Mayer problem of optimal control theory with Cartesian elements—position and velocity vectors—as state variables [1]. After applying the Pontryagin Maximum Principle and determining the maximum Hamiltonian, classical orbital elements are introduced through a canonical Mathieu transformation. The short periodic terms are then eliminated from the maximum Hamiltonian through an infinitesimal canonical transformation built through Hori method—a perturbation canonical method based on Lie series [2]. The new maximum Hamiltonian, resulting from the infinitesimal canonical transformation, describes the extremal trajectories associated with the long duration maneuvers.
Closed-form analytical solutions are obtained for the average canonical system governed by this new Hamiltonian. The separation of variables technique is applied to solve the Hamilton-Jacobi equation associated to the average canonical system [3]. For long duration maneuvers, the existence of conjugate points is investigated through the Jacobi condition and the envelope of extremals is obtained considering different configurations of initial and final orbits. Curves representing maneuvers with same fuel consumption for a specified duration, called isocost curves, are presented for several configurations of initial and final orbits. An iterative algorithm based on the first-order analytical solution is described for solving the two-point boundary value problem of going from an initial orbit to a final orbit.
For transfers between close coplanar elliptic orbits, a simplified solution is straightforwardly derived by linearizing the new Hamiltonian and the generating function obtained through Hori method. Since this solution is given by a linear system of algebraic equations in the initial value of the adjoint variables, the two-point boundary value problem is solved by simple techniques.
Some of the results presented in the paper are very similar to the ones obtained by Edelbaum [4] and Marec and Vinh [5] through different approaches based on the concept of mean Hamiltonian. We note that the theory described in the paper includes directly the periodic terms such that it can be applied to any transfer independently of its duration. On the other hand, since the theory is based on canonical transformations, a first-order solution for the Cartesian elements can also be easily obtained from Lie theorem [2].
2. Optimal Space Trajectories
A low-thrust limited-power propulsion system, or LP system, is characterized by low-thrust acceleration level and high specific impulse [1]. The ratio between the maximum thrust acceleration and the gravity acceleration on the ground, , is between and . For such system, the fuel consumption is described by the variable defined as where is the magnitude of the thrust acceleration vector , used as control variable. The consumption variable is a monotonic decreasing function of the mass of the space vehicle, where is the maximum power and is the initial mass. The minimization of the final value of the fuel consumption is equivalent to the maximization of .
Let us consider the motion of a space vehicle , powered by a limited-power engine in a Newtonian central gravity field. At time , the state of the vehicle is defined by the position vector , the velocity vector and the consumption variable . This set of variables provides a straightforward physical interpretation to the optimal thrust acceleration.
The problem of general optimal low-thrust limited-power transfer can be formulated as a Mayer problem of optimal control as follows [1]. It is proposed to transfer the space vehicle from the initial state at time to the final state at time , such that the final consumption variable is a minimum (Figure 1).
For coplanar transfer problem, the initial and final conditions are specified in terms of orbital elements—semimajor axis, eccentricity, and pericenter argument—introduced in the next section. It is also assumed that and the position of the vehicle in the initial orbit are specified. The state equations are where is the gravitational parameter. It is assumed that the thrust direction and magnitude are unconstrained [1].
According to the Pontryagin Maximum Principle [6], the optimal thrust acceleration must be selected from the admissible controls such that the Hamiltonian reaches its maximum. The Hamiltonian is formed using (2.3), where , and are the adjoint variables and dot denotes the dot product. Since the optimization problem is unconstrained, it follows that The optimal thrust acceleration is modulated [1]. Then, the optimal trajectories are governed by the maximum Hamiltonian , Since is an ignorable variable, it is determined by simple integration, and its adjoint is a first integral, whose value is obtained from the transversality conditions: . Thus,
This Hamiltonian can be decomposed into two parts: where is the undisturbed Hamiltonian and is the disturbing function concerning the optimal thrust acceleration. The undisturbed Hamiltonian plays a fundamental role in the canonical transformation theory described in the next section.
3. Transformation from Cartesian Elements to Orbital Elements
In order to obtain a first-order analytical solution for transfers between arbitrary elliptic coplanar orbits, the set of classical orbital elements is introduced through a canonical transformation by using the properties of generalized canonical systems [7].
Consider the extended two-body problem defined by the canonical system of differential equations governed by the undisturbed Hamiltonian , that is, the problem defined by the differential equations for the state variables and , which describe the classical two-body problem, and the differential equations for the adjoint variables and . These equations are given by where is the unit vector pointing radially outward of the moving frame of reference.
The general solution of the differential equations for the state variables and is well-known from the classical two-body problem [8]. In the coplanar case, it is given by where is the unit vector along circumferential direction of the moving frame of reference; is the semimajor axis, is the eccentricity, and is the true anomaly. The unit vectors and are written as function of the orbital elements as follows: where and are the unit vectors of a fixed frame of reference (Figure 1) and is the pericenter argument.
The general solution of the differential equations for the adjoint variables and is obtained by computing the inverse of the Jacobian matrix of the point transformation between the Cartesian elements and the orbital ones, defined by (3.2) and (3.3). This matrix is obtained through the variations of the orbital elements induced by the variations in the Cartesian elements, as described in the next paragraphs.
Let us consider the inverse of the point transformation defined by (3.2) and (3.3): where the eccentricity vector and the angular momentum vector , shown in Figure 2, are given, as function of the Cartesian elements and , by the following equations: Here, the symbol denotes the cross product. The true anomaly has been replaced by eccentric anomaly . These anomalies are related through the equation
Now, consider the variations in the Cartesian elements, and , given in the moving frame of reference by
The variations of the orbital elements—, , , and —induced by the variations in the Cartesian elements, and , are obtained straightforwardly from (3.5) through (3.8) and are given by where the variations of the vectors and are written as being the unit vector normal to the plane of the orbit. Here, and denote the radial and circumferential components of the velocity vector, respectively (see (3.3)).
From (3.2), (3.3), and (3.9) through (3.14), one gets the explicit form of the variations of the orbital elements—, , , and —induced by the variations in the Cartesian elements, and . At this point, it is useful to replace the variation in the eccentric anomaly by the variation in the mean anomaly , obtained from the well-known Kepler equation This variation is given by Thus, where is the mean motion.
Equation (3.17) can be put in the form where the matrix is inverse Jacobian matrix of the point transformation between the Cartesian elements and the orbital ones, defined by (3.2) and (3.3).
Following the properties of generalized canonical systems [7], the general solution of the differential equations for the adjoint variables and is given by with and expressed in the moving frame of reference. Thus, it follows from (3.17) through (3.19) that , , and so forth are known functions of the elliptic motion, which can be written explicitly in terms of eccentricity and mean anomaly through Lagrange series [8].
Equations (3.2), (3.3), and (3.20) define a canonical Mathieu transformation between the Cartesian elements and the orbital ones . The Hamiltonian is invariant with respect to this canonical transformation. Thus,
The new Hamiltonian function , defined by (3.21) and (3.22), describes the optimal low-thrust limited-power trajectories in a Newtonian central gravity field for transfers between arbitrary elliptic coplanar orbits. It should be noted that becomes singular for circular orbits and, thus, the theory developed in the paper cannot be applied to transfers involving orbits with very small eccentricities.
A first-order analytical solution for the canonical system described by the Hamiltonian function will be derived in the next sections through canonical transformation theory.
4. A First-Order Analytical Solution
In order to get an approximate formal solution to the problem of transfers between arbitrary elliptic coplanar orbits, a perturbation technique based on Lie series—Hori method—will be applied. For completeness, a brief description of Hori method [2] is presented in the appendix.
In this study, for low-thrust propulsion systems, it is assumed that the functions and have different order of magnitude: has zero order and has first order in a small parameter associated to the magnitude of the thrust acceleration.
Consider an infinitesimal canonical transformation built through Hori method, Following the algorithm described in the appendix, at order 0 one gets
Now, consider the undisturbed canonical system described by : general solution of which is given by The subscript 0 denotes the constants of integration.
Introducing the general solution defined by (4.4) in the equation of order 1 of the algorithm (A.6), one gets The functions and are written in terms of the constants of integration. According to (A.11), the mean value of must be calculated. After lengthy calculations, one gets Terms factored by have been omitted in equations above, since only transfers (no rendezvous) are considered [4]. But, it should be noted that mixed secular terms can arise in the new Hamiltonian in the case of rendezvous maneuvers. For a first-order solution, and .
The new Hamiltonian governs the optimal long duration transfers between arbitrary elliptic coplanar orbits. The general solution of the canonical system described by this Hamiltonian will be obtained through two canonical transformations.
First, consider the Mathieu transformation: defined by the following equations: The Hamiltonian function is invariant with respect to this transformation. Thus,
The canonical system governed by the new Hamiltonian function has three first integrals ( and are ignorable variables), (Constant should not be confused with the eccentric anomaly .) These first integrals play important role in determining the general solution of the canonical system governed by the new Hamiltonian function through Hamilton-Jacobi theory.
Let us consider the problem of determining a canonical transformation defined by a generating function such that and become the new generalized coordinates, Since the new Hamiltonian function is a quadratic form in the conjugate momenta, the separation of variables technique will be applied for solving the Hamilton-Jacobi equation [3].
Consider the transformation equations where is the generating function. The corresponding Hamilton-Jacobi equation is then given by Following the separation of variables technique, it is assumed that the generating function is equal to a sum of functions, each of which depends on a single old variable, that is,
Therefore, from (4.11) through (4.17), it follows that A complete solution of these equations is given by with . We note that is given as indefinite integral, since only its partial derivatives are needed (4.15), as shown in the next paragraphs.
Now, consider the differential equations for the conjugate momenta of the canonical system governed by the new Hamiltonian , that is, whose solution is very simple: where , , are arbitrary constants of integration.
Introducing the generating function , defined by (4.17) and (4.19), into the transformation equation (4.15), and taking into account the general solution defined by (4.21) for the conjugate momenta, one gets
From equations above, one finds the solution of the canonical system governed by the Hamiltonian for a given set of initial conditions: with the auxiliary constants , and defined as functions of the initial value of the adjoint variables (conjugate momenta) by The constants , , and in (4.22) can also be written as functions of the initial value of the adjoint variables as follows: The initial conditions for the state variables (generalized coordinates) are given by , and , and is obtained from .
Equation (4.23) through (4.30) represents the solution of the canonical system concerning the problem of optimal long duration low-thrust limited-power transfers between arbitrary elliptic coplanar orbits in a Newtonian central field.
For maneuvers with arbitrary duration, the periodic terms must be included [9]. Following Hori method and applying the initial conditions, one gets a first-order analytical solution: with given by (4.9). As mentioned before, these equations cannot be applied to transfers involving very small eccentricities.
A complete first-order analytical solution for the consumption variable is obtained through a different procedure as described in [9]. This procedure is based on the linearized solution for transfers between close orbits and it involves the generating function . The expression of can be put in a compact form as follows: with .
Equations (4.33) are in agreement with the ones obtained in [9] through Bogoliubov-Mitropolsky method, a noncanonical perturbation technique. But it should be noted that in [9] terms in appear instead of terms in ; the other terms are equivalent. The eccentric anomaly is computed from Kepler's equation with the mean anomaly given by where . The differential equation for the mean anomaly is derived from the undisturbed Hamiltonian and an additional term factored by , where dots denote the terms defined by (4.6).
For transfers between coaxial orbits, with no change in the pericenter argument, the analytical solution described in the preceding paragraphs simplifies. Since , some equations must be rewritten: with and .
Figures 3, 4, and 5 show a comparison between the complete first-order analytical solution, the average solution, and a numerical solution. The numerical solution has been obtained through integration of the set of canonical equations which describes the optimal trajectories (system of differential equations governed by new Hamiltonian function , defined by (3.21)–(3.22)). Two sets of initial conditions and transfer duration defined in Table 1 (in canonical units) are used in the comparison. In Table 2, final values of state variables are shown. A Runge-Kutta-Fehlberg method of orders 4 and 5, with step-size control, relative error tolerance of , and absolute error tolerance of , as described in [10, 11], has been used in the numerical integrations. Note that there exists an excellent agreement between the complete analytical solution and the numerical one. Figures 3 and 4 show that the difference between these solutions is lesser than for semimajor axis and eccentricity; but the analytical solution for pericenter argument loses accuracy for large values of time. Figure 5 shows the results for the consumption variable. On the other hand, by comparing the analytical or numerical solution with the average solution, we see that the amplitudes of the short periodic terms are small, but they can be significant for transfers with moderate duration. Note that the difference between the numerical and average solutions decrease as the transfer duration increases. The good agreement between the numerical and the complete analytical solution suggest that the latter can be used in the solution of the two-point boundary value problem of going from an initial orbit to a final orbit, as described in Section 6.
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
5. Analysis of Long Duration Transfers
In this section, a complete analysis of long duration transfers described by the Hamiltonian is presented. The investigation of conjugate points and the solution of the two-point boundary value problem are considered.
To perform this analysis, a new consumption variable is introduced [4], , and (4.23) and (4.24) can be put in the form where . By eliminating from these equations, one finds On the other hand, , thus Note that in equations above.
Now, consider (4.24)–(4.27) which define a two-parameter family of extremals in the phase space for a given initial phase point corresponding to an initial orbit. By eliminating the auxiliary variables and , and can be written as explicit functions of , that is, and . The conjugate points to the phase point are determined through the roots of the equation [12] Since does not depend on , . On the other hand, from (4.24) and (5.2), one finds Thus, the problem of determining the conjugate points reduces to the analysis of the roots of the following equation:
From (4.26) and (4.27), one finds the explicit form of , that is, The partial derivative is given by The auxiliary variable is reintroduced to simplify the expression.
Therefore, from (5.7) and (5.9), after some simplifications, it follows that the conjugate points are given by the roots of following equation [4]: As discussed in [13, 14], the conjugate points determined through (5.10) occur for transfers during which the direction of motion in the orbit reverses; that is, for transfers between direct orbits (no direction reversals) the extremals computed for the canonical system governed by the Hamiltonian are globally minimizing.
Some numerical results are presented through Figures 6–9. Figures 6 and 7 show the field of extremals and isocost curves for long duration transfers between coaxial orbits, considering two cases: in the first one, the space vehicle departs from a circular orbit and, in the second one, it departs from an elliptic orbit with eccentricity . In both cases, the semimajor axis is , such that the results are presented in canonical units. The extremals are plotted for different values of the constant and the isocost curves for different values of the new consumption variable . Only transfers between direct orbits are considered.
(a)
(b)
Figures 8 and 9 show the field of extremals and isocost curves for long duration transfers between noncoaxial orbits. Note that in the general case the family of extremals is described by two constants and . In the plots of Figure 8, has been chosen in order to illustrate the general aspect of the field of extremals and compare the results with the ones obtained in [4, 5, 13, 14]. The projection of extremals on the -plane shows clearly the nonexistence of conjugate points for the extremals computed from (4.6)–(4.11) in the case of transfers between direct orbits (no direction reversals).
Figure 9 shows a sample of isocost curves for different values of the new consumption variable and a fixed value of . Note that these curves lie on an isocost surface with complex geometry. For simplicity, the symmetric interval for is considered in this representation.
To complete the analysis, a brief description of an iterative algorithm for solving the two-point boundary value problem of going from an initial orbit to a final orbit is presented. The steps of this algorithm can be described as follows.(1) Guess a starting value of ,(2) Determine and through (4.9) and (4.26).(3) Compute through (4.27).(4) If , adjust the value of and repeat steps (2) and (3) until , where is a prescribed small quantity.(5) Compute and through (4.25) and (5.3).(6) Compute using equation , where (7) Compute successively , and through the following equations: The adjust of the value of can be obtained through the well-known Newton-Raphson method, that is, with the partial derivative computed from (5.9).
In the case of coaxial orbits, and the solution of the two-point boundary value problem involves only the solution of (4.37), which is given through step (6) with replaced by . The initial value of the adjoint variables is given as follows: is computed through step (7) and is given by
With the solution of the two-point boundary value problem, the optimal thrust acceleration required to perform the maneuver can be determined through (2.5) and (3.22), that is, The state and adjoint variables are computed from (4.23) through (4.30). For simplicity, prime denoting the average variables is omitted.
6. Solution of the Two-point Boundary Value Problem
In this section, an iterative algorithm based on the approximate analytical solution defined by (4.33) is briefly described for solving the two-point boundary value problem of going from an initial orbit to a final orbit at the prescribed final time .
For a given time, (4.33) can be represented as follows: where , and . Note that , and appear explicitly in the short periodic terms, but they also appear implicitly through , , and . Thus, functions , , are nonlinear in these variables.
So, the two-point boundary value problem can be stated as follows: find , and such that , and . This problem can be solved through a Newton-Raphson algorithm [10, 11] which uses the algorithm described in the preceding section (for long duration transfers) to generate a starting approximation of the adjoint variables.
Figure 10 shows a comparison between the solution of the two-point boundary value problem based on the average and complete solutions for a transfer between coaxial orbits with the departure orbit defined by the set of initial conditions , , and the arrival orbit defined by the set of terminal conditions , , , and two transfer durations and . All variables are represented in canonical units.
(a)
(b)
From results above, we note that average solution of the TPBVP does not represents a mean solution of the complete solution of the same TPBVP; that is, these solutions have different values for the initial adjoint variables, and, the solution of the TPBVP based on the average solution becomes closer to the solution based on the complete analytical one, as the transfer duration increases.
7. Transfers between Close Elliptical Coplanar Orbits
For transfers between close elliptical coplanar orbits, an approximate solution of the system of differential equations governed by the Hamiltonian , defined through (3.21) and (3.22), can be straightforwardly obtained by linearizing the functions and about a reference elliptic orbit defined by a nominal set of orbital elements , and . This solution is very similar to the complete solution given by (4.33), with , and replacing , , , and can be put in a suitable matrix form where denotes the imposed changes on orbital elements (state variable), , , is the matrix of initial values of the adjoint variable and is symmetric matrix. The adjoint variables are constant, and the matrix is given by where where is calculated by solving Kepler’s equation for , with . Here, is the time of pericenter passage. The overbar denotes the reference orbit.
Equations (7.1)–(7.3) represent an approximate first-order solution for optimal low-thrust limited-power transfers between close elliptic coplanar orbits. The initial value of the adjoint variables must be determined to satisfy the two-point boundary value problem of going from an initial orbit to a final orbit at the prescribed final time . In this case, the vector can be determined by simple techniques since (7.1) defines a linear system of algebraic equations [15–17]; no iterative method such as Newton-Raphson is needed.
Figure 11 shows a comparison between the complete (nonlinear) first-order solution defined by (4.33) and the approximate (linearized) first-order solution defined by (7.1)–(7.3) for three set of terminal orbits described in Table 3, with and 100 canonical units. Note that linearized solution provides a good approximation for the nonlinear first-order solution. This approximation is better for transfers with larger duration. For transfers with moderate duration corresponding to some revolutions, the amplitude of periodic terms decreases as the eccentricities of terminal orbits increase.
(a)
(b)
(c)
(d)
(e)
(f)
8. Conclusion
A first-order analytical solution for the problem of optimal low-thrust limited-power transfers between arbitrary elliptic coplanar orbits in a Newtonian central gravity field has been obtained through Hamilton-Jacobi theory and a perturbation method based on Lie series. A complete analysis of long duration transfers, which includes the investigation of conjugate points and the solution of the two-point boundary value problem, is presented. Results show the good agreement between the complete analytical solution and the one obtained through numerical integration of the set of canonical equations which describes the optimal trajectories for a transfer between coaxial orbits. On the other hand, the average solution provides a good approximation for large values of the time when compared to the complete analytical solution. For transfers between close coplanar elliptic orbits, a linearized solution is straightforwardly obtained by linearizing the new Hamiltonian function and the generating function about a reference elliptic orbit . Numerical results show the good agreement between the linearized solution and the complete (nonlinear) analytical solution for such transfers. Finally, we note that a first-order solution for the Cartesian elements can also be easily obtained from Lie theorem [2], since the theory presented in the paper is based on canonical transformations.
Appendix
Hori method is based on Lie theorem [2].
Let be a Hamiltonian written as a power series of a small parameter : where the undisturbed Hamiltonian describes an integrable system. Here, and are -vectors of generalized coordinates and conjugate momenta, respectively.
Consider an infinitesimal canonical transformation defined by a generating function , also developed as a power series of the small parameter , This transformation is such that the new dynamical system has some advantages for the solution. Here, prime denotes the new variables, and and are also -vectors. The new Hamiltonian , resulting from this canonical transformation, is also developed as a power series of the small parameter ,
The transformation equations are given explicitly by where , , , , denotes an arbitrary function and braces stand for Poisson brackets.
According to the algorithm of Hori method, the new Hamiltonian and the generating function are obtained, at each order of the small parameter , from the equations:(i)order 0: (ii)order 1: (iii)order 2: All functions in above equations are written in terms of the new set of canonical variables .The th order equation of the algorithm can be put in the form where is a function obtained from the preceding orders, involving , , , and , .
The determination of the functions and is based on the general solution of the undisturbed system and it is performed through an averaging principle. Following the integration theory proposed in [18], (A.8) can be put in the form with written in terms of the arbitrary constants of integration, , of the general solution of the canonical system (A.9). and are unknown functions.
Following Hori [2], suppose that the canonical transformation generated by the function is such that the time is eliminated from the new Hamiltonian . This is accomplished taking as the mean value of the function . Therefore, functions and are given through the equations where stands for the mean value of the function. It should be noted that the averaging process and the integration in (A.11) are performed considering the explicit dependence on time of the functions.
Acknowledgment
This research has been supported by CNPq under contract 305049/2006-2.