Research Article  Open Access
Sandro da Silva Fernandes, Francisco das Chagas Carvalho, "A FirstOrder Analytical Theory for Optimal LowThrust LimitedPower Transfers between Arbitrary Elliptical Coplanar Orbits", Mathematical Problems in Engineering, vol. 2008, Article ID 525930, 30 pages, 2008. https://doi.org/10.1155/2008/525930
A FirstOrder Analytical Theory for Optimal LowThrust LimitedPower Transfers between Arbitrary Elliptical Coplanar Orbits
Abstract
A complete firstorder analytical solution, which includes the short periodic terms, for the problem of optimal lowthrust limitedpower 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. Closedform analytical solutions are obtained for the average canonical system by solving the HamiltonJacobi 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 powerlimited 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 firstorder analytical solution, which includes the short periodic terms, for the problem of optimal lowthrust limitedpower 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 lowthrust 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.
Closedform analytical solutions are obtained for the average canonical system governed by this new Hamiltonian. The separation of variables technique is applied to solve the HamiltonJacobi 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 firstorder analytical solution is described for solving the twopoint 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 twopoint 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 firstorder solution for the Cartesian elements can also be easily obtained from Lie theorem [2].
2. Optimal Space Trajectories
A lowthrust limitedpower propulsion system, or LP system, is characterized by lowthrust 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 limitedpower 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 lowthrust limitedpower 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 firstorder 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 twobody 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 twobody 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 wellknown from the classical twobody 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 wellknown 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 lowthrust limitedpower 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 firstorder analytical solution for the canonical system described by the Hamiltonian function will be derived in the next sections through canonical transformation theory.
4. A FirstOrder 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 lowthrust 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 firstorder 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 HamiltonJacobi 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 HamiltonJacobi equation [3].
Consider the transformation equations where is the generating function. The corresponding HamiltonJacobi 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 lowthrust limitedpower 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 firstorder analytical solution: with given by (4.9). As mentioned before, these equations cannot be applied to transfers involving very small eccentricities.
A complete firstorder 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 BogoliubovMitropolsky 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 firstorder 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 RungeKuttaFehlberg method of orders 4 and 5, with stepsize 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 twopoint 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 twopoint 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 twoparameter 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 twopoint 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 wellknown NewtonRaphson method, that is, with the partial derivative computed from (5.9).
In the case of coaxial orbits, and the solution of the twopoint 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 twopoint 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 Twopoint 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 twopoint 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 twopoint boundary value problem can be stated as follows: find , and such that , and . This problem can be solved through a NewtonRaphson 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 twopoint 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 firstorder solution for optimal lowthrust limitedpower transfers between close elliptic coplanar orbits. The initial value of the adjoint variables must be determined to satisfy the twopoint 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 NewtonRaphson is needed.
Figure 11 shows a comparison between the complete (nonlinear) firstorder solution defined by (4.33) and the approximate (linearized) firstorder 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 firstorder 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 firstorder analytical solution for the problem of optimal lowthrust limitedpower transfers between arbitrary elliptic coplanar orbits in a Newtonian central gravity field has been obtained through HamiltonJacobi 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 twopoint 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 firstorder 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