# Low-Thrust Orbital Transfers in the Two-Body Problem

Authors: A. A. Sukhanov; A. F. B. A. Prado
Journal: Mathematical Problems in Engineering
 (2012)

---

## Abstract

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. The method 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.


---

## Body

## 1. 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 to any force field will be published in theCosmic 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.


## 2. Statement of the Problem

The spacecraft motion is described by the following:(2.1)ẋ=f(x,t)+g,
where x={r,v}  is the spacecraft state vector,  f=f(x,t)={v,fv}, fv=fv(x,t) is the acceleration caused by the external forces,(2.2)g={0,α}.
In the considered case of the two-body problem,(2.3)fv=-μr3r.
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(2.4)J=12∫0Tα2dt,
where α=|α|. Minimal value of the functional (2.4) gives minimal propellant consumption [13].Letqi,qf 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.


## 3. Transfer between Given State and Given Orbit

The boundary values for the problem formulated in Section2 may be written as(3.1)x0=yi0,xT=yfT,
where yi,yf are state vectors of the initial and final orbits. The case will be considered in this section when vector yi0 is given and vector yfT is not given.


### 3.1. 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):(3.2)ξ̇=Fξ+g,
where(3.3)ξ=ξ(t)=x(t)-yi(t),
is state vector of the linearized motion and(3.4)F=∂f∂yi.
Vector f and matrix F in (3.2) and (3.4) are calculated in the initial orbit, this is why matrix F is a function of time even if vector f does not depend on time explicitly. The Hamiltonian for (3.2) and performance index (2.4) is [14](3.5)H=-α22+ptFξ+Pvtα+pt,
where p={pr,pv} is a vector of the costate variables corresponding to (3.2), pv is Lawden’s primer vector, pt is costate variable corresponding to the additional equation ṫ=1 making the system autonomous. Vector p satisfies the following costate variational equation:(3.6)ṗt=-(∂H∂ξ)t=-ptF,
variable pt satisfies the following equation:(3.7)pṫ=-∂H∂t=-pòḞξ.
Let us consider sixth-order matrixΨ=Ψ(t) which is a general solution to (3.6) with initial value Ψ(0)=I. Matrix Ψ is costate transition matrix given by(3.8)Ψ=∂yi0/∂yi.
Dividing matrix Ψ into two 6×3-dimensional submatrices as follows:(3.9)Ψ=[ΨrΨv],
the costate variables may be represented as(3.10)p=Ψtβ,pv=Ψvtβ,
where β is a constant 6-dimensional vector. The Hamiltonian (3.5) reaches its maximum if(3.11)α=pv=Ψvtβ.
Solution to (3.2) is given by the Cauchy formula(3.12)ξ=ξ(t)=ξ0+∫0tΦ(t,τ)gdτ,
where(3.13)Φ(t1,t2)=∂yi(t1)∂yi(t2)
is the state transition matrix. Due to (3.1) and (3.3) ξ0=0 in (3.12). Using (2.2), (3.9), and (3.11) and relations(3.14)Φ(t,τ)=Φ(t,0)Φ(0,τ)=Φ(t,0)Φ-1(τ,0),Φ=Φ(t,0)=Ψ-1,
the solution (3.12) can be represented as(3.15)ξ=ΦSβ,
where(3.16)S=S(t)=∫0tΨvΨvtdt
is matrix of sixth order (also see [9]). As follows from (3.3), (3.11), and (3.15), in order to find optimal thrust vector and the state vector of the transfer trajectory, it is sufficient to obtain vector β. Nonspecified state vector yfT may be found using transversality condition which in the considered case is(3.17)pT=p(T)=(∂qf∂yfT)tσ,
where σ is an arbitrary constant 5-dimensional vector. Due to the closeness of the initial and final orbits, the approximate equality ∂qf/∂yf≈U is fulfilled, where(3.18)U=U(t)=∂qi∂yi(t),
is a 5×6-dimensional matrix. Then the condition (3.17) takes the form(3.19)pT=UTtσ.
Equations (3.10) and (3.19) and the last relation of (3.14) give(3.20)β=(UTΦT)tσ.
The closeness of the initial and final orbits makes it possible to use the following approximate linear relation:(3.21)Δq=qf-qi=UTξT.
Equations (3.15), (3.20), and (3.21) give(3.22)Δq=Wσ,
where(3.23)W=UTΦTST(UTΦT)t
is 5-dimensional matrix. Using (3.8), (3.9), (3.13), (3.14), (3.16), and (3.18), we obtain(3.24)UTΦT=∂qi∂yiT∂yiT∂yi0=∂qi∂yi0=U0,(3.25)UTΦTST(UTΦT)t=∂qi∂yi0∫0T∂yi0∂vi(∂yi0∂vi)tdt(∂qi∂yi0)t,
where vi=vi(t) is velocity vector in the initial orbit. (3.23), (3.25) give(3.26)W=U0STU0t=∫0TQQtdt,
where(3.27)Q=Q(t)=-∂qi∂vi
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)β=U0tW-1Δq.
Now due to (3.3), (3.11), (3.15), and (3.28), the optimal thrust vector and the spacecraft state vector are(3.29)α=ΨvtU0tW-1Δq,(3.30)x(t)=yi(t)+ΦSU0tW-1Δq.
Equations (3.1) and (3.30) give the state vector of entry into the final orbit(3.31)yfT=yiT+ΦTSTU0tΔq.
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)J=12ΔqtW-1Δq.
Equations (3.29), (3.30), and (3.32) give the solution to the problem.Figure 1
Transfer between given orbits.Note that matricesU, W,Φ, Ψ are calculated in [9, 15, 16] in an explicit form which makes the suggested method analytical.


### 3.2. Transfer between Arbitrary Orbits

Let us consider transfer between arbitrary orbits and divide the time intervalT into n subintervals defined by instants  t0=t1,…,tn-1,tn=T. Also we assume that n-1 intermediate reference orbits between the initial and final orbits are specified somehow and q1,…,qn-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:(3.33)qj=qj+jn(qf-qi),(j=1,…,n-1).
Let us define vectors(3.34)Δqj=qj-qj-1,(j=1,…,n),
where q0=qi,qn=qf (here subscript “0” is the number of the reference orbit). Equations (3.21) and (3.34) give the following equality:(3.35)∑j=1nΔqj=Δq.
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 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(3.36)Jj=12ΔqjtWj-1Δqj,(j=1,…,n),
where similarly to (3.26), (3.27) the following definitions are used:(3.37)Wj=∫tj-1tjQjQjtdt,Qj=∂qj-1∂vj-1(j=1,…,n),
where (41) vj is velocity vector of the jth reference orbit. Performance index of the whole problem is(3.38)J=∑j=1nJj.
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 Δqj,j=1,…,n taking into account (3.35). Let us introduce the helping function(3.39)L=J-λt(∑j=1nΔqj-Δq),Figure 2
Transfer between arbitrary orbits.whereλ is a Lagrange multiplier. Necessary conditions of a minimum for the functional (3.38) are(3.40)(∂L∂Δqj)t=Wj-1Δqj-λ=0(j=1,…,n).
Thus,(3.41)Δqj=Wjλ(j=1,…,n),
and (3.35) and (3.41) give(3.42)λ=W-1Δq,
where(3.43)W=∑j=1nWj.
Now new values of vectors Δqj can be found from (3.41) and (3.42) as follows:(3.44)Δqj=WjW-1Δq(j=1,…,n).
New reference orbits are defined by the new values of the orbital elements given by(3.45)qj+1=qj+Δqj(j=0,…,n-1).
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 tj-1≤t≤tj(j=1,…,n) as follows:(3.46)α(t)=ΨjvtUjtW-1Δq,(3.47)x(t)=yj-1(t)+ΦjSjUj0tW-1Δq,
where(3.48)Sj=Sj(t)=∫tj-1tΨjvΨjvtdt,Uj=∂qj-1∂yj-1,
matrix W is given by (3.43), matrices Φ=Φj(t,0),Ψjv=Ψjv(t) are calculated in the j-1st reference orbit.Multiplying (3.40) by Δqj and taking into account (3.35), (3.36), and (3.38), we obtain(3.49)J=12λtΔq.
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).


### 3.3. Calculation Procedure

Letyj0,yj1 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 tj-1, tj 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 y00=yi(0) is given) and the respective initial state vector of the transfer trajectory is x0=y00.(2)
Vectoryj0 is calculated for j=1 using the following equation, that is, similar to (3.30):(3.50)yj0=yj-11+Φj1Sj1Uj0tWj-1Δqj,
where Sj1=Sj(tj),Uj0=Uj(0)  matrix Ψjv=Ψv(tj) is calculated in the j-1st reference orbit.(3)
Step  2 is repeated forj=2,…,n-1 and for j=n (3.50) gives state vector yfT in the final orbit.(4)
New vectorsΔqj 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 amountn 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.


## 3.1. 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):(3.2)ξ̇=Fξ+g,
where(3.3)ξ=ξ(t)=x(t)-yi(t),
is state vector of the linearized motion and(3.4)F=∂f∂yi.
Vector f and matrix F in (3.2) and (3.4) are calculated in the initial orbit, this is why matrix F is a function of time even if vector f does not depend on time explicitly. The Hamiltonian for (3.2) and performance index (2.4) is [14](3.5)H=-α22+ptFξ+Pvtα+pt,
where p={pr,pv} is a vector of the costate variables corresponding to (3.2), pv is Lawden’s primer vector, pt is costate variable corresponding to the additional equation ṫ=1 making the system autonomous. Vector p satisfies the following costate variational equation:(3.6)ṗt=-(∂H∂ξ)t=-ptF,
variable pt satisfies the following equation:(3.7)pṫ=-∂H∂t=-pòḞξ.
Let us consider sixth-order matrixΨ=Ψ(t) which is a general solution to (3.6) with initial value Ψ(0)=I. Matrix Ψ is costate transition matrix given by(3.8)Ψ=∂yi0/∂yi.
Dividing matrix Ψ into two 6×3-dimensional submatrices as follows:(3.9)Ψ=[ΨrΨv],
the costate variables may be represented as(3.10)p=Ψtβ,pv=Ψvtβ,
where β is a constant 6-dimensional vector. The Hamiltonian (3.5) reaches its maximum if(3.11)α=pv=Ψvtβ.
Solution to (3.2) is given by the Cauchy formula(3.12)ξ=ξ(t)=ξ0+∫0tΦ(t,τ)gdτ,
where(3.13)Φ(t1,t2)=∂yi(t1)∂yi(t2)
is the state transition matrix. Due to (3.1) and (3.3) ξ0=0 in (3.12). Using (2.2), (3.9), and (3.11) and relations(3.14)Φ(t,τ)=Φ(t,0)Φ(0,τ)=Φ(t,0)Φ-1(τ,0),Φ=Φ(t,0)=Ψ-1,
the solution (3.12) can be represented as(3.15)ξ=ΦSβ,
where(3.16)S=S(t)=∫0tΨvΨvtdt
is matrix of sixth order (also see [9]). As follows from (3.3), (3.11), and (3.15), in order to find optimal thrust vector and the state vector of the transfer trajectory, it is sufficient to obtain vector β. Nonspecified state vector yfT may be found using transversality condition which in the considered case is(3.17)pT=p(T)=(∂qf∂yfT)tσ,
where σ is an arbitrary constant 5-dimensional vector. Due to the closeness of the initial and final orbits, the approximate equality ∂qf/∂yf≈U is fulfilled, where(3.18)U=U(t)=∂qi∂yi(t),
is a 5×6-dimensional matrix. Then the condition (3.17) takes the form(3.19)pT=UTtσ.
Equations (3.10) and (3.19) and the last relation of (3.14) give(3.20)β=(UTΦT)tσ.
The closeness of the initial and final orbits makes it possible to use the following approximate linear relation:(3.21)Δq=qf-qi=UTξT.
Equations (3.15), (3.20), and (3.21) give(3.22)Δq=Wσ,
where(3.23)W=UTΦTST(UTΦT)t
is 5-dimensional matrix. Using (3.8), (3.9), (3.13), (3.14), (3.16), and (3.18), we obtain(3.24)UTΦT=∂qi∂yiT∂yiT∂yi0=∂qi∂yi0=U0,(3.25)UTΦTST(UTΦT)t=∂qi∂yi0∫0T∂yi0∂vi(∂yi0∂vi)tdt(∂qi∂yi0)t,
where vi=vi(t) is velocity vector in the initial orbit. (3.23), (3.25) give(3.26)W=U0STU0t=∫0TQQtdt,
where(3.27)Q=Q(t)=-∂qi∂vi
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)β=U0tW-1Δq.
Now due to (3.3), (3.11), (3.15), and (3.28), the optimal thrust vector and the spacecraft state vector are(3.29)α=ΨvtU0tW-1Δq,(3.30)x(t)=yi(t)+ΦSU0tW-1Δq.
Equations (3.1) and (3.30) give the state vector of entry into the final orbit(3.31)yfT=yiT+ΦTSTU0tΔq.
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)J=12ΔqtW-1Δq.
Equations (3.29), (3.30), and (3.32) give the solution to the problem.Figure 1
Transfer between given orbits.Note that matricesU, W,Φ, Ψ are calculated in [9, 15, 16] in an explicit form which makes the suggested method analytical.


## 3.2. Transfer between Arbitrary Orbits

Let us consider transfer between arbitrary orbits and divide the time intervalT into n subintervals defined by instants  t0=t1,…,tn-1,tn=T. Also we assume that n-1 intermediate reference orbits between the initial and final orbits are specified somehow and q1,…,qn-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:(3.33)qj=qj+jn(qf-qi),(j=1,…,n-1).
Let us define vectors(3.34)Δqj=qj-qj-1,(j=1,…,n),
where q0=qi,qn=qf (here subscript “0” is the number of the reference orbit). Equations (3.21) and (3.34) give the following equality:(3.35)∑j=1nΔqj=Δq.
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 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(3.36)Jj=12ΔqjtWj-1Δqj,(j=1,…,n),
where similarly to (3.26), (3.27) the following definitions are used:(3.37)Wj=∫tj-1tjQjQjtdt,Qj=∂qj-1∂vj-1(j=1,…,n),
where (41) vj is velocity vector of the jth reference orbit. Performance index of the whole problem is(3.38)J=∑j=1nJj.
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 Δqj,j=1,…,n taking into account (3.35). Let us introduce the helping function(3.39)L=J-λt(∑j=1nΔqj-Δq),Figure 2
Transfer between arbitrary orbits.whereλ is a Lagrange multiplier. Necessary conditions of a minimum for the functional (3.38) are(3.40)(∂L∂Δqj)t=Wj-1Δqj-λ=0(j=1,…,n).
Thus,(3.41)Δqj=Wjλ(j=1,…,n),
and (3.35) and (3.41) give(3.42)λ=W-1Δq,
where(3.43)W=∑j=1nWj.
Now new values of vectors Δqj can be found from (3.41) and (3.42) as follows:(3.44)Δqj=WjW-1Δq(j=1,…,n).
New reference orbits are defined by the new values of the orbital elements given by(3.45)qj+1=qj+Δqj(j=0,…,n-1).
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 tj-1≤t≤tj(j=1,…,n) as follows:(3.46)α(t)=ΨjvtUjtW-1Δq,(3.47)x(t)=yj-1(t)+ΦjSjUj0tW-1Δq,
where(3.48)Sj=Sj(t)=∫tj-1tΨjvΨjvtdt,Uj=∂qj-1∂yj-1,
matrix W is given by (3.43), matrices Φ=Φj(t,0),Ψjv=Ψjv(t) are calculated in the j-1st reference orbit.Multiplying (3.40) by Δqj and taking into account (3.35), (3.36), and (3.38), we obtain(3.49)J=12λtΔq.
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).


## 3.3. Calculation Procedure

Letyj0,yj1 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 tj-1, tj 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 y00=yi(0) is given) and the respective initial state vector of the transfer trajectory is x0=y00.(2)
Vectoryj0 is calculated for j=1 using the following equation, that is, similar to (3.30):(3.50)yj0=yj-11+Φj1Sj1Uj0tWj-1Δqj,
where Sj1=Sj(tj),Uj0=Uj(0)  matrix Ψjv=Ψv(tj) is calculated in the j-1st reference orbit.(3)
Step  2 is repeated forj=2,…,n-1 and for j=n (3.50) gives state vector yfT in the final orbit.(4)
New vectorsΔqj 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 amountn 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.


## 4. Other Types of the Transfer

### 4.1. 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., vectoryi0 is not given and vector yfT 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 yn-10 of the start from n-1st reference orbit can be found for a given state vector yfT of the arrival to the final orbit, and so forth, until vector yi0 is found. The equations necessary to solve the problem considered here can be easily derived from the equations given in Section 3.


### 4.2. 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 vectoryi0 and in the first iteration of the calculation procedure described in Section 3.3 the final state vector yfT may be found. In the second iteration of the calculation procedure for this state vector a new value of the vector yi0 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.


## 4.1. 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., vectoryi0 is not given and vector yfT 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 yn-10 of the start from n-1st reference orbit can be found for a given state vector yfT of the arrival to the final orbit, and so forth, until vector yi0 is found. The equations necessary to solve the problem considered here can be easily derived from the equations given in Section 3.


## 4.2. 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 vectoryi0 and in the first iteration of the calculation procedure described in Section 3.3 the final state vector yfT may be found. In the second iteration of the calculation procedure for this state vector a new value of the vector yi0 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.


## 5. 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 vectorqf 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 qi,qf,qj(j=1,…,n-1), m×6-dimensional matrix Uj, and m-order matrices Wj,Sj. Nongiven orbital elements are determined using transversality condition (3.17) and the respective conditions for vectors qj.


## 6. 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(6.1)Bα=0,
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 Wj,Sj and vector αj from (41), (3.48), and (3.46) replaced by(6.2)Wj=∫tj-1tjQjPQjtdt,Sj=Sj(t)=∫tj-1tΨjvPΨjvtdt,αj=PΨvtUjtWj-1Δqj,
where third-order matrix(6.3)P=I-Bt(BBt)-1B,
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 Wj given by (6.2) may be singular, while in the case of no constraint on the thrust direction matrices Wj  given by (41) are nonsingular in any interval of integration [9]. As is shown in [18] nonsingularity of all matrices Wj given by (6.2) is a sufficient condition of feasibility of the transfer with constraint (6.1).


## 7. 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 radiiai = 1 and af = 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.Table 1
Comparison of the solutions to the original and linearized problems.
Transfer durationin periods of the initial orbitSolution to theSolution to theDifference, %original problemlinearized problemNo. of complete orbitsJNo. of complete orbitsJ2091,012×10-381,053×10-34,140184,997×10-4165,193×10-43,9100461,991×10-4422,070×10-44.0200930,995×10-4841,034×10-43.9As is seen in Table1, 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.


## 8. 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(8.1)q={rπ,rα,i,Ω,ω},whererπ,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.


### 8.1. Transfer between Elliptical Orbits with High Mutual Inclination

Transfer between two orbits given byqi={7,30,50,80,-60},  qf={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/s2 is shown in Figure 4; Figure 5 shows angles φ, Ψ. Performance index and total ΔV for the transfer are J=44.42m2s-3, ΔV = 10.05 km/s.Figure 3
Transfer between two elliptic orbits with high mutual inclination.Figure 4
Jet acceleration for transfer between two elliptic orbits with a high mutual inclination.Figure 5
Angles of the thrust direction transfer between two elliptic orbits with a high mutual inclination.


### 8.2. Transfer to an Orbit with Given Perigee and Apogee Radii

Transfer in timeT=400hrs 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: qi={7,20,0,0,0},qf={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 m2s−3, ΔV = 2.82 km/s.Figure 6
Transfer trajectory to an orbit with given perigee and apogee radii.Figure 7
Jet acceleration for transfer to orbit with given perigee and apogee distances.Figure 8
Angleφ for transfer to orbit with given perigee and apogee distances.


### 8.3. Transfer to an Orbit with Given Energy

A transfer to an orbit with a given energy, namely, to a hyperbolic orbit given only byC3 = 1 km2/s2, is considered here. Initial orbit with elements qi={7,20,0,0,0}  is taken, the transfer duration is T=1000hrs. 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 m2 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]).


### 8.4. Constrained Thrust Direction

A transfer between the orbits given by the elementsqi={7,20,0,0,0},qf={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 m2s−3, ΔV=5.39 km/s.Figure 9
Transfer trajectory to the hyperbolic orbit with given energy.Now let us consider the following constraint on the thrust direction: the thrust is always orthogonal to the spacecraft position vector, that is,B=rt in (6.1). The projector (6.3) in this case is P=I-rrt/r2. 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 m2s −3, ΔV = 5.52 km/s.Acceleration value versus time for the unconstrained and the constrained thrust direction is shown in Figure10, 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.Figure 10
Jet acceleration for transfer to hyperbolic orbit with given energy.Transfer trajectories without (a) and with (b) constraint on the thrust direction.
(a)(b)Jet acceleration for transfer without (a) and with (b) constraint on the thrust direction.
(a)(b)Angles of the thrust direction for transfer without (a) and with (b) constraint on the thrust direction.
(a)(b)


## 8.1. Transfer between Elliptical Orbits with High Mutual Inclination

Transfer between two orbits given byqi={7,30,50,80,-60},  qf={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/s2 is shown in Figure 4; Figure 5 shows angles φ, Ψ. Performance index and total ΔV for the transfer are J=44.42m2s-3, ΔV = 10.05 km/s.Figure 3
Transfer between two elliptic orbits with high mutual inclination.Figure 4
Jet acceleration for transfer between two elliptic orbits with a high mutual inclination.Figure 5
Angles of the thrust direction transfer between two elliptic orbits with a high mutual inclination.


## 8.2. Transfer to an Orbit with Given Perigee and Apogee Radii

Transfer in timeT=400hrs 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: qi={7,20,0,0,0},qf={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 m2s−3, ΔV = 2.82 km/s.Figure 6
Transfer trajectory to an orbit with given perigee and apogee radii.Figure 7
Jet acceleration for transfer to orbit with given perigee and apogee distances.Figure 8
Angleφ for transfer to orbit with given perigee and apogee distances.


## 8.3. Transfer to an Orbit with Given Energy

A transfer to an orbit with a given energy, namely, to a hyperbolic orbit given only byC3 = 1 km2/s2, is considered here. Initial orbit with elements qi={7,20,0,0,0}  is taken, the transfer duration is T=1000hrs. 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 m2 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]).


## 8.4. Constrained Thrust Direction

A transfer between the orbits given by the elementsqi={7,20,0,0,0},qf={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 m2s−3, ΔV=5.39 km/s.Figure 9
Transfer trajectory to the hyperbolic orbit with given energy.Now let us consider the following constraint on the thrust direction: the thrust is always orthogonal to the spacecraft position vector, that is,B=rt in (6.1). The projector (6.3) in this case is P=I-rrt/r2. 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 m2s −3, ΔV = 5.52 km/s.Acceleration value versus time for the unconstrained and the constrained thrust direction is shown in Figure10, 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.Figure 10
Jet acceleration for transfer to hyperbolic orbit with given energy.Transfer trajectories without (a) and with (b) constraint on the thrust direction.
(a)(b)Jet acceleration for transfer without (a) and with (b) constraint on the thrust direction.
(a)(b)Angles of the thrust direction for transfer without (a) and with (b) constraint on the thrust direction.
(a)(b)


## 9. 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.(2)
The original nonlinear problem is replaced in the method by a linearized problem solution to which it is not optimal for the original problem.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, orbit-to-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 numbern 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.



---

*Source: 905209-2012-04-10.xml*