Low-Thrust Orbital Transfers in the Two-Body Problem

Low-thrust transfers between given orbits within the two-body problem are considered; the thrust is assumed power limited. A simple method for obtaining the transfer trajectories based on the linearization of the motion near reference orbits is suggested. Required calculation accuracy can be reached by means of use of a proper number of the reference orbits. The method may be used in the case of a large number of the orbits around the attracting center; no averaging is necessary in this case. The suggested method also is applicable to the cases of partly given final orbit and if there are constraints on the thrust direction. Themethod gives an optimal solution to the linearized problem which is not optimal for the original nonlinear problem; the difference between the optimal solutions to the original and linearized problems is estimated using a numerical example. Also examples illustrating the method capacities are given.


Introduction
In this paper spiral transfers with power-limited low thrust between two given orbits are considered.The transfer trajectory may include a large number of orbits around the attracting center if the transfer is performed in a strong gravity field; this situation takes place, for example, near Earth.This makes optimization of the transfer more difficult.
There are various mathematical methods for optimization of multirevolution orbital transfers 1-7 , most of them are based on averaging of motion.However, the methods have some defects: they are complicated 3, 5 or have a limited application, that is, are applicable only to the circular or neighboring orbits 1, 6, 7 or only to the planar or coaxial orbits 2, 4, 6, 7 .This paper suggests a simple mathematical method for calculation of the orbital transfers within the two-body problem.A paper describing a generalization of the method Mathematical Problems in Engineering to any force field will be published in the Cosmic Research journal.The method is based on a linearization of motion near some reference orbits; in this respect the method is similar to the Method of Transporting Trajectory MTT for solving the two-point boundary value problem TPBVP suggested in 8 and developed in 9-12 .The method suggested here makes it possible obtaining transfer trajectories of the following types: i transfer between a given position in the initial orbit and an optimally determined position in the final orbit this case takes place, for example, if the spacecraft launches from a given orbit at a given time and phasing of the final orbit is not necessary ; ii transfer between the initial orbit and a given position in the final orbit with obtaining an optimal position of the launch this situation takes place, for instance, in the transfer to a given point of geostationary orbit if the time of launch from the initial orbit may be selected arbitrarily ; iii transfer between two given orbits with obtaining optimal launch and arrival positions this is a classical case of the orbital transfer .
The method suggested in this paper does not need any averaging of motion.Required calculation accuracy is reached by means of use of a proper number of the reference orbits: the bigger is the number, the shorter are intervals of linearization and the higher is accuracy.The method does not put any limits on the number of orbits around the attracting center and on the number of the reference orbits.The suggested method is applicable also to the case of a partly given final orbit e.g., only semimajor axis and eccentricity or only orbital energy may be given and to the case of constraints imposed on the thrust direction.
However, the method gives an optimal solution to the linearized problem and this solution is not optimal for the initial nonlinear problem.Note that this nonoptimality is not caused by the linearization errors, but is a principal sequence of the substitution of the original problem for the linearized one.A comparison of the optimal solutions to the initial and linearized problems made for a numerical example is given in the paper.Other numerical examples illustrate the method capacities.

Statement of the Problem
The spacecraft motion is described by the following: is the acceleration caused by the external forces, g {0, α}.

2.2
In the considered case of the two-body problem, Let us assume for simplicity that the electric power of the propulsion system is constant.This case takes place, for example, if a nuclear power source is used or if the transfer is performed near a planet with nearly circular orbit.Then the performance index for the power-limited thrust is where α |α|.Minimal value of the functional 2.4 gives minimal propellant consumption 13 .
Let q i , q f be 5-dimensional vectors of orbital elements defining the initial and final orbits.The problem is to find transfer trajectory between the initial and final orbits in given time T in which minimal value of the performance index 2.4 is reached.

Transfer between Given State and Given Orbit
The boundary values for the problem formulated in Section 2 may be written as where y i , y f are state vectors of the initial and final orbits.The case will be considered in this section when vector y i0 is given and vector y fT is not given.

Transfer between Neighboring Orbits
First let us assume that the initial and final orbits are close to each other.Then the equation of motion 2.1 can be linearized near the initial orbit as follows Figure 1 : The closeness of the initial and final orbits makes it possible to use the following approximate linear relation: where is a 5 × 3-dimensional matrix; as is seen from 3.18 and 3.27 , matrix Q is a part of matrix U. Using 3.20 , 3.22 , and 3.24 , vector β can be found as follows:

3.28
Now due to 3.3 , 3.11 , 3.15 , and 3.28 , the optimal thrust vector and the spacecraft state vector are

3.30
Equations 3.1 and 3.30 give the state vector of entry into the final orbit

3.31
Final orbit

Transfer trajectory
Reference orbits Initial orbit Substituting 3.29 into 2.4 and taking into account definitions 3.16 , 3.26 , the minimal value of the performance index can be found as follows:

3.32
Equations 3.29 , 3.30 , and 3.32 give the solution to the problem.Note that matrices U, W, Φ, Ψ are calculated in 9, 15, 16 in an explicit form which makes the suggested method analytical.

Transfer between Arbitrary Orbits
Let us consider transfer between arbitrary orbits and divide the time interval T into n subintervals defined by instants t 0 t 1 , . . ., t n−1 , t n T .Also we assume that n−1 intermediate reference orbits between the initial and final orbits are specified somehow and q 1 , . . ., q n−1 are 5-dimentional vectors of elements of the reference orbits see Figure 2 .These vectors may be given, for instance, in the following way: q j q j j n q f − q i , j 1, . . ., n − 1 .

3.33
Let us define vectors Δq j q j − q j−1 , j 1, . . ., n , 3.34 where q 0 q i , q n q f here subscript "0" is the number of the reference orbit .Equations 3.21 and 3.34 give the following equality: Δq j Δq.

3.35
Let us divide the transfer trajectory into n arcs corresponding to the time subintervals and assume that the jth arc begins in the j-1st reference orbit and ends in the jth one.Assuming Mathematical Problems in Engineering number n big enough to make the j-1st and jth reference orbits close to each other for all j 1, . . ., n, the results of the Section 3.1 may be applied to each pair of the reference orbits.Now the problem is to find the reference orbits providing optimality of the transfer.Due to 3.32 the performance index for the jth transfer arc is where similarly to 3.26 , 3.27 the following definitions are used: where 41 v j is velocity vector of the jth reference orbit.Performance index of the whole problem is

3.38
In order to find the transfer trajectory that gives minimum value for J, it is sufficient to find intermediate reference orbits providing a minimum for 3.38 .Thus, function 3.38 should be minimized with respect to vectors Δq j , j 1, . . ., n taking into account 3.35 .Let us introduce the helping function where λ is a Lagrange multiplier.Necessary conditions of a minimum for the functional 3.38 are W j .

3.43
Now new values of vectors Δq j can be found from 3.41 and 3.42 as follows: Δq j W j W −1 Δq j 1, . . ., n .

3.44
New reference orbits are defined by the new values of the orbital elements given by q j 1 q j Δq j j 0, . . ., n − 1 .

3.45
Assuming optimal locations of the reference orbits found and applying 3.29 and 3.30 to the jth arc of the transfer trajectory, we obtain optimal thrust vector and the trajectory state vector in the time interval t j−1 ≤ t ≤ t j j 1, . . ., n as follows:

3.49
Due to 3.42 and 3.49 minimal value of the performance index in the considered case is given by 3.32 with matrix W defined by 3.43 .

Calculation Procedure
Let y 0 j , y 1 j be the state vectors of the jth reference orbit at the beginning and at the end of the jth time subinterval i.e., at times t j−1 , t j resp. .Then, the solution to the problem considered here may be obtained by means of the following iterative calculation procedure.
1 n−1 intermediate reference orbits are specified somehow, for example, using 3.33 .
The launch position in the initial orbit is specified in the case considered here i.e., state y 0 0 y i 0 is given and the respective initial state vector of the transfer trajectory is x 0 y 0 0 . 2 Vector y 0 j is calculated for j 1 using the following equation, that is, similar to 3.30 : where S 1 j S j t j , U j0 U j 0 matrix Ψ jv Ψ v t j is calculated in the j − 1st reference orbit.3 Step 2 is repeated for j 2, . . ., n − 1 and for j n 3.50 gives state vector y fT in the final orbit.
4 New vectors Δq j are calculated using 3.44 and new reference orbits with elements given by 3.45 are found.
5 Performance index is calculated using 3.32 and steps 2-4 are repeated until decrement ΔJ of the performance index gets smaller than a given parameter ε > 0.
As soon as |ΔJ| < ε, the thrust vector α and the state vector x of the transfer trajectory may be calculated at each time subinterval using 3.46 and 3.47 .
The suggested method is approximate, although any desired accuracy can be reached by means of selecting an appropriate amount n of subintervals.It can be shown that if n −→ ∞ then the solution converges to a limit which is an accurate solution to the linearized problem.

Transfer between Given Orbit and Given State Vector
Now let us consider the case when the launch position can be selected in the initial orbit arbitrarily and the position of the entry into the final orbit is given i.e., vector y i0 is not given and vector y fT is given in 3.1 .In this case the method described in Section 3 should be applied in the backward direction with decreasing time, that is, vector y 0 n−1 of the start from n − 1st reference orbit can be found for a given state vector y fT of the arrival to the final orbit, and so forth, until vector y i0 is found.The equations necessary to solve the problem considered here can be easily derived from the equations given in Section 3.

Transfer between Two Given Orbits
Let us assume that neither launch position in the initial orbit nor arrival position in the final orbit are given and should be determined in an optimal way during solving the transfer problem.This case may be solved using the suggested method in the following way.
A first guess for the launch position should be given somehow.This position defines the vector y i0 and in the first iteration of the calculation procedure described in Section 3.3 the final state vector y fT may be found.In the second iteration of the calculation procedure for this state vector a new value of the vector y i0 may be found as described in Section 4.1, and so forth, that is, odd iterations of the calculation procedure use the case described in Section 3 and even iterations use the case described in Section 4.1.

Partly Given Final Orbit
The suggested method can be also used in the case of partly given elements of the final orbit, that is, if vector q f has dimension m < 5.For instance, only energy of the final orbit m 1 or semimajor axis and eccentricity m 2 may be given.In this case the method described above is applied with m-dimensional vectors q i , q f , q j j 1, . . ., n − 1 , m × 6-dimensional matrix U j , and m-order matrices W j , S j .Nongiven orbital elements are determined using transversality condition 3.17 and the respective conditions for vectors q j .

Constraint Imposed on the Thrust Direction
Here the case when a constraint is imposed on the thrust direction is considered.This constraint may be caused by specific features of the spacecraft design or of its attitude control system.Let us assume the constraint given by Bα 0, 6.1 where B B x, t is a matrix of dimension 1 × 3 i.e., B is a row or 2 × 3 in this case the thrust direction is given and the problem is to find optimal thrust value .In this case, as is shown in 17, 18 , the suggested optimization method is also applicable with matrices W j , S j and vector α j from 41 , 3.48 , and 3.46 replaced by where third-order matrix is a projector onto the constraining set given by 6.1 .Note that rank of matrix P is less than 3.This is why matrices W j given by 6.2 may be singular, while in the case of no constraint on the thrust direction matrices W j given by 41 are nonsingular in any interval of integration 9 .As is shown in 18 nonsingularity of all matrices W j given by 6.2 is a sufficient condition of feasibility of the transfer with constraint 6.1 .

On Optimality of the Method
One of the necessary conditions of optimality of the thrust vector is constancy of the Hamiltonian in the whole time interval of the transfer 14 .It can be shown that the Hamiltonian of the linearized problem given by 3.5 is constant in the solution given by the suggested method, although the Hamiltonian of the original nonlinear problem calculated in the solution of the linear problem is not constant.Thus, the method described in the paper gives the solution which is not optimal in the original problem.
A comparison of the solutions to the problem of the orbital transfer in the original formulation and linearized one was performed.Transfer between coplanar circular orbits of radii a i 1 and a f 4 with gravitational parameter of the primary body μ 1 was considered.Solution to the original problem was provided by Petukhov.He used the mathematical developed by his method for solving two-point boundary value problem 19 ; optimal solution to the orbital transfer problem was obtained by Petukhov by means of variation of the arrival point in the final orbit.Results of the comparison are given in Table 1.
As is seen in Table 1, the difference between the values of the performance index for the optimal solutions to the original and linearized problems is small and practically does not depend on the transfer duration.Although in more complicated transfers the difference may be bigger.

Illustrative Examples
This section demonstrates potentialities of the suggested method by means of examples of transfers in the Earth's sphere of influence.Orbits are given by the orbital elements q {r π , r α , i, Ω, ω}, where r π , r α are radii of perigee and apogee in thousands of kilometers Mm , i is the orbital inclination, Ω is the longitude of the ascending node, ω is the argument of the perigee.Angular elements are given in degrees.Direction of the thrust vector is given in the examples by angles between the projection of the thrust vector onto the instantaneous orbital plane angle ϕ and between the thrust vector and orbital plane angle Ψ .In all examples given below the optimal positions of departure from the initial orbit and arrival to the final orbit were obtained i.e., the classical case of orbital transfer described in Section 4.2 was considered .The changes ΔV of the velocity by means of the thrust also are given in the examples.These changes were calculated using numerical integration of the magnitude of the thrust vector given by 3.46 or 6.2 .The number of the time subintervals was taken n 5000 in all examples.

Transfer between Elliptical Orbits with High Mutual Inclination
Transfer between two orbits given by q i {7, 30, 50, 80, −60}, q f {40, 80, 80, −80, 70} with T 400 hour is considered here.The transfer trajectory is shown in Figure 3 in projections onto the equator plane xy and the polar plane xz.Dashed lines show node lines of the initial and final orbits.The jet acceleration value divided by g 9.8066 m/s 2 is shown in Figure 4; Figure 5 shows angles ϕ, Ψ. Performance index and total ΔV for the transfer are J 44.42 m 2 s −3 , ΔV 10.05 km/s.

Transfer to an Orbit with Given Perigee and Apogee Radii
Transfer in time T 400 hrs to a partly given orbit, namely, to an orbit with given only perigee and apogee radii, is considered here.Optimal transfer is planar in this case.Only perigee and apogee radii of the initial orbit are specified, because the other initial orbital elements are not important and may be taken equal to zero.Elements of the initial and final orbits are taken as follows: q i {7, 20, 0, 0, 0}, q f {40, 80, 0, 0, 0}.Transfer orbit is shown in Figure 6; as is seen in this figure optimal is coaxial collocation of the final orbit with respect to the initial orbit.Figure 7 shows respective propulsion acceleration value divided by g, and angle ϕ is  shown in Figure 8 Ψ 0 in the planar transfer .Performance index and total ΔV are J 3.01 m 2 s −3 , ΔV 2.82 km/s.

Transfer to an Orbit with Given Energy
A transfer to an orbit with a given energy, namely, to a hyperbolic orbit given only by C 3 1 km 2 /s 2 , is considered here.Initial orbit with elements q i {7, 20, 0, 0, 0} is taken, the transfer duration is T 1000 hrs.The transfer trajectory which is obviously planar in this case and the jet acceleration are shown in Figures 7 and 8.The dot at the end of the spiral trajectory in Figure 7 marks the entry into the hyperbolic orbit.Corresponding performance index and total ΔV are J 4.85 m 2 s −3 , ΔV 4.93 km/s.The thrust direction is tangential in this case; this confirms nonoptimality of the solution to the linearized problem of the orbital transfer, because in the original nonlinear problem the optimal thrust is not tangential see, e.g., 20 .

Constrained Thrust Direction
A transfer between the orbits given by the elements q i {7, 20, 0, 0, 0}, q f {50, 100, 60, 60, 60}, in a time T 400 hours, is considered here.The transfer trajectory is shown in Figure 9.The performance index and total ΔV are J 10.83 m2 s −3 , ΔV 5.39 km/s.Now let us consider the following constraint on the thrust direction: the thrust is always orthogonal to the spacecraft position vector, that is, B r t in 6.1 .The projector 6.3 in this case is P I−rr t /r 2 .The transfer trajectory for the constrained thrust direction visually does not differ from the one for the unconstrained direction shown in Figure 9. Performance index and total ΔV in the case of the constrained thrust direction are J 11.44 m 2 s −3 , ΔV 5.52 km/s.
Acceleration value versus time for the unconstrained and the constrained thrust direction is shown in Figure 10, the transfer trajectories are shown in Figure 11, the Jet acceleration for the transfers are shown in Figure 12 and the angles of the thrust direction for the transfers are shown in Figure 13.

Conclusion
The mathematical method for calculation of low-thrust orbital transfers presented in this paper has two essential disadvantages, as follows.
1 The method is applicable to the power-limited thrust, while the existing thrusters are close to ones with constant or given exhaust velocity.These disadvantages are compensated by simplicity of the method and its analytical form at each iteration, and also by the wide applicability of the method: it works well in the case of a big difference between the initial and final orbits, for a very high number of orbits around the attracting center; also it is applicable for different transfer types such as point-to-orbit, orbitto-point, and orbit-to-orbit transfers , in the cases of partly given final orbit and of a constraint imposed on the thrust direction.Despite the fact that the method is based on linearization of motion, any necessary accuracy of calculations may be reached by means of augmentation of the number n of the reference orbits.
The suggested method may be used at early phases of the mission design when a high optimization accuracy is not needed and at the same time massive calculations are necessary for selection of a best mission scheme.

Nomenclature r, v:
Position and velocity vectors t: Current time t 0, t T : Initial and final instants of the transfer I: Unit matrix α: Jet acceleration vector thrust vector α: |α| μ: Gravitational parameter of the attracting center ϕ: Angle between the projection of the thrust vector onto the orbital plane and the velocity vector ψ: Angle between the thrust vector and the orbital plane Subscripts "0" and "T ": Values of the parameters at the time instants 0 and T if another meaning of the "0" subscript is not stipulated , superscript "t" denotes transposition.

Figure 3 :
Figure 3: Transfer between two elliptic orbits with high mutual inclination.

Figure 4 :Figure 5 :
Figure 4: Jet acceleration for transfer between two elliptic orbits with a high mutual inclination.

Figure 6 :Figure 7 :
Figure 6: Transfer trajectory to an orbit with given perigee and apogee radii.

10 Figure 8 :
Figure 8: Angle ϕ for transfer to orbit with given perigee and apogee distances.

Figure 9 :
Figure 9: Transfer trajectory to the hyperbolic orbit with given energy.

Figure 10 :Figure 11 :
Figure 10: Jet acceleration for transfer to hyperbolic orbit with given energy.

Figure 12 :
Figure 12: Jet acceleration for transfer without a and with b constraint on the thrust direction.

Figure 13 :
Figure 13: Angles of the thrust direction for transfer without a and with b constraint on the thrust direction.

Table 1 :
Comparison of the solutions to the original and linearized problems.