Orbital transfers: optimization methods and recent results

A wide variety of techniques have been employed in the past for 
optimizing orbital transfers, which represent the trajectories that 
lead a spacecraft from a given initial orbit to a specified final 
orbit. This paper describes several original approaches to 
optimizing impulsive and finite--thrust orbital transfers, and 
presents some very recent results. First, impulsive transfers 
between Keplerian trajectories are considered. A new, analytical 
optimization method applied to these transfers leads to conclusions 
of a global nature for transfers involving both ellipses and escape 
trajectories, without any limitation on the number of impulses, and 
with possible constraints on the radius of closest approach and 
greatest recession from the attracting body. A direct optimization 
technique, termed direct collocation with nonlinear programming 
algorithm, is then applied to finite--thrust transfers between 
circular orbits. Lastly, low--thrust orbital transfers are optimized 
through the joint use of the necessary conditions for optimality and 
of the recently introduced heuristic method referred to as particle 
swarm optimization. This work offers a complete description and 
demonstrates the effectiveness of the distinct techniques applied to 
optimizing orbital transfer problems of different nature.


1.
Introduction. On January 2, 1959, the Soviet Luna 2 spacecraft was the first space vehicle to perform orbital maneuvers, that were midcourse trajectory corrections while approaching the Moon. In the succeeding decades orbital maneuvering became a common task either for satellites orbiting the Earth or for spacecraft dedicated to planetary exploration.
As the vehicle weight is a crucial issue for space missions, the minimization of propellant consumption for performing a specified orbital maneuver is desirable. Hence, the optimization of orbital maneuvers has the main objective of minimizing propellant expenditure and has been a subject of great interest for many decades in the last century. The related studies started since 1920s with the pioneering work by Hohmann [40], and continued in the subsequent decades. Significant advances are dated back to the 1950s, when also the modern optimal control theory began developing, by gradually assuming its current form, due to the researches of some eminent scientists, such as Bliss [9], Leitmann [61,62], Cicala [14], Belmann [62,5], Miele [62,66], Pontryagin [62,77], Bryson [12], and Vinh [92]. In the same period Lawden conceived and introduced the primer vector theory [59], which is concerned with the application of the first-order necessary conditions for optimality (arising from the calculus of variation) to space trajectories. Conditions on the primer, which represents the adjoint variable conjugate to the velocity vector, can be used to describe an optimal trajectory, either with finite-thrust periods or velocity impulses. A velocity impulse approximates the effect of the propulsive thrust, by assuming that an instantaneous velocity change can occur (while the spacecraft position remains unchanged). The impulsive approximation has application in the case of high-thrust rockets and spacecraft and represents a very accurate approximation [81,94] of the finite-thrust maneuvers that a vehicle performs once it is in orbit, under the assumption that gravitational losses during propulsion can be neglected. When applied to impulsive trajectories, the primer vector yields the direction and the position of the impulsive changes of velocity, which must occur tangentially at the apse points of a Keplerian arc. Lawden [59,60] conjectured also the optimality of a specific type of singular, finite-thrust trajectory, the so-called Lawden's spiral. However, if no constraint holds for the thrust level, the Lawden's spiral turns out to be non-optimal, as demonstrated by several researchers [57,80,4]. This means that a locally optimal trajectory is unequivocally composed of coasting (null-thrust, Keplerian) arcs, separated by impulsive changes of velocity, which are to be applied tangentially at the apse points. Definitely, these arguments suffice to prove that any finite-thrust trajectory can be equaled or outperformed by an impulsive trajectory. Significant contributions to the investigation on impulsive trajectories were provided by Hoelker and Silber [39], Shternfeld [85], Edelbaum [23], Ting [88,89], Lawden [59], Barrar [2], Moyer [68], Marec [63], Battin [3], Palmore [70], and Prussing [78]. Section 2 of this paper covers impulsive transfers. With reference to transfers between two elliptic orbits, this research is intended to prove that the Hohmann transfer [40] and the bielliptic transfer [39,85] with the maximum apoapse radius of the intermediate orbit are globally optimal. To do this, a method is employed that utilizes ordinary calculus in conjunction with a simple graphical procedure aimed at eliminating sub-optimal impulses. The graphical procedure described in this paper is analogous to that used by Hazelrigg [34]; nevertheless, all the fundamental results attained by Hazelrigg through Green's theorem here are simply deduced by using elementary properties of the functions of a single variable. The final choice between Hohmann and bielliptic transfer depends on the orbital parameters (i.e. on the apoapse and periapse radii) of the two orbits. This study addresses a direct, analytical method to make this choice in the general case of ellipses of arbitrary periapse and apoapse radii. Global optimal transfers involving elliptic orbits and escape trajectories with possible constraints on the radius of closest approach or greatest recession are also investigated.
Of course, spacecraft can employ only a finite thrust level to perform maneuvers. In the last decades several researches have been concerned with the optimization of finite-thrust spacecraft trajectories, which was pursued through either direct [32,26,27,8,84,18], indirect [67,10,11,55,64,79], or hybrid [45,97,30] techniques. Section 3 is focused on finite-thrust transfers between two coplanar circular orbits of specified initial and final radii. The overall trajectory is assumed to be composed of three arcs: (i) an initial thrust arc, where the (maximum) constant thrust is employed, (ii) an intermediate coast arc, with no propulsion, and (iii) a final thrust arc (with maximum thrust again). The optimization has the purpose of minimizing propellant consumption, by choosing the optimal thrust direction time history, and the optimal time intervals associated with each flight phase. This kind of orbital transfer resembles the (impulsive) Hohmann transfer, in the sense that a coast arc separates two propulsive phases. Thus, finite-thrust transfers optimization has also the natural objective of evaluating the performance degradation due to the use of finite-thrust, in comparison with the corresponding optimal impulsive transfer. Different propulsive parameters are considered, in order to evaluate the accuracy of the impulsive approximation as a function of the spacecraft propulsive characteristics. The optimization method employed in this work is the direct collocation with nonlinear programming (DCNLP) algorithm [32,26]. This technique converts the optimization problem of interest into a nonlinear programming problem, through discretization of the state and control variables and the use of (high-order) quadrature rules for the system dynamics equations. This direct approach has been successfully applied to space trajectory optimization problems in the past [36,17,87], and was proven to be reliable and numerically accurate. However, a reasonable guess must be provided for the DCNLP algorithm to converge. This task is here accomplished by employing the results obtained for impulsive orbital transfers.
More and more frequently, space vehicles are equipped with propulsive systems that supply only modest thrust levels. This occurs, for instance, when spacecraft are thrusted by electric propulsion. The main advantage of using these propulsive systems is due to a higher value of their specific impulse (i.e. to a higher value of the effective exhaust velocity of the ejected mass). Optimizing continuous lowthrust orbit transfers is a challenging problem, due to two circumstances: (i) the propulsive thrust exerts a low control authority over a long time interval, and (ii) no suitable guess for (local) optimization methods can be derived from the results attained for impulsive orbital transfers. In fact, the impulsive approximation is not suitable for modeling low-thrust trajectories, especially when continuous thrust is used. In the 1960s Edelbaum [24] proposed some analytical solutions concerning low-thrust transfers. More recently, significant contributions in this field are due to Kluever [56], Kechichian [48,49,50], Ilgen [44], Jenkin [46], and Scheel and Conway [83]. Section 4 considers the optimization of continuous low-thrust transfers between two circular orbits. The problem consists in minimizing the time of flight -and simultaneously propellant consumption, since continuous thrust is employed -by determining the optimal thrust direction time history. To achieve this goal, a recently introduced heuristic technique, referred to as particle swarm optimization (PSO), is utilized. This methodology is inspired by the the dynamics of bird flocks while searching for food [19]. In fact, the swarming algorithm emulates the (stochastic) cooperative behavior of the individuals that compose a swarm. No guess is needed, because the initial population is randomly generated at the first iteration of the process. Each particle is associated with a position vector and a velocity vector at a given iteration and represents a possible solution to the problem. At the end of the process, the best particle (i.e. the best solution with reference to the objective function) is selected. The PSO technique requires translating the optimal control problem into a parameter optimization problem, if possible with a fairly small number of parameters. The use of the necessary conditions for optimality allow representing the control variable through a reduced number of parameters, which are primarily the initial values of the adjoint variables conjugate to the equations of motion. A simple approach for equality constraint handling, which is often challenging when evolutionary techniques are used, is described and employed. This work on orbital transfers is intended to: (i) describe three types of orbit transfer optimization problems that arise in the context of space missions, (ii) propose three completely different, original approaches recently employed to solve these problems, (iii) illustrate the related (analytical or numerical) results, with the final purpose of corroborating the suitability of the solution techniques that the paper addresses, and (iv) discuss the validity of the impulsive approximation in modeling finite-thrust orbit transfer maneuvers.
2. Impulsive orbital transfers. More than 80 years ago Hohmann [40] conjectured that in an inverse square gravitational field the minimum fuel impulsive transfer between coplanar circular orbits is the bitangent elliptic transfer performed through two tangential impulsive changes of velocity, which occur at the apse points of the transfer orbit. Such a transfer was termed the "Hohmann transfer" and for a relatively long time it was assumed to be globally optimal in terms of total velocity change. However, since 1959 Hoelker and Silber [39], Shternfeld [85], and Edelbaum [23] independently proposed the three-impulse bielliptic transfer, which was shown to be more economical than the Hohmann transfer when the initial orbit and the final orbit are sufficiently far apart from each other, so that a midcourse impulse can occur sufficiently far outside the outer orbit. Later, Ting [88] proved the local optimality of the Hohmann transfer and stated that the optimal orientation of the axes of the initial and final (elliptical) orbits is always coaxial and aligned (i.e. with periapses on the same side). In addition, Ting [88] demonstrated that for any orbital transfer at the intersection of two non-coplanar orbits there exists a corresponding coplanar transfer with lower cost. This circumstance implies that if an orbit is specified only by its specific energy and angular momentum, the minimum cost transfer is unequivocally coplanar, i.e. corresponds to initial, final, and transfer orbits all belonging to the same plane. Lastly, Ting [89] proved also that orbital transfers performed through four or more impulses are never optimal. Proofs of the global optimality of the Hohmann transfer in the class of two-impulse transfers were first supplied by Lawden [59] using the calculus of variations, and Barrar [2]. Then Moyer [68] demonstrated the global optimality of Hohmann and bielliptic transfers from a circular to an elliptic orbit without restrictions on the number of impulses. Still later, Marec [63] employed a hodograph analysis to prove the optimality of the Hohmann transfer, whereas Battin [3] achieved the same result by using Lagrange multipliers. Most recently, Palmore [70] and Prussing [78] presented elementary proofs of the global optimality of the Hohmann transfer in the class of two-impulse transfers. In his valuable paper [34], Hazelrigg described a quite interesting approach to the analysis of global optimal transfers without any limitation on the number of impulses. His method was based on the use of Green's theorem and led to results of a global nature regarding impulsive transfers involving any kind of Keplerian trajectory.
In all of the above researches, as well as in this study, the initial and final trajectories are specified through their specific energy and angular momentum only. The total time required for the transfer is assumed unspecified, as well as the orientation of the initial and final coplanar trajectories (time-open, angle-open transfers). This is the starting point of the present research. The first section of this work describes a simple, original approach proposed by the author [73] with the intent of determining the global optimal impulsive transfers between two arbitrary coplanar Keplerian trajectories, without any restriction on the number of impulses and with possible constraints on the radius of closest approach or greatest recession (i.e., with a minimum allowed periapse radius and a maximum allowed apoapse radius).  [73], whereas for parabolic trajectories it is trivial because on parabolic trajectories the velocity vanishes as the distance from the attracting body tends to infinity. Lastly, the statement (3) will be proven in this section, by considering the limiting case of an impulse at periapse as the periapse radius tends to zero. Of course, the application of an impulse at periapse of a rectilinear trajectory is concretely infeasible because gravity tends to infinity. However, the statement (3) refers to a limiting situation and means that arbitrarily modest impulses at periapse (capable of changing the trajectory energy) can be applied if no constraint holds on the radius of closest approach. This section is concerned with a convenient representation of Keplerian trajectories and locally optimal transfers (performed through impulsive variations of velocity) as a starting point for all the subsequent considerations leading to the identification of global optimal transfers. In addition, this section provides the expressions for evaluating the magnitude of the impulses that perform a locally optimal transfer.
A conic trajectory is defined through its periapse radius and through the parameter = 1/ , where: • for elliptic orbits is the apoapse radius • for circular orbits = ⇒ = 1 • for parabolic trajectories → ∞ and therefore = 0 • for hyperbolic orbits < 0 and | | represents the periapse distance from the vacant focus As a result, in the ( , )-plane, portrayed in Fig. 1, each point corresponds to a specific conic. Similarly to the previous work made by Hazelrigg [34], optimal orbital transfers can be conveniently represented in this plane. Due to property (1), in the ( , )-plane a sequence of impulsive transfers that satisfy the necessary conditions for (local) optimality may be represented as a cascade of connected horizontal and vertical segments. Moreover, in the same plane, due to property (2), an infinitesimal impulse at infinite distance is associated with a hyperbolic arc (as the energy remains unchanged): where denotes the gravitational parameter of the attracting body. Different regions in the ( , )-plane correspond to distinct types of Keplerian trajectories: • the region > 0 corresponds to elliptic orbits • the axis = 0 corresponds to parabolic trajectories • the region < 0 corresponds to hyperbolic trajectories • the axis = 0 corresponds to rectilinear trajectories, of elliptic ( > 0), parabolic ( = 0), and hyperbolic type ( < 0) • in the region ( > 0, > 0), the hyperbolic arc = 1/ corresponds to circular orbits and represents the boundary of all the feasible states • in the region ( > 0, < 0), the hyperbolic arc = −1/ represents the boundary of all the feasible states, because the condition + 1 = 0 yields an infinite value for the energy Figure 1. Representation of the Keplerian trajectories in the ( , )-plane Impulsive transfers that satisfy the necessary conditions for (local) optimality are represented as: • horizontal segments (impulses at apoapse) and vertical segments (impulses at periapse) in the region ( > 0, > 0, ≤ 1) (elliptic and circular orbits) • vertical segments (impulses at periapse) and hyperbolic arcs (infinitesimal impulses at infinite distance, characterized by a constant value for the energy and associated with the relationship (1)) in the region ( > 0, < 0, > −1) (hyperbolic trajectories) • horizontal segments (infinitesimal impulses at infinite distance, characterized by = 0) on the line = 0 (parabolic trajectories) • vertical segments (infinitesimal impulses at periapse) on the line = 0 (rectilinear trajectories) • starting from the line = 0: horizontal segments if > 0 (impulses at apoapse of a rectilinear trajectory of elliptic type), hyperbolic arcs if < 0 (infinitesimal impulses at infinite distance on a rectilinear trajectory of hyperbolic type) Using the vis-viva equation arising from conservation of energy, it is straightforward to obtain the magnitudes of the tangential impulses to enter or depart a trajectory: 1. impulses at periapse modify the apoapse radius; the impulse Δ that changes the initial apoapse radius (= 1/ ) to the value (= 1/ ) has the following magnitude: 2. impulses at apoapse modify the periapse radius; the impulse Δ that changes the initial periapse radius to the value has the following magnitude: Each velocity change is parallel to the instantaneous velocity vector and has the same direction if the orbit energy is to be increased or the opposite direction if the orbit energy is to be reduced. In the region > 0, from state I to state F a possible sequence of ten impulses that satisfy the necessary conditions for (local) optimality is portrayed in Fig. 2(a). In the region of hyperbolic trajectories ( < 0), impulses at periapse modify , and the relationship (2) holds again. In contrast, Eq. (3) is meaningless because no apoapse exists on a hyperbolic orbit. As previously remarked, among hyperbolic trajectories optimal transfers are thus represented as vertical segments (impulses at periapse) and hyperbolic arcs (infinitesimal impulses at infinite distance). Two possible orbital transfers satisfying the necessary conditions for (local) optimality are portrayed in Fig. 2(b): (i) transfer ABCD is performed through two impulses at periapse, divided by an infinitesimal impulse at infinite distance: as the first impulse at infinity must be unequivocally performed while leaving from the attracting body, the second cannot be performed, and the sequence ABCD turns out to be infeasible; (ii) transfer LMNO is performed through two infinitesimal impulses divided by an impulse at periapse: the first impulse is performed at infinity while approaching the attracting body, the third is performed at infinity while leaving from it. This means that orbital transfers can involve at most two infinitesimal impulses at infinite distance. Transfers similar to that associated to path ABCD of Fig. 2 Three basic properties concerning (locally) optimal transfers can be easily derived. First, using Eqs.
i.e. the impulse required for changing the initial state (related to the initial orbit) into the final state (related to the final orbit) is the same if the two states are interchanged, provided that the orbital transfers are locally optimal. This circumstance implies that the total cost of an optimal orbital transfer can be evaluated regardless of which is the initial and which is the final state. Second, the following additive properties hold: where < 1 < or < 1 < and < 1 < or < 1 < . This means that, in in terms of total velocity change, any orbital transfer performed by applying two subsequent tangential impulses at periapse (apoapse) is equivalent to the application of a single impulse. By induction, one can deduce that a vertical (horizontal) line can represent a sequence of an arbitrary number of subsequent impulses at periapse (apoapse), which sequence is equivalent to applying a single impulse. Finally, from Eq. (2) one obtains: i.e. in the limit as → 0, the velocity impulse required for changing to the value vanishes. This means that on a rectilinear trajectory an arbitrary change of the orbit energy can be achieved by means of an infinitesimal impulse at periapse, represented as a vertical segment lying on the axis = 0. 2.2. Globally optimal transfers between elliptic orbits. To determine the fundamental features of the global optimal impulsive transfers between two elliptic orbits without any restriction on the number of impulses, globally optimal transfers in the class of two-impulse transfers must first be identified. This goal can be achieved by applying a simple analytical procedure based on ordinary calculus, which leads to the same results achieved by Hazelrigg [34] through Green's theorem.
2.2.1. Two-impulse globally optimal transfers. In the region ( > 0, > 0), including all the possible elliptic orbits, four mutual positions exist for the initial and the final states, corresponding to the terminal (i.e. initial and final) orbits, and represented as points I and F in Fig. 3. However, as the total cost of an optimal orbital transfer can be evaluated regardless of which is the initial and which is the final state, only the cases portrayed in Figures 3(a) and 3(b) are to be examined. 1. Case A ( Fig. 3(a)): > , > : I ≡ B and F ≡ D. The two paths IAF and ICF are compared, to prove that Δ > Δ . Using Eqs. (2)-(3), one obtains: where What will be demonstrated is that both terms included in the square brackets of Eq. (7) are positive. To do this, the first derivatives of 1 ( ) and 2 ( ) with respect to and are considered: Eqs. (9) become Then, As > and > , the relationships (12) imply: and, due to Eq. (13), one obtains: and these inequalities finally lead to the desired result: Hence, the optimal path from I to F is composed of two segments: IC and CF, the former corresponding to an impulse at apoapse, the latter corresponding to an impulse at periapse. In addition, due to the previous considerations about interchangeability of I and F, the optimal path for the case portrayed in Fig. 3(d) is composed of the segments IC and CF (corresponding to an impulse at periapse followed by an impulse at apoapse). where Also in this case, what is to be demonstrated is that both terms included in the square brackets of Eq. (16) are positive. The same steps employed for case A lead to the desired result, i.e.: Hence, the optimal path from I to F is composed of two segments: IB and BF, the former corresponding to an impulse at apoapse, the latter corresponding to an impulse at periapse. In addition, due to the previous considerations about interchangeability of I and F, the optimal path for the case portrayed in Fig. 3(c) is composed of the segments IB and BF (corresponding to an impulse at periapse followed by an impulse at apoapse). Figures 4(b) and 4(c) illustrate the optimal transfers related to the cases of Figures 3(b) and 3(c), respectively. This analysis suffices to identify the globally optimal impulsive transfer in the class of two-impulse orbital transfers: the Hohmann transfer. In addition, it emerges that the terminal ellipses are optimally oriented when they are coaxial and aligned (i.e. with periapses on the same side), as illustrated in Fig. 4. This result is consistent with the considerations made by Ting [88], Hazelrigg [34], and Winn [93].

2.2.2.
Globally optimal transfers with arbitrary number of impulses. The above analysis is extremely useful for achieving conclusions of a global nature regarding impulsive transfers without any restriction on the number of impulses. In fact, with reference to an arbitrary sequence of (local) optimal impulses, such as the one portrayed in Fig. 5, a simple graphical method allows the selection of the global optimal path from I to F. The elimination of sub-optimal paths in Fig. 5 can be made by taking into account the above results in the ( , )-plane: • from 1 to 3, path 123 is outperformed by the path 1A3, so 123 is eliminated • starting from 3, 4 is reached passing through A, but A is reached directly from 1, so path 123A can be eliminated • from A to 5, path A45 is outperformed by path AB5, so A45 is eliminated • starting from 5, 6 is reached passing through B, but B is reached directly from A, so path A45B can be eliminated • from B to 7, path B67 is outperformed by path BC7, so B67 is eliminated • from 7 to 9, path 789 is outperformed by path 7D9, so 789 is eliminated • from C to D, path C7D is outperformed by path CED, so C7D is eliminated Finally, the optimal transfer is composed of five consecutive impulses at periapse (I1, 1A, AB, BC, CE) followed by three consecutive impulses at apoapse (ED, D9, 9F). However, due to the additive properties (5), in terms of total velocity change the first five impulses are equivalent to applying a single impulse at periapse (IE), whereas the latter three impulses are equivalent to applying a single impulse at apoapse (EF). Therefore, this straightforward graphical method, which can be employed for any arbitrary sequence of (locally optimal) impulses in the shaded region of Fig ≥ min{ , }. In fact, the same graphical procedure can be applied to the impulses represented in Fig. 6 (where denotes the minimum allowed value for , corresponding to the maximum allowed value of the apoapse radius, ( ) ): • from 1 to 3, path 123 is outperformed by path 1A3, so 123 is eliminated • starting from 3, 4 is reached passing through A, but A is reached directly from 1, so path 123A can be eliminated • from A to 5, path A45 is outperformed by path AB5, so A45 is eliminated • from 6 to 7, path 6C7 is outperformed by path 6D7, so 6C7 is eliminated Definitely, the optimal transfer is performed through three consecutive impulses at periapse (I1, 1A, AB) followed by three consecutive impulses at apoapse (B5, 56, 6D) and, lastly, by two impulses at periapse (D7, 7F). However, due to the additive properties (7) and (8), in terms of total velocity change such a transfer is equivalent to applying three impulses: the first at periapse (represented as IB), the second at apoapse (represented as BD), the third at periapse again (represented as DF). Hence, the method leads to the identification of an alternative class of possible global optimal transfers: the three-impulse bielliptic transfers; each of them corresponds to a specific value of , which is related to the apoapse radius of the intermediate orbit.
is constrained to the following range: The biparabolic transfer is the limiting case of the bielliptic transfer as the midcourse impulse location tends to infinity, i.e. = = 0. The total velocity change corresponding to each bielliptic transfer, also referred to as the characteristic velocity of the transfer, is the sum of the following three terms: 1. Impulse at periapse IB 2. Impulse at apoapse BD 3. Impulse at periapse DF It is worth mentioning that Hohmann and biparabolic transfers can be found as special cases from the above expressions (19)-(21): (i) the Hohmann transfer corresponds to = min{ , } (and either Δ ( , ) = 0 or Δ ( , ) = 0); (ii) the biparabolic transfer corresponds to = 0 (and Δ ( , ) = 0, i.e. the midcourse impulse becomes of infinitesimal size). Eqs. (19)-(21) identify a specific class of transfers including Hohmann and all the bielliptic transfers. The graphical method employed for determining the optimal path in Fig. 6 yields the same result for any value of in the range ≤ ≤ min{ , }. Hence, the method at hand is not able to make a comparison between the total costs corresponding to the paths IBDF and IEF of Fig. 6. This is due to the fact that the inequality Δ ( , ) + Δ ( , ) + Δ ( , ) < Δ ( , ) cannot be established through the above results on two-impulse transfers. As a consequence, the method cannot be used to choose the optimal path between IBDF and IEF. However, this final choice is addressed in the next subsection through an original, general approach capable of identifying the global optimal transfer between two elliptic orbits of arbitrary periapse and apoapse radii.

Hohmann transfer versus bielliptic transfer.
The global optimal transfer between two elliptic orbits belongs to the class of three-impulse transfers whose magnitudes are given by Eqs. (19)- (21). For specified values of , , , and the globally optimal transfer corresponds to the optimal choice of , which is related to the apoapse radius of the intermediate orbit, i.e. to the midcourse impulse location. To address this issue, the following dimensionless auxiliary variables are introduced: .
These variables are constrained to the following ranges: Moreover, the following inequalities hold: Eqs. (23)-(24) imply: Using the above definitions, the relationships (19)- (21) can be rewritten as follows: where Hence, a generic transfer corresponding to the value has the following (normalized) cost : Now, two transfers are compared: the first, denoted with the superscript 1, is the transfer with maximum apoapse radius for the intermediate orbit, i.e. = , whereas the second, denoted with the superscript 2, is a generic three-impulse transfer (with ≤ ≤ min{ , }). The following function is defined as the difference of the two related costs: The inspection of the right-hand side of Eq. (30) reveals that is independent of and : = ( , , ). Now, the following property of represents an immediate consequence of its definition: { < 0 ⇔ transfer 1 has a lower cost than transfer 2 > 0 ⇔ transfer 2 has a lower cost than transfer 1 Therefore, the bielliptic transfer with the maximum apoapse radius of the intermediate orbit is the global optimal transfer if < 0. However, the optimal value of (i.e., of ) remains to be determined if > 0. By definition, this optimal value * minimizes (2) ( , , , ): * = arg min The function (2) is positive definite, and (1) is independent of ; as a result: * = arg max i.e. the optimal value * is such that is maximized wherever > 0. To determine * , the partial derivative ∂ /∂ can be considered. This derivative has two distinct expressions (for < 1 and > 1), both independent of : What will be demonstrated is that the function increases monotonically with respect to (i.e. ∂ /∂ ) for any pair of values ( , ) such that > 0, i.e. for any point of the ( , )-plane where > 0. In the ( , )-plane the equations ∂ /∂ = 0 (cf. Eqs. (33)) are employed to obtain explicitly the functioñ =˜ ( ) in the regions of interest, i.e. for < 1 and > 1. Then, using these expressions (not reported for the sake of brevity), the function ( , , ) is evaluated over the curves ∂ /∂ = 0. In this context, becomes a function of and only: . The ranges where and vary are different in the regions˜ > 1 and˜ < 1: The surfaces corresponding to ( , ) are illustrated in Figures 7(a) and 7(b), which refer to˜ > 1 and˜ < 1, respectively. In both cases turns out to be negative in the regions defined by Eqs. (34). This means that < 0 over the curves ∂ /∂ = 0; hence, the region where > 0 is unequivocally included in the region where ∂ /∂ > 0. As a consequence: i.e., with respect to the function is monotonically increasing for any pair of values ( , ) such that > 0. This property holds regardless of the value of (> 0). As a result, is maximized when assumes its maximum value, i.e. * = min( , ). In summary, the globally optimal value * , which identifies the globally optimal transfer, is * = Therefore, in the ( , )-plane the curves = 0, which can be plotted once is specified, separate the regions where the Hohmann transfer is globally optimal from the regions where the bielliptic transfer with the maximum apoapse radius of the intermediate orbit is globally optimal. The final choice between these two transfers turns out to depend only on three parameters: , min( , ), and , i.e. the ratios / , max{ / , / }, and ( ) / . A schematic illustration of the curves = 0 (represented as solid lines) is portrayed in Fig. 8. With reference to distinct values of , several curves corresponding to = 0 are shown in Figures 9(a) and 9(b), corresponding to the cases > 1 and < 1, respectively. For each curve, the value of is reported in the respective inset. As previously remarked, the biparabolic transfer is the limiting case as the midcourse impulse location tends to infinity, i.e = 0 ⇔ = 0. The method described in this section includes previous results regarding optimal transfers between two circular orbits as special cases. Circular orbits are such that Hence, = / . Two cases must be distinguished: 1.
> , i.e. > 1. As Hence, in the ( , )-plane the geometrical locus related to transfers between circular orbits (with > ) is the hyperbola corresponding to * = 1, which is represented as dashed curve in Fig. 8. If = 0 (i.e. = 0), one obtains the following well known result concerning exterior biparabolic transfers: • the Hohmann transfer is globally optimal if 1 < < , where is the value of at point A of Fig. 9; it follows that ≈ 11.939 • the biparabolic transfer is globally optimal if > Even though not reported for the sake of brevity, an analytical expression of exists, because represents one of the roots of the following third degree equation (cf. [73]): Hence, in the ( , )-plane the geometrical locus related to transfers between circular orbits (with < ) is the segment corresponding to * = 1, which is represented as dashed line (LM) in Fig. 8. If = 0 (i.e. = 0), one obtains the following well known result concerning interior biparabolic transfers: • the Hohmann transfer is globally optimal if < < 1, where is the value of at point B of Fig. 9; it follows that ≈ 0.08376 • the biparabolic transfer is globally optimal if < Also in this case an analytical expression of exists, because represents one of the roots of the following equation (cf. [73]): If the generic 1 is a root of Eq. (38), then the value 1/ 1 is a root of Eq. (40): this property can be easily verified by letting = 1/ 1 in Eq. (40), under the assumption that 1 satisfies Eq. (38). This circumstance implies that the two values and are reciprocal, i.e. = 1.
regardless of the remaining data of the initial and final orbits (provided that > 1). Conversely, ( ) (≈ 15.582) identifies the maximum meaningful value of for which vanishes. This implies that in the region ( ≥ , > 1, ≤ 1) the bielliptic transfer is unequivocally globally optimal if > ( ) , regardless of the remaining data of the initial and final orbits, provided that This result is well known for circular orbits and here is proven to hold for the more general case of elliptic orbits of arbitrary periapse and apoapse radii.
( ) represents one of the roots of the following third degree equation (cf. [73]): In addition, with reference to Fig. 8, the value of at point S, , tends to ( ) in the limiting situation as → 1. It follows that ( ) ≈ 0.06418 and one can conclude that an interior bielliptic transfer is unequivocally globally optimal if < ( ) regardless of the remaining data of the initial and final orbits.  9 it is apparent that the bielliptic transfer is never optimal if < < , where and are the values of at points D and E of Fig. 9(b) and 9(a), respectively. These two values represent respectively the maximum and the minimum value of over the curves ( = 0, , ) = 0, which constitute the geometrical loci where the interior and exterior biparabolic transfers have the same cost as the Hohmann transfer. The two values and can be derived analytically [73] and turn out to be = 1/9 and = 9. In conclusion, the method described for selecting the global optimal transfer can be employed whether

2.3.
Globally optimal transfers between hyperbolic trajectories. As mentioned in Section 2, in the region of hyperbolic trajectories ( < 0) two types of locally optimal impulses exist: • impulses at periapse, capable of changing (and therefore the trajectory energy ) and represented as vertical segments in the ( , )-plane • infinitesimal impulses at infinite distance, corresponding to hyperbolic arcs of constant energy in the ( , )-plane As previously stated, orbital transfers between hyperbolic trajectories can involve at most two infinitesimal impulses at infinite distance, separated by an impulse at periapse, because transfers similar to that corresponding to the sequence ABCD of Fig. 2(b) are infeasible. With reference to Fig. 11, the determination of the global optimal transfer between two hyperbolic trajectories (represented as points I and F) is equivalent to the identification of the optimal impulse at periapse that allows changing the trajectory energy from the value to . This is due to the fact that, starting from I, the points located on the same hyperbolic arc (i.e. H, A, and E in Fig. 11) can be reached by means of infinitesimal impulses (without any additional cost). Hence, to determine the globally optimal transfer from I to F, only the optimal location of the periapse impulse needs to be found. What will be demonstrated is that the globally optimal transfer between two hyperbolic trajectories requires the application of a finite impulse at the minimum allowed periapse radius, denoted with ( ) . To do this, the costs associated with two generic segments AB and EG (portrayed in Fig. 11), such that > , can be compared. From Eq. (2) one obtains: Using Eq. (1), the values of in A, B, E, and G can be expressed as functions of , , , and ( > ≥ 0): To prove that Δ > Δ , the function (Δ − Δ ) / √ 2 is considered. The relationships (44) and (45) yield: The auxiliary function ( ) is now introduced as and Eq. (46) can be rewritten as The derivative d /d is given by: As > , Eqs. (48)- (49) lead to the desired result: Since A and E are arbitrary, Eq. (2.3) proves that the minimum magnitude impulse corresponds to the minimum periapse radius, i.e. to the closest allowed approach to the attracting body. The above considerations lead to the identification of the globally optimal transfer from I to F, which is composed of the following three impulses: 1. One infinitesimal impulse at infinite distance while approaching the attracting body, corresponding to the arc IH 2. One finite impulse at the minimum periapse radius, corresponding to the vertical segment HL 3. One infinitesimal impulse at infinite distance while leaving from the attracting body, corresponding to the arc LF Any other orbital transfer corresponds to a greater cost in terms of total velocity change. Moreover, due to the fact that the total cost of an optimal orbital transfer can be evaluated regardless of which is the initial and which is the final trajectory, the globally optimal transfer from F to I corresponds to the path FLHI. Figure 11. Globally optimal transfer between two hyperbolic trajectories 2.4. Globally optimal transfers involving ellipses and escape trajectories. This section addresses the determination of the globally optimal transfer between an elliptic orbit and an escape trajectory, in the presence of constraints on the radius of closest approach to the attracting body and on the maximum apoapse radius of intermediate orbits: Due to these constraints in the ( , )-plane the shaded regions portrayed in Fig. 12 are not accessible, i.e. no point representing an intermediate trajectory can be located inside these regions. With reference to Fig. 12, to determine the globally optimal path from I to F, corresponding respectively to the initial elliptic orbit and to the final hyperbolic trajectory, the transfers IHF, IEADF, and IELGF are compared. However, due to the additive property (5) and to the fact that the points located on the same hyperbolic arc (i.e. C, G, H, D, and F in Fig. 12) can be reached by means of infinitesimal impulses, the paths EH, EAD, and ELG can be compared instead of IHF, IEADF, and IELGF. Figure 12. Globally optimal transfers from an elliptic orbit to an escape trajectory As a first step, the inequality Δ > Δ will be demonstrated. To do this, the function (Δ − Δ ) / √ 2 is considered: Using Eq. (1), the values of in D and H can be expressed as functions of (> 0), (≥ 0), and (≥ 0): and Eq. (52) becomes where ( ) To prove that (Δ − Δ ) > 0, the derivative (d /d ) is considered. The condition (d /d ) > 0 yields: After squaring and rearranging Eq. (55), one obtains The variable is constrained to the interval [0, 1], due to the fact that 0 ≤ ≤ 1. The inequality (56) is satisfied for any value of because the function on the left-hand side of (56) is greater than 1 ∀ ∈]0, 1], whereas for the right-hand side it follows that: 0 < /( + ) ≤ 1 (the degenerate case = 0 and = 0 corresponds to the transfer to a parabolic trajectory: in this context the paths EH and AD are equivalent). Therefore, the derivative (d /d ) is positive and this circumstance, through Eq. (53), leads to the desired result: Since A is arbitrary, Eq. (2.4) proves that a transfer similar to IEADF is not globally optimal. This means that starting from an ellipse with specified periapse radius , the globally optimal transfer to a hyperbolic trajectory does not include any intermediate elliptic orbit with a periapse radius greater than . Now, the path EH is compared with ELG, by defining the following function:

MAURO PONTANI
where the value of in G is − /( + ). After introducing the expressions of and (written in terms of , , and ) into Eq. (57), one obtains: The auxiliary function ( ) is now defined as and Eq. (58) can be simply rewritten as To establish whether ( The condition (d /d ) > 0 yields: As a consequence, since L is arbitrary, the transfer IEHF (equivalent to IHF) is globally optimal if < . In fact, if the inequality < holds, then from Eqs. (60) and (62) it follows that: < 0, implying that Δ < Δ ∀ L such that < . On the contrary, if > then Δ < Δ (i.e. > 0), and the optimal value of remains to be determined. This optimal value, denoted with * , minimizes Δ . Such a requirement is equivalent to the maximization of , which is defined by Eq. (57), because Δ is independent of and Δ is positive definite: * = arg min Δ ⇒ * = arg max (63) provided that > 0 (i.e. > ). To determine * , the derivative of with respect to is considered, under the assumption that > . From Eqs. (58), (61), and (62) it is evident that Hence, the optimal value * coincides with the minimum allowed value for the periapse radius, i.e. Thus, one concludes that starting from an elliptic orbit two possible globally optimal transfers to a hyperbolic trajectory exist, depending on the values of and : 1. If < , the transfer IHF is globally optimal. It is composed of one impulse at periapse (segment IH), and one infinitesimal impulse at infinite distance while leaving from the attracting body (arc HF) 2. If > , the transfer IEBCF is globally optimal. It is composed of: one impulse at periapse (segment IE), one impulse at the maximum apoapse radius (segment EB), one impulse at the minimum periapse radius (segment BC), and finally one infinitesimal impulse at infinite distance while leaving from the attracting body (arc CF) This result is partially in contrast with that derived by Hazelrigg [34], who identified only the path IEBCF as globally optimum. In addition, the optimal choice between the two options 1 and 2 turns out to be independent of ( ) . Of course, the cases = 0 (no constraint on the radius of greatest recession), ( ) = 0 (no constraint on the radius of closest approach), and = 0 (transfer from an elliptic orbit to a parabolic trajectory) can be deduced as special cases.
As a final remark, due to the interchangeability of I and F (with regard to the total cost of the transfer), the globally optimal paths from F to I are FHI (if < ) and FCBEI (if > ).
3. Finite-Thrust orbital transfers. Space vehicles actually employ finite-thrust to perform orbit transfer maneuvers. Their optimization is aimed again at minimizing propellant consumption. In this section, this goal is achieved by determining the optimal thrust direction time history, and the time intervals associated with each flight phase, under the assumption that the overall transfer trajectory is composed of three arcs: 1. Maximum-thrust arc, in which the thrust level is set to its maximum value, denoted with (time interval [ 0 , 1 ]) 2. Coast (null-thrust) arc, in which no thrust is used (time interval [ 1 , 2 ]) 3. Maximum-thrust arc, in which the thrust level is set to again (time interval [ 2 , 3 ]) The problem of minimizing propellant consumption is formulated as an optimal control problem, where time-dependent continuous variables and time-independent parameters are involved. The time-varying state and control variables include the spacecraft position and velocity, and the thrust-pointing-angle, respectively. Unknown parameters for the problem at hand are represented by the time intervals associated with each flight phase. As finite-thrust transfer problems are not amenable to a closed-form analytical solution, numerical algorithms must be employed.
Several researches have been focused on optimal finite-thrust transfers. Three different types of (deterministic) numerical approaches have been employed to optimize such transfers: • Direct methods, based on converting the optimization problem into a nonlinear programming problem (often involving a large number of parameters), without the use of the analytical necessary conditions for optimality. In this class of methods, relevant contributions are due to Enright and Conway, who successfully employed two distinct types of direct methods: (i) direct collocation with nonlinear programming [32,26] (DCNLP), and (ii) direct transcription with nonlinear programming [27]. This latter technique was employed by Betts [8] also. Differential inclusion has been investigated by Seywald [84] and by Coverstone-Carrol and Williams [18].
• Indirect methods, based on applying the necessary conditions for optimality, which arise from the calculus of variations. Several techniques have been proposed. Miele [67] employed the gradient-restoration algorithm to optimizing aeroassisted orbital transfers. Brown et al. [10] applied a shooting method (in conjunction with some simplifying analytical developments) to generate multiple-burn transfers. A shooting approach was also proposed by Brusch and Vincent [11] and Kluever and Pierson [55], who addressed the optimization of Earth-Moon transfers. McAdoo et al. [64] developed OPBURN, a tool that combines two indirect techniques. Lastly, Redding [79] used an indirect method for optimizing transfers from low Earth orbit to geosynchronous orbit with a large number of burns. • Hybrid methods, which combine direct and indirect techniques, usually by replacing some difficult conditions that arise from the calculus of variations with some simplified relation. Ilgen [45] developed the HYTOP algorithm, which retains the boundary conditions deriving from the calculus of variations, and employs the equinoctial orbital elements. Zondervan et al. [97] proposed a hybrid approach based on directly choosing the time intervals associated with each flight phase, without defining a switching function. Most recently, Kluever [30] used a hybrid algorithm that incorporates multiple shooting.
These three types of techniques belong to the class of "deterministic" methods, as opposed to "stochastic" methods. Due to their theoretical foundations, direct and indirect algorithms possess specific features, which are extensively dealt with in the scientific literature. With regard to convergence, in general direct techniques are more robust, because they are occasionally capable of converging to the desired result even in the presence of a poor guess. Conversely, indirect methods are more numerically accurate and do not need any parametrization of the time-varying quantities involved in the problem of interest. Both direct and indirect algorithms are relatively fast and accurate in yielding an optimal solution if a suitable guess is provided. The implementation of deterministic techniques can require considerable analytical developments, as well as the definition of specific routines for finding some mathematical entities that are involved in the solution process (e.g., the Jacobian matrix of the objective function). However, the two intrinsic features of deterministic approaches that represent their main limitations are: (i) the need of a starting guess; (ii) the locality of results. In fact, both direct and indirect methods are local in nature, in the sense that the optimal solution they are able to detect depends on the available first attempt guess "solution". The numerical result found by deterministic techniques usually lies in the neighborhood of the guess. In this section a simple approach is described to generate a suitable guess for finite-thrust transfers, starting from the results obtained for impulsive transfers. This process for guess generation is possible for the fact that the trajectory structure resembles the Hohmann transfer, in the sense that a coast arc divides two thrust arcs (while in the Hohmann transfer two impulses are applied, with an intermediate coast arc). With a suitable guess available, the DCNLP algorithm is then employed to attain the desired optimal transfer for a variety of cases. This numerical optimization method is chosen because of its numerical accuracy, reliability, and robustness, and also due to the fact that the optimization problem can be translated into a nonlinear programming problem involving an acceptable number of parameters. A detailed description of the DCNLP algorithm is contained in Appendix 1.
This section considers different interior and exterior transfers, corresponding to distinct propulsive parameters and distinct terminal orbits. Finite-thrust optimization has thus the ultimate intent of determining the performance penalty due to the use of finite-thrust in place of impulsive changes of velocity, with reference to different propulsive characteristics of the orbiting spacecraft and distinct terminal orbits. Definitely, this investigation is also aimed at evaluating the accuracy of the impulsive approximation in modeling finite-thrust transfers.
3.1. Problem definition. At the initial time 0 (= 0) the spacecraft is placed in a circular orbit of radius 1 . The final, coplanar, circular orbit has radius 2 . Therefore, if , , and represent respectively the radius, the velocity (magnitude), and the flight path angle (with respect to the local horizontal), the initial conditions (at 0 ) and the final conditions (at ) are given by: where denotes the gravitational parameter of the attracting body. The spacecraft is assumed to employ a constant thrust during both the thrust arcs, and therefore the thrust-to-mass ratio ( / ) has the following expression: (67) whereas ( / ) equals 0 in the cost arc ( 1 ≤ ≤ 2 ). In Eq. (88) 0 is the initial spacecraft mass, and therefore 0 represents the thrust-to-mass ratio at 0 , whereas is the effective exhaust velocity of the propulsive system. It is straightforward to derive the expression for the final-to-initial mass ratio, / 0 : It is worth noticing that Eq. (68) implies that minimizing the term ( 3 − 2 + 1 − 0 ) will also maximize the ratio / 0 , i.e. it will also minimize the propellant required for the transfer maneuver. This means that the following objective function can be assumed for the problem at hand: The spacecraft equations of motion (also referred to as the state equations hence forward) involve the radius , the velocity , and the flight path angle : cos + sin (70) The angle is the thrust-pointing-angle (relative to the velocity direction). The state vector is .
whereas the control vector includes only: . = . During the coast phase propulsion is not employed, and the dynamics equations actually do not need to be integrated. In fact, in [ 1 , 2 ] the trajectory is a Keplerian conic (more specifically, an ellipse), and only a single parameter suffices to define the spacecraft dynamics, i.e. the time evolution of the state components. What will be demonstrated is that the eccentric anomaly variation, denoted with Δℰ, can be conveniently assumed as this parameter. To do this, first the semimajor axis (SMA) and eccentricity of the Keplerian trajectory are to be determined. From the vis viva equation, arising from the conservation of energy, one obtains the SMA : Second, the eccentricity is easily derived from the expression of the (specific) angular momentum ℎ, which is given by ℎ = ( 1 ) ( 1 ) cos[ ( 1 )] and also equals √ (1 − 2 ) along Keplerian elliptic trajectories. As a result: Third, the true anomaly at 1 , 1 , can be calculated by using the expressions for the radius and the radial component of the velocity along an elliptic trajectory: It is then straightforward to obtain the eccentric anomaly at 1 , ℰ 1 , from the true anomaly 1 : cos ℰ 1 = cos 1 + 1 + cos 1 and sin ℰ 1 = sin 1 √ 1 − 2 1 + cos 1 (74) The SMA and the eccentricity remain unchanged in the time interval [ 1 , 2 ]. In contrast, the eccentric anomaly varies and its value at 2 , ℰ 2 , is simply given by ℰ 2 = ℰ 1 +Δℰ. The corresponding true anomaly, 2 , can be determined by employing the counterparts of Eqs. (74): The radius at 2 is whereas the velocity at 2 can be derived again through the vis viva equation: Then, the flight path angle at 2 can be obtained from the expression of the radial component of velocity along an ellipse: Lastly, the coast duration Δ can be easily calculated through Kepler's law: These relations demonstrate that the only parameter Δℰ suffices to derive the state components at 2 (at the beginning of the second thruts arc) from their values at 1 (at the end of the first thrust arc). As discussed in the Introduction section, any finite-thrust transfer can be equaled or outperformed by an impulsive transfer. This circumstance means that the evaluation of the performance associated with the impulsive Hohmann transfer yields the theoretical minimum propellant consumption, i.e. the best performance attainable from a real propulsive system, under the assumption that the ratio 2 / 1 is constrained to the range [0.08376, 11.939] (as discussed in Section 2). The performance penalty due to the use of finite-thrust in place of impulsive changes of velocity can be evaluated by introducing the following function ℘: where ( / 0 ) is given by Eq. (68) and represents the final-to-initial mass ratio when finite-thrust is employed, whereas ( / 0 ) is the same ratio associated with the Hohmann transfer between the same terminal orbits. This latter ratio can be calculated by means of the Tsiolkovsky's law: where Δ 1 and Δ 2 are the two velocity impulses that allow performing the Hohmann transfer. For terminal circular orbits of radii 1 and 2 , it is relatively straightforward to derive the expressions for Δ 1 and Δ 1 from Eqs. (19)- (21): where . = 2 / 1 . Now, the function , related to the final-to-initial mass ratio for a Hohmann transfer, is introduced as and from Eqs. (81)- (82), this function turns out to depend only on and : = ( , ). Some contour lines of this function are portrayed in Fig. 13, for different values of . In general, a minimum admissible ratio ( / 0 ) is to be assumed, due to feasibility of the spacecraft (in the sense that propellant can represent only a portion of the spacecraft, not its entire mass). The curves shown in Fig. 13 are useful to determine the range of admissible values for (in dimensionless units), once , 1 , and ( / 0 ) have been specified. For instance, if = 10, √ / 1 = 1 (in dimensionless units), and ( / 0 ) = 0.7, then the effective exhaust velocity must be greater than ≈ 1.5 (cf. Fig. 13) for the transfer to be feasible.

3.2.
Method of solution. The optimization of finite-thrust orbit transfers is pursued through the use of the DCNLP algorithm, which is a deterministic direct optimization method. Similarly to other deterministic methods, the DCNLP technique requires defining a first attempt reasonable (approximate) "guess" solution to converge to an optimal solution. While in general this task can be challenging, with regard to the problem at hand this work proposes a simple method to generate guess solutions for different terminal orbits. The procedure for guess generation takes advantage of the trajectory structure, which resembles that of a Hohmann transfer. For specified terminal orbits and exhaust velocity , the derivation of Δ 1 and Δ 2 (by means of Eqs. (82)) finally leads to the definition of a guess for finite-thrust transfers. First, as the velocity impulses are optimally oriented when they are applied tangentially at the apse points, the thrust direction for finite-thrust transfers is assumed aligned with the velocity, and therefore in both the thrust arcs = 0 (for exterior transfers, i.e.
Equating these expressions with the values Δ 1 and Δ 2 calculated for the Hohmann transfer leads to obtaining the guess values of Δ 1 and Δ 2 : Lastly, for the parameter Δℰ the following guess is adopted: Δℰ = , because the Hohmann transfer (that approximates the finite-thrust transfer) is characterized by a transfer angle equal to 180 degrees. This step concludes the procedure for guess generation, based on the results attained for impulsive transfers.
Starting from a suitable guess, the DCNLP algorithm converts the optimization problem into a nonlinear programming problem, where the unknown parameters are primarily state and control parameters, and secondarily parameters depending on the application under consideration. The basic features of the DCNLP technique are the following (a detailed description is in Appendix 1): • the continuous problem is discretized in time, so only the values of the state variables and of the control variables at discrete times are employed by the algorithm • the state equations are translated into nonlinear algebraic equations by means of high-order quadrature rules (the highly accurate Gauss-Lobatto fifth-order quadrature rules in this research) • the resulting nonlinear programming problem is solved by a numerical optimizer, which does not utilize explicitly the necessary conditions for optimality arising from the calculus of variations For the transcription of the optimal control problem into a discrete problem each thrust arc is partitioned into 10 subarcs ( 1 = 10 and 2 = 10). Under the assumption that the state vector has components and the control vector has components, the former is represented through (2 ( 1 + 2 ) + 2 ) parameters and the latter through (4 ( 1 + 2 ) + 2 ) parameters. Additional parameters to be optimized are represented by the thrust arcs durations, Δ 1 and Δ 2 , and by the eccentric anomaly variation, Δℰ. As = 3 and = 1 for the problem at hand, the parameters to be optimized are The translation of the optimal control problem into a nonlinear programming problem through collocation yields 2 ( 1 + 2 ) equality constraints, which replace the dynamics equations. The final conditions for orbit insertion, represented by Eqs. (66), constitute additional equality constraints to satisfy.
In the end, the numerical solver has to determine the optimal values of 211 parameters, while holding 2 ( 1 + 2 ) + 3 = 123 equality constraints. The optimal parameters associated with the discrete values of the state are finally employed to interpolate the related components through fifth-degree polynomials (cf. Appendix 1).

3.3.
Numerical results. The problem of optimizing finite-thrust orbit transfers is solved by employing canonical (dimensionless) units: the radius of the initial orbit represents the distance unit (DU), whereas the time unit (TU) is such that = 1DU 3 /TU 2 . As an immediate consequence, the function , defined by Eq. (83), coincides with ( / 0 ) . If the minimum value ( / 0 ) is set to 0.15, the admissible values of (in DU/TU) are located above the two curves portrayed in Fig. 14. This means that for a given value of , there exists a minimum acceptable value of , , which corresponds to a Hohmann transfer with final-to-initial mass ratio equal to 0.15.
In the context of finite-thrust transfer optimization, several values of are considered ( = 0.2, 0.4, 0.6, 0.8, 2, 4, 6, 8, 10). For each value of , some admissible values of are chosen. With reference to Fig. 14, each point corresponds to a finitethrust transfer between two circular orbits with specified values of (= 2 / 1 ) and (in DU/TU). For each of these points four values of the thrust level are considered: 0 = 1, 0.75, 0.5, 0.25 (in DU/TU 2 ). As 20 points are taken into account (cf. Fig. 14), in total 80 finite-thrust orbital transfers have been optimized through the DCNLP algorithm.
The numerical results generated through optimization are reported in Table 1 in terms of the function ℘ (that measures the performance penalty with respect to the Hohmann transfer). Figures 15 through 18 refer to two cases of optimal finite-thrust transfers, and report the related state variables and optimal control laws, as well as  Table 1 it is apparent that in general the performance degradation is modest (never greater than 4.35 × 10 −2 ) and decreases as the thrust level increases and as the value of approaches 1. In most cases the performance penalty is negligible (℘ < 10 −3 for 63 out of 80 cases).   These numerical results definitely prove that the impulsive approximation is highly accurate for a relatively wide range of values of the thrust level 0 and of the effective exhaust velocity , for constrained to [0.2, 10].
4. Low-Thrust orbital transfers. In recent years, innovative propulsive systems have been conceived and realized. Electric propulsion, solar-electric propulsion and solar sails are the most common types of these new, alternative solutions for spacecraft propulsion. Unlike chemical propulsion, they are able to supply only modest thrust levels over a long time interval. While chemical propulsion is still necessary for ascending rockets departing from the Earth surface, low-thrust propulsion can be profitably used for performing space maneuvers, such as orbital transfers. The main advantage of using these propulsive systems is due to a higher value of their specific impulse. As remarked in the Introduction section, optimizing continuous low-thrust orbit transfers is a challenging problem. First, only a low control authority is available and is distributed over a long time interval. Second, the results obtained for impulsive transfers cannot be employed to generate first attempt solutions for studying low-thrust orbit transfers.
Several stochastic methods have been recently introduced that do not require the definition of a guess solution to solve an optimal control problem, unlike deterministic techniques (such as the DCNLP algorithm). Stochastic methods are also referred to as heuristics, meta-heuristics, or evolutionary algorithms, and are inspired by natural phenomena. Evolutionary computation techniques exploit a population of individuals, representing possible solutions to the problem of interest. The optimal solution is sought through cooperation and competition among individuals. The most popular class of these techniques is represented by the genetic algorithms (GA) [31], which model the evolution of a species, based on Darwin's principle of survival of the fittest. Differential evolution algorithms represent alternative stochastic approaches with some analogy with genetic algorithms, in the sense that new individuals are generated from old individuals and are eventually preserved after comparing them with their parents. Ant colony optimization [25] is another method, inspired by the behavior of ants, whereas the simulated annealing algorithm [25] mimics the equilibrium of large numbers of atoms during an annealing process.
The particle swarm optimization (PSO) technique, which is the methodology being addressed in this paper, was first introduced by Eberhart and Kennedy [19,51] in 1995 and belongs to the category of swarm intelligence methods [25,52]. It mimics the unpredictable motion of bird flocks while searching for food, taking advantage of the mechanism of information sharing that affects the overall behavior of a swarm [19,51,52,41,43,42,22,21,13,15,91,71,65,54]. The initial population that composes the swarm is randomly generated at the first iteration of the process. Each particle is associated with a position vector and a velocity vector at a given iteration. More specifically, the position vector includes the values of the unknown parameters of the problem, whereas the velocity vector determines the position update. Each particle represents a possible solution to the problem, and corresponds to a specific value of the objective (or fitness) function. At the end of the process, the best particle (i.e. the best solution with reference to the objective function) is selected. Both the position and the velocity vector are updated in a single iteration. For each particle the formula for velocity update includes three terms with stochastic weights; one of these terms is the so-called social component, related to the collective best position ever visited by a portion of the particles that form the swarm. A detailed description of the PSO algorithm is in Appendix 2.
A number of options are available for implementing the PSO technique. First of all, different values for the stochastic weights have been proposed [21,13,15], and they seem to affect both the algorithm convergence and the capability to detect the global optimum. In addition, two different versions of the particle swarm exist [25,52]: the global version, where the collective best position (associated with the social term in the velocity updating expression) is selected by considering the entire swarm, and the local version, where for each particle the collective best position is selected among the particles located in a proper neighborhood of the particle itself. Albeit less computationally efficient, in the scientific literature [25,52,43,42,22,21,13] the local version of the PSO algorithm, based on the definition of the neighbors of each particle, has been reported to be occasionally capable of avoiding local minima. Additional improvements based on the application of evolutionary operators to the particle swarm methodology are also reported by several researchers [37,47].
The basic version of the particle swarm algorithm appears as very intuitive and is extremely easy to program. Although computationally expensive with respect to gradient-based methods [91], in the scientific literature [1,20,33] the particle swarm technique is reported to be more efficient when compared to genetic algorithms, due to a reduced number of function evaluations. Despite its promising features and the vast number of papers devoted to this technique, most researchers concentrated on topological [28] and multimodal mathematical problems [65,72], and only a limited number of applications appear of practical interest. In the scientific literature several papers describe the use of the particle swarm methodology in the context of chemical processes [16,69]. Fourie and Groenwold [29] employed the PSO technique for shape and size optimization in structural engineering. Most recently Khurana et al. [53] applied the swarming theory to airfoil shape optimization. Bessette and Spencer [6,7] successfully employed the particle swarm technique for space trajectory optimization, focusing on the optimization of time-limited orbital transfers [6] and impulsive interplanetary trajectories [7]. In both papers they claim that the PSO method outperforms the differential evolution algorithm (DE), with regard to convergence speed as well as to reliability. This statement is somewhat in contrast with what is affirmed by other researchers, as Spaans and Mooji [86], who claim the superiority of differential evolution methods with respect to the particle swarm algorithm. Vasile et al. [90] investigated several algorithms for global optimization (including PSO and DE), employing different settings for each of them and using the success rate as the main index for the performance evaluation of each method. In their paper an improved algorithm based on differential evolution is proposed. Nevertheless in the conclusions the authors remark on the importance of appropriately coupling the evolutionary technique with the problem structure, stating that the optimal choice of a method of solution is definitely problem dependent. Zhu et al. [96,95] combined DE and PSO for satellite constellation design [96] and for determining the globally optimal low-thrust trajectory for asteroid exploration [95]. Rosa Sentinella and Casalino [82] proposed a hybrid optimization procedure that runs three different optimizers -based on GA, DE, and PSO -in parallel. They applied this method to the optimization of multiple-impulse rendezvous trajectories and of Earth-to-Mars round trip missions. Most recently, Pontani and Conway [74,75,76] employed a simple version of the PSO algorithm for optimizing impulsive transfers and rendezvous, for the optimization of finite-thrust trajectories, and for the determination of Lyapunov and Halo (periodic) orbits in the circular restricted three-body problem.
This section addresses the optimization of low-thrust orbital transfers between two circular, coplanar orbits by employing a basic, general purpose global version of the PSO algorithm, according to the general settings suggested in References [41,43,42]. This section is intended to show that the PSO methodology, even in its simplest formulation and without interacting with other algorithms, is capable of effectively solving orbit transfer problems with great numerical accuracy.
4.1. Problem definition. The problem of optimizing low-thrust orbital transfers between two coplanar circular orbits consists in determining the thrust direction time history that minimizes the time of flight while satisfying the final conditions for injection into the target orbit.
At the initial time 0 (= 0) the spacecraft is placed in a circular orbit of radius where denotes the gravitational parameter of the attracting body. The spacecraft is assumed to employ a constant low-thrust during the entire time of flight, and therefore the thrust-to-mass ratio ( / ) has the following expression: In Eq. (88) 0 is the initial spacecraft mass, and therefore 0 represents the thrustto-mass ratio at 0 , whereas is the effective exhaust velocity of the propulsive system. It is worth noticing that Eq. (88) implies that minimizing the time of flight will also maximize the remaining mass, i.e. it will also minimize the propellant required.
The spacecraft equations of motion (also referred to as the state equations hence forward) involve the two components of velocity and the radius : The angle is the thrust-pointing-angle (relative to the local horizontal). The state vector is .
where .  (87). The necessary conditions for optimality include the following set of adjoint equations for * (in this section the superscript "*" will denote the optimal value of the respective variable henceforth):˙ * In addition, the optimal control * can be expressed as a function of the costate through the Pontryagin minimum principle: * = arg min ⇒ cos * = − * 2 √ ( * 1 ) 2 + ( 2 1 ) 2 and sin * = − √ ( * 1 ) 2 + ( 2 1 ) 2 (96) Lastly, as the final time is unspecified, the following transversality condition must hold: The necessary conditions for optimality allow translating the optimal control problem into a two-point boundary-value problem involving Eqs. (86) and (89) To do this, one has first to recognize the special structure of the costate equations (93)- (95), which are homogeneous in . This circumstance implies that if an optimization algorithm is capable of finding some initial value of such that then the same proportionality (98) holds between and the optimal * at any , due to homogeneity of Eqs. (93)- (95). Moreover, the control u can be written as a function of through Eq. (96) and one can recognize that u coincides with the optimal control * : This circumstance implies that if the conditions (98) hold then the final conditions (87) are fulfilled at the minimum final time * . In contrast, the transversality condition is violated, because the value of ( * ), due to Eq. (97), turns out to be Therefore, provided that the proportionality condition holds, the optimal control * can be determined without considering the transversality condition, which in fact is ignorable in this context. Thus, this condition is discarded, in order to reduce the number of equality constraints considered by the PSO algorithm, with the intent of improving its performance. Once a costate proportional to the optimal costate * has been determined, it can be suitably scaled. In fact, from Eq. (101) it is apparent that the proportionality coefficient is simply given by = − ( * ). In summary, the method of solution employed for this problem is based on the following points: • the control is expressed as a function of the costate through Eq. isfied to the desired accuracy) have been determined, (= − ( * )) is calculated • the multiplier is scaled by the coefficient to obtain * , which fulfills also the transversality condition (97) The problem reduces to the determination of four unknown parameters ( 1 ( 0 ), 4.3. Numerical results. Each particle that forms the swarm is given a position vector and a velocity vector (for further details, cf. Appendix 2). At a given iteration the position vector of particle includes the values of the unknown parameters. For the problem at hand: Three equality constraints (related with the final conditions for orbit injection) are involved in the problem at hand. The control is expressed as a function of the costate through the necessary conditions for optimality (96), and therefore the satisfaction of the final conditions for orbit injection suffices to ensure the optimality of the solution. Hence, the following fitness function can be considered by the PSO algorithm:˜ For this example the PSO algorithm uses a population of 50 particles and is run for 500 iterations ( = 50 and = 500); is set to 1 (for = 1, 2, 3). Figure 19 illustrates the objective evolution as a function of the iteration index and points out that the PSO algorithm produces accurate results also after 400 iterations in both cases. The two transfers of interest are completed in the following minimum times:  Figure 20 portrays the optimal trajectories, and the corresponding optimal control time histories and adjoint variables for the two cases taken into account. Figure 19. Objective evolution as a function of the iteration index It is worth noticing that the constraint reduction allows arbitrarily defining the search space for the initial values of the Lagrange multipliers. This means that they can be sought in the interval −1 ≤ ( 0 ) ≤ 1 by the PSO algorithm, and only a posteriori their correct values (fulfilling also the transversality condition (97)) can be recovered, as discussed in the previous section. With reference to impulsive transfers, the analysis included in this paper leads to the determination of the number and magnitudes of the impulsive changes of velocity that perform the global optimal transfer between two arbitrary Keplerian trajectories, with possible constraints on the radius of closest approach or greatest recession from the attracting body. With regard to ellipse-to-ellipse transfers, the Hohmann transfer and the bielliptic transfer with maximum apoapse radius of the midcourse trajectory are both globally optimal. This property is proven through a direct approach, which is based only on ordinary calculus, in conjunction with a simple graphical procedure. The final selection between these two transfers depends on the apoapse and periapse radii of the terminal orbits. A direct, analytic procedure capable of taking into account all of the ellipse-to-ellipse transfers is proposed. However, in general, the Hohmann transfer turns out to be unequivocally globally optimal if 1/9 < / < 9 (where and denote the periapse radii of the final and the initial orbit, respectively), whereas it is outperformed by an arbitrary bielliptic transfer if / > 15.582 or / < 0.06418. Actually, analytical expressions exist for all the limiting values that determine unequivocally the global optimality of either the Hohmann transfer or the bielliptic transfer. Globally optimal transfers involving escape trajectories have also been investigated. The globally optimal transfer between two hyperbolic trajectories is performed through a single (finite) impulse, which occurs at the minimum periapse radius (i.e. at the radius of allowed closest approach to the attracting body). Lastly, the global optimal transfer from an elliptic orbit to an escape trajectory is performed either through a two-impulse sequence or through a four-impulse transfer.
Section 3 covers finite-thrust transfers. The direct collocation with nonlinear programming is successfully applied to optimizing several cases of orbit transfers between coplanar, circular orbits. Since the trajectory structure resembles a Hohamnn transfer, the results obtained for optimal impulsive transfers are profitably employed -using a simple procedure derived through analysis -to generate first attempt approximate solutions, which are needed by the direct algorithm. The numerical results unequivocally demonstrate that the impulsive approximation is highly accurate, when applied to modeling orbital transfers, for a wide range of values of the spacecraft propulsive characteristics. The numerical results also show the effectiveness and reliability of the numerical optimization method, which does not utilize the analytical necessary conditions for optimality, but instead is based on converting the optimal control problem into a nonlinear programming problem. Section 4 is concerned with the optimization of continuous, low-thrust orbital transfers. Unlike the finite-thrust transfers considered in Section 3, the use of lowthrust propulsion prevents employing the results obtained for impulsive transfers. The problem is particularly challenging, since a low control authority is distributed over a long time interval. This circumstance has the unfavorable consequence that a considerable number of parameters would be needed if the direct collocation with nonlinear programming algorithm would be used. In addition, providing a suitable guess would not be straightforward. This paper describes the use of the particle swarm optimization technique, a recently introduced heuristic method that does not need any starting guess. This technique is successfully applied to optimizing low-thrust orbit transfer. In this context, the necessary conditions for optimality arising from the calculus of variations come back into play. In fact, they are employed to express the time-varying control variable as a function of the costate, with the purpose of defining a low-dimensional parameter set. The numerical results, attained for two test cases, prove that the swarming methodology is well suited for treating continuous, low-thrust orbit transfer optimization problems. In fact, the method appears reliable and accurate, despite its intuitiveness and simplicity. However, as alternative stochastic techniques, the particle swarm algorithm can occasionally encounter difficulties when dealing with equality constraints. This study employs a simple approach for treating equality constraints, and show the way of reducing their number, with regard to the problem of interest. In the presence of unsatisfactory results several options are available for improving the performance attainable by the swarming algorithm. First, the number of particles and iterations can be increased to achieve more satisfactory results. Secondly, the fine tuning of the weighting coefficients in the formula for velocity update could yield improved results, at least with regard to specific cases. In addition, two versions of the particle swarm algorithm exist (i.e. the local version and the global version), with their respective advantages and drawbacks. The local version can be employed if the global version is suspected to have yielded a local minimum. However, swarming theory is a fast developing discipline and one can conclude that the choice of a proper implementation of PSO algorithm is (still) problem dependent.  Fig. 21(a)). The method of solution at hand requires the discretization of the continuous variables, i.e. the state and the control components, which are represented by a set of parameters. Their number corresponds to the degree of the polynomial that approximates each state variable in each subarc. In this research, the fifth-degree Gauss-Lobatto quadrature rule [35] is employed, because it is characterized by an improved accuracy when compared to lower degree rules, such as Simpson's rule. Hence, in each subarc six coefficients ( , , , , , ) for each state component are required: ( ) = + + 2 + 3 + 4 + 5 ⇒˙ = + 2 + 3 2 + 4 3 + 5 4 The state components and their time derivatives must be continuous at ( = 1, . . . , ), and therefore: where denotes the -th component of the state equations right-hand-side. For a given quadrature rule, in the subarc [ −1 , ] the collocation points are special points that ensure the maximum numerical accuracy [38]. This property holds iḟ In summary, implicit integration by means of the Gauss-Lobatto quadrature rules requires satisfying the relations (105) and (107) If the state has components and the control has components, respectively (2 + ) and (4 + ) discrete values for the state and the control vector are assumed as unknown parameters for the DCNLP algorithm.
Once the optimal values of the unknown parameters have been determined, they are employed to evaluate { −1 , , }. Finally, each component of the state vector is interpolated by employing a fifth-degree polynomial function of time (cf. Eq. (104)), with coefficients obtained by solving the linear system (110)-(115).