On the Maximal Orbit Transfer Problem

Marian Mureşan

Assume that a spacecraft is in a circular orbit and consider the problem of finding the largest possible circular orbit to which the spacecraft can be transferred with constant thrust during a set time, so that the variable parameter is the thrust-direction angle β. Also assume that there is only one center of attraction at the common center of the two circular orbits. Finally, assume normalized values for all constants and variables.

This article is divided into five sections: the orbit transfer problem, equations of motion, the optimal control problem, necessary conditions for the Mayer problem, and a dynamic approach to the maximal orbit transfer problem using Mathematica’s built-in Manipulate function.

The Earth-Mars orbit transfer problem is timely, given the successful flight and smooth landing of the American Curiosity rover on Mars.

The Orbit Transfer Problem

For the orbit transfer problem, assume that:

There is a unique center of attraction.

Initially the spacecraft moves in a circular trajectory around the center of attraction.

The spacecraft moves with a constant thrust from a rocket engine operating in the time interval [0,b].

The spacecraft moves to the largest possible circular orbit around the center of attraction.

The orbit transfer trajectory is coplanar with the two circular orbits and the center of attraction.

All these assumptions are stated in [, p. 66]. Here is a sketch of a solution to the problem with some notation. The blue curve is the orbital transfer trajectory, while the red and green curves are the initial lower circular orbit and the final upper circular orbit.

Module[
{f,P,Q},
f=BezierFunction[{{3,0},{3.3,1.6},{1,5},{-2.3,4.5},{4.5 Cos[5 Pi/6],4.5 Sin[5 Pi/6]}}];
P=f[0.5];
Q=P[[2]](P[[2]]-3.75)/P[[1]]+P[[1]];
Show[
Graphics[{
{Red,Circle[{0,0},3,{0,3Pi/4}]},
{Green,Circle[{0,0},4.5,{Pi/6,5Pi/6}]},
Arrowheads[.023],
Arrow[{P+0.4{Cos[ Pi/4-0.13],Sin[ Pi/4-0.13]},P+0.4{Cos[ Pi/4-0.135],Sin[ Pi/4-0.135]}}],
Arrow[{0.55{Cos[ Pi/3+0.31],Sin[ Pi/3+0.31]},0.55{Cos[ Pi/3+0.32],Sin[ Pi/3+0.32]}}],    
Arrow[{{0,0},{3,0}}], 
Arrow[{{0,0 },P}],  
Arrow[{{0,0},4.5{Cos[5 Pi/6],Sin[5 Pi/6]}}], 
Arrow[{{P[[2]](P[[2]]-3.538)/P[[1]]+P[[1]], 3.538},{Q,3.75}}],
Arrow[{{Q,3.75},{Q,3.75}+0.12P}],
Arrow[{P,{Q,3.75}+0.12P}],   
{Red,Arrow[{P,P+{0.8,.65}}]}, 
PointSize[0.012], Point[{3,0}],
Point[{4.5 Cos[5 Pi/6],4.5 Sin[5 Pi/6]}],
Point[P],
Circle[P,0.4,{Pi-0.15,Pi/4-.14}],
Circle[{0,0},0.55,{0,Pi/3+0.32}],
Text[Style[Row[{Style["r",Italic],"(0)"}],12],{3.45,0.0}],
Text[Style[Row[{Style["r",Italic],"(",Style["b",Italic],")"}],12],4.85{ Cos[5 Pi/6],Sin[5 Pi/6]}],
Text[Style["v",Italic,12],P+{-0.65,-.12}],
Text[Style["u",Italic,12],P+{-1.4,0.49}],
Text[Style["T",Italic,Red,12],P+{0.55,0.63}],
Text[Style["β",12],P+{-0.1,0.55}],
Text[Style["θ",12],f[0.25]/4 ],
Text[Style["center of attraction",Italic,12],{0.0,-0.3}],
Text[Style["spacecraft",Italic,12],P+{.85,0.1}],
Text[Style[Row[{Style["r",Italic],"(",Style["t",Italic],")"}],12],P/2+{0.4,0}],
Text[Style["initial\norbit",Italic,12],{2.1,1}],
Text[Style["final\norbit",Italic,12],{3.3,2.3}],
{Red,Disk[{0,0},0.06]}
}],
 ParametricPlot[f[x],{x,0,1},PlotStyle->{Blue,Thickness->0.005}],ImageSize->{500,300}
]
]

The notation from [, pp. 66–68], [], or [] is:

t is time in the given interval [0,b], which is called the horizon.

r(t) is the radial distance from the center of attraction to the spacecraft; r(t) increases as fuel is burned; r(0)=(Subscript[r, 0]) is the initial distance; r(b) is the final and maximal distance.

θ(t) is the polar angle, measured counterclockwise from the straight line connecting the center of attraction with the position of the spacecraft at t=0.

u(t) is the radial velocity component.

v(t) is the tangential velocity component.

β(t) is the thrust-direction angle; it is the control variable.

(Subscript[m, 0]) is the initial mass of the spacecraft with propellant included; (Subscript[m, 0])-m't is the time-dependent mass, which decreases due to the constant fuel consumption rate m'>0.

T is the thrust, also assumed to be constant.

μ is the gravitational constant.

Equations of Motion

The equations of motion of the spacecraft consistent with the above assumptions, according to [, p. 67] and [], are

r'(t)=u(t),

u'(t)=(((v^2)(t))/(r(t)))-(μ/((r^2)(t)))+((T sin β(t))/((Subscript[m, 0])-m't)),

v'(t)=-((u(t)v(t))/(r(t)))+((T cos β(t))/((Subscript[m, 0])-m't)),

θ'(t)=((v(t))/(r(t))).

The associated boundary conditions are

r(0)=(Subscript[r, 0]),

u(0)=0,

v(0)=(Sqrt[(μ/(Subscript[r, 0]))]) ,

θ(0)=0,

u(b)=0,

v(b)=(Sqrt[(μ/(r(b)))]).

The system of nonlinear differential equations () to () with the boundary value conditions () to (), the control function β, and the maximizing condition

max r(b)

form the optimal control problem to be solved, assuming that the state functions r, u, v, and θ and the control function β are sufficiently smooth. Conditions (), (), (), and () guarantee that the trajectory of the spacecraft is tangent to the two circular orbits.

The Optimal Control Problem

The goal is to maximize r(b), the radius of the orbit transfer at the endpoint in time, so the cost functional is determined by

Λ[r,u,v,θ,β]=r(b).

Thus the horizon is [0,b] with b>0. This is a Mayer optimal control problem (see Ch. 4 in []).

Since the differential equations () to () with conditions () to () and the cost functional () are not time dependent, the optimal control problem is equivalent to either of the following two problems:

differential equations () to () with conditions () to (), a given r(b), and b finite and arbitrary, with optimality condition to minimize b

differential equations () to () with conditions () to (), a given r(b), with the optimality condition to minimize the fuel consumption (∫_0^b)((Subscript[m, 0])-m't)dt

Theorem 1

Under the hypotheses of Filippov’s theorem (theorem 9.2.i of []), the optimal control problem () to () with conditions () to () and the maximizing functional () and () has an absolute maximum in the nonempty set Ω of admissible pairs.

Necessary Conditions for a Mayer Problem

For brevity, here is an abbreviated version of theorem 4.2.i in []: Let the Mayer problem be expressed as

Λ[x,u]=g(a,x(a),b,x(b)) (cost functional),

((dx)/(dt))=f(t,x(t),u(t)), t∈[a,b],almost everywhere (a.e.) (differential constraint),

e[x]=(a,x(a),b,x(b))∈B⊂(ℝ^(1+n+1+n)) (boundary conditions),

(t,x(t))∈A,t∈[a,b] (time and state constraint),

u(t)∈U(t),t∈[a,b] (control constraint).

A pair (x(t),u(t)), a≤t≤b, is said to be admissible (or feasible) provided that x:[a,b]->(ℝ^n) is absolutely continuous [], u:[a,b]->(ℝ^p) is measurable, and x and u satisfy () a.e. Let Ω be the class of admissible pairs (x,u). The goal is to find the minimum of the cost functional () over Ω, that is, to find an element ((SuperStar[x]),(u^⋆))∈Ω so that -∞<Λ[(x^⋆),(u^⋆)]≤Λ[x,u] for all (x,u)∈Ω. Introduce the variables λ=((Subscript[λ, 1]),…,(Subscript[λ, n])), called multipliers, and an auxiliary function H(t,x,u,λ), called the Hamiltonian, defined on T×U×(ℝ^n) by

H(t,x,u,λ)=(∑_(i=1)^n)(Subscript[λ, i])(Subscript[f, i])(t,x,u).

Define

M(t,x,λ)=(inf_(u∈U(t)))H(t,x,u,λ).

Further necessary assumptions:

There exists an element ((x^⋆),(u^⋆))∈Ω such that -∞<Λ[(x^⋆),(u^⋆)]≤Λ[x,u] for all (x,u)∈Ω.

A is closed in (ℝ^(1+n)).

The set S={(t,x,u)|(t,x)∈A,u∈U(t)} is closed in (ℝ^(1+n+p)).

f∈(C^1)(S,ℝ).

Notation:

(Subscript[f, it])=((∂(Subscript[f, i]))/(∂t)), (Subscript[f, i,(Subscript[x, j])])=((∂(Subscript[f, i]))/(∂(Subscript[x, j]))), (Subscript[H, (Subscript[x, j])])=((∂H)/(∂(Subscript[x, j])))=(∑_(i=1)^n)(Subscript[λ, i]) (Subscript[f, i,(Subscript[x, j])]), (Subscript[H, t])=((∂H)/(∂t))=(∑_(i=1)^n)(Subscript[λ, i]) (Subscript[f, i,t]), (Subscript[H, (Subscript[λ, j])])=((∂H)/(∂(Subscript[λ, j])))=(Subscript[f, j]).

The graph {(t,(x^⋆)(t))|a≤t≤b} of the optimal trajectory (x^⋆) belongs to the interior of A.

U does not depend on time and is a closed set.

The endpoint e((x^⋆))=(a,(x^⋆)(a),b,(x^⋆)(b)) of the optimal trajectory (x^⋆) is a point of B, where B has a tangent variety B' (of some dimension δ, 0≤δ≤2n+2) whose vectors are denoted by

h=(a,(Subscript[ξ, 1]),b,(Subscript[ξ, 2])),(Subscript[ξ, 1])=((Subscript[ξ, 1]^1),…,(Subscript[ξ, 1]^n)),(Subscript[ξ, 2])=((Subscript[ξ, 2]^1),…,(Subscript[ξ, 2]^n)),

or by

h=(da,d(Subscript[x, 1]),db,d(Subscript[x, 2])),d(Subscript[x, 1])=(d(Subscript[ξ, 1]^1),…,d(Subscript[ξ, 1]^n)),d(Subscript[x, 2])=(d(Subscript[ξ, 2]^1),…,d(Subscript[ξ, 2]^n)).

Theorem 2

Assume the above eight hypotheses and let ((x^⋆),(u^⋆)) be an optimal pair for the Mayer problem () and (). Then the optimal pair ((x^⋆),(u^⋆)) necessarily has the following properties:

((d(Subscript[λ, i]))/(dt))=-(Subscript[H, (Subscript[x, i])])(t,(x^⋆)(t),(SuperStar[u])(t),λ(t)),i=1,…,n,t∈[a,b] (a.e.).

((dM)/(dt))=(d/(dt))M(t,(x^⋆)(t),λ(t))=(Subscript[H, t])(t,(x^⋆)(t),(u^⋆)(t),λ(t)),t∈[a,b] (a.e.).

((Subscript[λ, 0])(Subscript[g, a])-M(a))da+(∑_(i=1)^n)((Subscript[λ, 0])(Subscript[g, (Subscript[x, 1]^i)])+(Subscript[λ, i])(a))d(Subscript[x, 1]^i)+((Subscript[λ, 0])(Subscript[g, b])+M(b))db+(∑_(i=1)^n)((Subscript[λ, 0])(Subscript[g, (Subscript[x, 2]^i)])-(Subscript[λ, i])(b))d(Subscript[x, i]^2)=0,

for every vector h=(da,d(Subscript[x, 1]),db,d(Subscript[x, 2]))∈B'.

From () and (a) of theorem 2, the Hamiltonian and the equations for the multipliers for () to () are

H=(Subscript[λ, 1])u(t)+(Subscript[λ, 2])((((v^2)(t))/(r(t)))-(μ/((r^2)(t)))+((T sin β(t))/((Subscript[m, 0])-m't)))+(Subscript[λ, 3])(-((u(t)v(t))/(r(t)))+((T cos β(t))/((Subscript[m, 0])-m't)) )+(Subscript[λ, 4])((v(t))/(r(t))) ,

(Subscript[λ, 1])'(t)=-((∂H)/(∂r))= -(Subscript[λ, 2])(t)(-(((v^2)(t))/((r^2)(t)))+((2μ)/((r^3)(t))))-(Subscript[λ, 3])(t)(((u(t)v(t))/((r^2)(t))) )-(Subscript[λ, 4])(t)(-((v(t))/((r^2)(t))) ),

(Subscript[λ, 2])'(t)=-((∂H)/(∂u))=-(Subscript[λ, 1])(t)+(Subscript[λ, 3])(t)((v(t))/(r(t))),

(Subscript[λ, 3])'(t)=-((∂H)/(∂v))=-(Subscript[λ, 2])(t)((2v(t))/(r(t)))+(Subscript[λ, 3])(t)((u(t))/(r(t)))-(Subscript[λ, 4])(t)(1/(r(t))),

(Subscript[λ, 4])'(t)=0.

From () and (), (Subscript[λ, 4])≡0 and thus () to () become

H=(Subscript[λ, 1])u(t)+(Subscript[λ, 2])((((v^2)(t))/(r(t)))-(μ/((r^2)(t)))+((T sin β(t))/((Subscript[m, 0])-m't)))+(Subscript[λ, 3])(-((u(t)v(t))/(r(t)))+((T cos β(t))/((Subscript[m, 0])-m't)) ),

(Subscript[λ, 1])'(t)= -(Subscript[λ, 2])(t)(-(((v^2)(t))/((r^2)(t)))+((2μ)/((r^3)(t))))-(Subscript[λ, 3])(t)(((u(t)v(t))/((r^2)(t))) ),

(Subscript[λ, 2])'(t)=-(Subscript[λ, 1])(t)+(Subscript[λ, 3])(t)((v(t))/(r(t))),

(Subscript[λ, 3])'(t)=-(Subscript[λ, 2])(t)((2v(t))/(r(t)))+(Subscript[λ, 3])(t)((u(t))/(r(t))).

Furthermore, from (b) in theorem 2,

0=((∂H)/(∂β))=(T/((Subscript[m, 0])-m't))((Subscript[λ, 2])(t) cos β(t)-(Subscript[λ, 1])(t) sin β(t)) ⟹ β(t)=arctan (((Subscript[λ, 2])(t))/((Subscript[λ, 3])(t))).

Thus the control function β is determined by the multipliers (Subscript[λ, 2]) and (Subscript[λ, 3]).

Based on (), note that the polar angle θ is determined by v and r.

From the transversality relation (d) in theorem 2 (i.e. equation ()),

(Subscript[λ, 1])(b)=1.

This yields a system of six nonlinear differential equations (), (), (), (), (), and () in the variables  r,  u,  v,  (Subscript[λ, 1]), (Subscript[λ, 2]), and (Subscript[λ, 3]) with six bilocal conditions (), (), (), (), (), and ().

As mentioned earlier, the variables β and θ follow.

The next section implements a dynamical approach to the maximal orbit transfer problem.

A Dynamic Approach to the Maximal Orbit Transfer Problem

The function MaximalRadiusOrbitTransfer dynamically shows the maximal radius orbit transfer between two coplanar circular orbits so that their centers are located at a single center of attraction. Here thrust is the constant thrust of the engine, dmr is the decreasing mass rate due to the constant propellant flow rate, b is the final time, m0 is the initial mass of the spacecraft including the propellant, μ is the gravitational constant, r0 is the initial radius, u0 is the initial radial velocity, ub is the final radial velocity, v0 is the initial tangential velocity, and k is the number of thrust vectors.

Clearly the problem is nonlinear and, to the author’s knowledge, no closed-form solution has been found. The possibility of obtaining a solution through a numerical method remains, as implied by theorem 1. The accuracy of the results depends sensitively on the initial values. The Method option is needed for Mathematica 9 or lower; for faster processing, remove it in Mathematica 10 or higher.

MaximalRadiusOrbitTransfer[thrust_,dmr_,b_,r0_,μ_,m0_,u0_,ub_,k_]:=Module[{v0,t,r,u,v,λr,λv,λu,sol,R,V,θ,θmax},
(* data for the flight *)
v0=Sqrt[μ/r0];       (* initial tangential velocity; tangent to the inner circle *)
sol=NDSolve[{             (* solution of the problem *)
r'[t]==u[t],
u'[t]==v[t]^2/r[t]-μ/r[t]^2+thrust/(m0-dmr t) λu[t]/Sqrt[λu[t]^2+λv[t]^2],v'[t]==-((u[t] v[t])/r[t])+thrust/(m0-dmr t) λv[t]/Sqrt[λu[t]^2+λv[t]^2],λr'[t]==λu[t] (v[t]^2/r[t]^2-2.0 μ/r[t]^3)-λv[t] (u[t] v[t])/r[t]^2,
λu'[t]==-λr[t]+λv[t] v[t]/r[t],
λv'[t]==-2.0 λu[t] v[t]/r[t]+λv[t] u[t]/r[t],
r[0]==r0,
u[0]==u0,
v[0]==v0,
λr[b]==1.0,
u[b]==ub,
v[b]==Sqrt[μ/r[b]]},
{r,u,v,λr,λu,λv},{t,0,b},Method->{"Shooting","StartingInitialConditions"->{r[0]==r0,u[0]==u0,v[0]==v0,λr[0]==1.0,λu[0]==-0.5,λv[0]==0.0}}];
R[t_]=Evaluate[r[t]/.sol][[1]];
V[t_]=Evaluate[v[t]/.sol][[1]];
(* the orbit transfer figure *) 
θ[time_]:=NIntegrate[Interpolation[Table[{n b/k,V[n b/k]/R[n b/k]},{n,0,k}],InterpolationOrder->1][t],{t,0,time}];
θmax=180 θ[b]/N[Pi];
Show[
Graphics[{
Text[Style[Row[{"final radius = ",R[b]}],12],{0,.36r0}],
Text[Style[Row[{"final θ = ",Superscript[θmax," ∘"]}],12],{0,.19r0}],
{EdgeForm[Thin],Red,Disk[{0,0},0.04]},
{Red,Circle[{0,0},r0,{0,2 Pi}]},
{Green,Circle[{0,0},R[b],{0,2 Pi}]},
{Red,Arrowheads[0.02],Arrow[Table[{R[n b/k] {Cos[θ[n b/k]],Sin[θ[n b/k]]},R[n b/k] {Cos[θ[n b/k]],Sin[θ[n b/k]]}+thrust {Sin[θ[n b/k]],Cos[θ[n b/k]]}},{n,0,k}]]}
}],
ListLinePlot[Table[R[n b/k] {Cos[θ[n b/k]],Sin[θ[n b/k]]},{n,0,k}],PlotStyle->{Blue,Thick}],
Axes->False,PlotRange->All,ImageSize->{400,400},AxesOrigin->{0,0}]
]

Manipulate[
Quiet@MaximalRadiusOrbitTransfer[thrust,dmr,b,1.0,1.0,1.0,0.0,0.0,k],
{{thrust,0.1405},0.14,0.15,0.0001,Appearance->"Labeled"},
{{dmr,0.07485,"decreasing mass rate"},0.07,0.08,0.001,Appearance->"Labeled"},
{{b,3.317291,"final time interval"},3.2,3.5,0.01,Appearance->"Labeled"},
{{k,20,"number of thrust arrows"},15,50,1,Appearance->"Labeled"},
SaveDefinitions->True,TrackedSymbols->{thrust,dmr,b,k},SynchronousUpdating->False
]

Acknowledgments

The author expresses his deep gratitude to Prof. Dr., Dr. H. C., Heiner H. Gonska of the University of Duisburg-Essen, Faculty of Mathematics, located in Duisburg, Germany, for his kind invitation and warm hospitality. The invitation was funded by a grant from the Center of Excellence for Applications of Mathematics supported by DAAD. The author also thanks Aida Viziru for her rapid and professional help.

References

A. E. Bryson and Y.-C. Ho, Applied Optimal Control: Optimization, Estimation, and Control, New York: Halsted Press, 1975.

S. Alfano and J. D. Thorne, Constant-Thrust Orbit-Raising Transfer Charts, Report PL-TR-93-1010, July 1993. www.dtic.mil/dtic/tr/fulltext/u2/a269088.pdf.

S. Alfano and J. D. Thorne, “Circle-to-Circle Constant-Thrust Orbit Raising,” The Journal of the Astronautical Sciences, 42(1), 1994 pp. 35–45. www.centerforspace.com/downloads/files/pubs/JAS.V42.N01.pdf.

L. Cesari, Optimization—Theory and Applications, Problems with Ordinary Differential Equations, Applications of Mathematics, Vol. 17, New York: Springer, 1983.

M. Mureşan, A Concrete Approach to Classical Analysis, New York: Springer, 2009.

R. Vinter, Optimal Control, Systems & Control: Foundations & Applications, Boston: Birkhäuser, 2000.