Low-Thrust Orbit Transfer#

Mathematical Formulation#

Problem Statement#

Find the optimal thrust direction \(\mathbf{u}(t) = [u_1, u_2, u_3]^T\) and thrust efficiency parameter \(\tau\) that maximize the final spacecraft mass:

(1)#\[J = -w(t_f)\]

Subject to the Gauss variational equations in modified equinoctial elements:

\[\frac{dp}{dt} = \frac{2p}{q}\sqrt{\frac{p}{\mu}} \delta_2\]
\[\frac{df}{dt} = \sqrt{\frac{p}{\mu}} \left[ \sin L \cdot \delta_1 + \frac{1}{q}((q+1)\cos L + f) \delta_2 - \frac{g}{q}(h \sin L - k \cos L) \delta_3 \right]\]
\[\frac{dg}{dt} = \sqrt{\frac{p}{\mu}} \left[ -\cos L \cdot \delta_1 + \frac{1}{q}((q+1)\sin L + g) \delta_2 + \frac{f}{q}(h \sin L - k \cos L) \delta_3 \right]\]
\[\frac{dh}{dt} = \sqrt{\frac{p}{\mu}} \frac{s^2 \cos L}{2q} \delta_3\]
\[\frac{dk}{dt} = \sqrt{\frac{p}{\mu}} \frac{s^2 \sin L}{2q} \delta_3\]
\[\frac{dL}{dt} = \sqrt{\frac{\mu p}{p^3}} q^2 + \sqrt{\frac{p}{\mu}} \frac{1}{q}(h \sin L - k \cos L) \delta_3\]
\[\frac{dw}{dt} = -\frac{T(1 + 0.01\tau)}{I_{sp}}\]

Where:

  • \(q = 1 + f \cos L + g \sin L\)

  • \(s^2 = 1 + h^2 + k^2\)

  • \(\delta_i = \Delta g_i + \Delta T_i\) (gravitational + thrust accelerations)

Thrust and Gravitational Accelerations#

Thrust accelerations:

\[\Delta T_1 = \frac{g_0 T (1 + 0.01\tau)}{w} u_1\]
\[\Delta T_2 = \frac{g_0 T (1 + 0.01\tau)}{w} u_2\]
\[\Delta T_3 = \frac{g_0 T (1 + 0.01\tau)}{w} u_3\]

Gravitational perturbations (J2, J3, J4):

\[\Delta g_1, \Delta g_2, \Delta g_3\]
computed from Earth’s gravitational harmonics

Terminal Constraints#

\[p(t_f) = p_{tf}\]
\[\sqrt{f^2(t_f) + g^2(t_f)} = e_{tf}\]
\[\sqrt{h^2(t_f) + k^2(t_f)} = \tan(i_{tf}/2)\]
\[f(t_f) h(t_f) + g(t_f) k(t_f) = 0\]
\[-10 \leq g(t_f) h(t_f) - k(t_f) f(t_f) \leq 0\]

Boundary Conditions#

  • Initial conditions: \(p(0) = 21837080.052835\) ft, \(f(0) = 0\), \(g(0) = 0\), \(h(0) = -0.25396764647494\), \(k(0) = 0\), \(L(0) = \pi\), \(w(0) = 1\)

  • Final conditions: Determined by terminal constraints

  • Control bounds: \(u_1^2 + u_2^2 + u_3^2 = 1\) (unit thrust vector)

  • Parameter bounds: \(-50 \leq \tau \leq 0\)

Physical Parameters#

  • Specific impulse: \(I_{sp} = 450\) s

  • Earth gravitational parameter: \(\mu = 1.407645794 \times 10^{16}\) ft³/s²

  • Standard gravity: \(g_0 = 32.174\) ft/s²

  • Thrust magnitude: \(T = 4.446618 \times 10^{-3}\) lbf

  • Earth radius: \(R_E = 20925662.73\) ft

  • Gravitational harmonics: \(J_2 = 1082.639 \times 10^{-6}\), \(J_3 = -2.565 \times 10^{-6}\), \(J_4 = -1.608 \times 10^{-6}\)

State Variables#

  • \(p(t)\): Semi-latus rectum (ft)

  • \(f(t), g(t)\): Eccentricity vector components

  • \(h(t), k(t)\): Inclination vector components

  • \(L(t)\): True longitude (rad)

  • \(w(t)\): Spacecraft mass (normalized)

Control Variables#

  • \(u_1(t), u_2(t), u_3(t)\): Thrust direction unit vector components

Static Parameter#

  • \(\tau\): Thrust efficiency parameter

Notes#

This problem uses modified equinoctial elements to avoid singularities in classical orbital elements and includes Earth’s gravitational harmonics (J2, J3, J4) for high-fidelity orbital mechanics.

Running This Example#

cd examples/low_thrust_orbit_transfer
python low_thrust_orbit_transfer.py

Code Implementation#

examples/low_thrust_orbit_transfer/low_thrust_orbit_transfer.py#
  1import casadi as ca
  2import numpy as np
  3from initial_guess_generator import generate_initial_guess
  4
  5import maptor as mtor
  6
  7
  8# ============================================================================
  9# Physical Constants and Parameters
 10# ============================================================================
 11
 12# Propulsion parameters
 13ISP = 450.0
 14G0 = 32.174
 15T_THRUST = 4.446618e-3
 16
 17# Gravitational parameters
 18MU = 1.407645794e16
 19RE = 20925662.73
 20J2 = 1082.639e-6
 21J3 = -2.565e-6
 22J4 = -1.608e-6
 23
 24
 25# ============================================================================
 26# Scaling Factors
 27# ============================================================================
 28
 29P_SCALE = 1e7
 30L_SCALE = np.pi
 31T_SCALE = 1e4
 32
 33
 34# ============================================================================
 35# Initial and Target Conditions
 36# ============================================================================
 37
 38# Initial modified equinoctial elements
 39PTI = 21837080.052835
 40FTI = 0.0
 41GTI = 0.0
 42HTI = -0.25396764647494
 43KTI = 0.0
 44LTI = np.pi
 45WTI = 1.0
 46
 47# Target conditions
 48PTF = 40007346.015232
 49EVENT_FINAL_9 = 0.73550320568829
 50EVENT_FINAL_10 = 0.61761258786099
 51EVENT_FINAL_11 = 0.0
 52EVENT_FINAL_12_LOWER = -10.0
 53EVENT_FINAL_12_UPPER = 0.0
 54
 55# Numerical tolerance
 56EQ_TOL = 0.001
 57
 58
 59# ============================================================================
 60# Helper Functions
 61# ============================================================================
 62
 63
 64def cross_product(a, b):
 65    return ca.vertcat(
 66        a[1] * b[2] - a[2] * b[1],
 67        a[2] * b[0] - a[0] * b[2],
 68        a[0] * b[1] - a[1] * b[0],
 69    )
 70
 71
 72def dot_product(a, b):
 73    return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]
 74
 75
 76# ============================================================================
 77# Problem Setup
 78# ============================================================================
 79
 80problem = mtor.Problem("Low-Thrust Orbit Transfer")
 81phase = problem.set_phase(1)
 82
 83
 84# ============================================================================
 85# Variables
 86# ============================================================================
 87
 88# Time variable (scaled)
 89t_s = phase.time(initial=0.0, final=(50000.0 / T_SCALE, 100000.0 / T_SCALE))
 90
 91# State variables (modified equinoctial elements)
 92p_s = phase.state(
 93    "p_scaled",
 94    initial=PTI / P_SCALE,
 95    boundary=(10.0e6 / P_SCALE, 60.0e6 / P_SCALE),
 96)
 97f = phase.state("f", initial=FTI, boundary=(-0.20, 0.20))
 98g = phase.state("g", initial=GTI, boundary=(-0.10, 1.0))
 99h = phase.state("h", initial=HTI, boundary=(-1.0, 1.0))
100k = phase.state("k", initial=KTI, boundary=(-0.20, 0.20))
101L_s = phase.state(
102    "L_scaled",
103    initial=LTI / L_SCALE,
104    boundary=(np.pi / L_SCALE, 20 * np.pi / L_SCALE),
105)
106w = phase.state("w", initial=WTI, boundary=(0.0, 2.0))
107
108# Control variables (thrust direction components)
109u1 = phase.control("u1", boundary=(-1.0, 1.0))
110u2 = phase.control("u2", boundary=(-1.0, 1.0))
111u3 = phase.control("u3", boundary=(-1.0, 1.0))
112
113# Parameter (thrust efficiency factor)
114tau = problem.parameter("tau", boundary=(-50.0, 0.0))
115
116
117# ============================================================================
118# Mathematical Formulation
119# ============================================================================
120
121# Convert scaled variables to physical units
122t = t_s * T_SCALE
123p = p_s * P_SCALE
124L = L_s * L_SCALE
125
126# Numerical safeguard
127eps = 1e-12
128
129# Modified equinoctial element calculations
130q = 1.0 + f * ca.cos(L) + g * ca.sin(L)
131r = p / ca.fmax(q, eps)
132alpha2 = h * h - k * k
133X = ca.sqrt(h * h + k * k)
134s2 = 1 + X * X
135
136# Position vector components
137r1 = r / s2 * (ca.cos(L) + alpha2 * ca.cos(L) + 2 * h * k * ca.sin(L))
138r2 = r / s2 * (ca.sin(L) - alpha2 * ca.sin(L) + 2 * h * k * ca.cos(L))
139r3 = 2 * r / s2 * (h * ca.sin(L) - k * ca.cos(L))
140rvec = ca.vertcat(r1, r2, r3)
141
142# Velocity vector components
143sqrt_mu_p = ca.sqrt(MU / ca.fmax(p, eps))
144v1 = (
145    -(1.0 / s2)
146    * sqrt_mu_p
147    * (ca.sin(L) + alpha2 * ca.sin(L) - 2 * h * k * ca.cos(L) + g - 2 * f * h * k + alpha2 * g)
148)
149v2 = (
150    -(1.0 / s2)
151    * sqrt_mu_p
152    * (-ca.cos(L) + alpha2 * ca.cos(L) + 2 * h * k * ca.sin(L) - f + 2 * g * h * k + alpha2 * f)
153)
154v3 = (2.0 / s2) * sqrt_mu_p * (h * ca.cos(L) + k * ca.sin(L) + f * h + g * k)
155vvec = ca.vertcat(v1, v2, v3)
156
157# Reference frame vectors
158rv = cross_product(rvec, vvec)
159rvr = cross_product(rv, rvec)
160norm_r = ca.sqrt(ca.fmax(dot_product(rvec, rvec), eps))
161norm_rv = ca.sqrt(ca.fmax(dot_product(rv, rv), eps))
162
163ir = rvec / norm_r
164ith = rvr / (norm_rv * norm_r)
165ih = rv / norm_rv
166
167# Earth-fixed reference calculations
168en = ca.vertcat(0.0, 0.0, 1.0)
169enir = dot_product(en, ir)
170in_vec = en - enir * ir
171norm_in = ca.sqrt(ca.fmax(dot_product(in_vec, in_vec), eps))
172in_normalized = in_vec / norm_in
173
174# Gravitational perturbation calculations
175sin_phi = rvec[2] / norm_r
176sin_phi_clamped = ca.fmax(ca.fmin(sin_phi, 1.0 - eps), -1.0 + eps)
177cos_phi = ca.sqrt(1.0 - sin_phi_clamped**2)
178
179# J2 perturbation
180P2 = 0.5 * (3.0 * sin_phi_clamped**2 - 1.0)
181Pdash2 = 3.0 * sin_phi_clamped
182r_safe = ca.fmax(r, RE / 100.0)
183deltagn_j2 = -MU * cos_phi / (r_safe * r_safe) * (RE / r_safe) ** 2 * Pdash2 * J2
184deltagr_j2 = -MU / (r_safe * r_safe) * 3.0 * (RE / r_safe) ** 2 * P2 * J2
185
186# J3 perturbation
187P3 = 0.5 * (5.0 * sin_phi_clamped**3 - 3.0 * sin_phi_clamped)
188Pdash3 = 0.5 * (15.0 * sin_phi_clamped**2 - 3.0)
189deltagn_j3 = -MU * cos_phi / (r_safe * r_safe) * (RE / r_safe) ** 3 * Pdash3 * J3
190deltagr_j3 = -MU / (r_safe * r_safe) * 4.0 * (RE / r_safe) ** 3 * P3 * J3
191
192# J4 perturbation
193P4 = (1.0 / 8.0) * (35.0 * sin_phi_clamped**4 - 30.0 * sin_phi_clamped**2 + 3.0)
194Pdash4 = (1.0 / 8.0) * (140.0 * sin_phi_clamped**3 - 60.0 * sin_phi_clamped)
195deltagn_j4 = -MU * cos_phi / (r_safe * r_safe) * (RE / r_safe) ** 4 * Pdash4 * J4
196deltagr_j4 = -MU / (r_safe * r_safe) * 5.0 * (RE / r_safe) ** 4 * P4 * J4
197
198# Total gravitational perturbations
199deltagn = deltagn_j2 + deltagn_j3 + deltagn_j4
200deltagr = deltagr_j2 + deltagr_j3 + deltagr_j4
201delta_g = deltagn * in_normalized - deltagr * ir
202
203# Gravitational acceleration components
204DELTA_g1 = dot_product(ir, delta_g)
205DELTA_g2 = dot_product(ith, delta_g)
206DELTA_g3 = dot_product(ih, delta_g)
207
208# Thrust acceleration components
209w_safe = ca.fmax(w, eps)
210DELTA_T1 = G0 * T_THRUST * (1.0 + 0.01 * tau) * u1 / w_safe
211DELTA_T2 = G0 * T_THRUST * (1.0 + 0.01 * tau) * u2 / w_safe
212DELTA_T3 = G0 * T_THRUST * (1.0 + 0.01 * tau) * u3 / w_safe
213
214# Total acceleration components
215delta1 = DELTA_g1 + DELTA_T1
216delta2 = DELTA_g2 + DELTA_T2
217delta3 = DELTA_g3 + DELTA_T3
218
219
220# ============================================================================
221# Dynamics
222# ============================================================================
223
224# Gauss variational equations for modified equinoctial elements
225sqrt_p_mu = ca.sqrt(ca.fmax(p, eps) / MU)
226q_safe = ca.fmax(ca.fabs(q), eps)
227
228pdot = 2 * p / q_safe * sqrt_p_mu * delta2
229fdot = (
230    sqrt_p_mu * ca.sin(L) * delta1
231    + sqrt_p_mu * (1.0 / q_safe) * ((q + 1.0) * ca.cos(L) + f) * delta2
232    - sqrt_p_mu * (g / q_safe) * (h * ca.sin(L) - k * ca.cos(L)) * delta3
233)
234gdot = (
235    -sqrt_p_mu * ca.cos(L) * delta1
236    + sqrt_p_mu * (1.0 / q_safe) * ((q + 1.0) * ca.sin(L) + g) * delta2
237    + sqrt_p_mu * (f / q_safe) * (h * ca.sin(L) - k * ca.cos(L)) * delta3
238)
239hdot = sqrt_p_mu * s2 * ca.cos(L) / (2.0 * q_safe) * delta3
240kdot = sqrt_p_mu * s2 * ca.sin(L) / (2.0 * q_safe) * delta3
241Ldot = (
242    sqrt_p_mu * (1.0 / q_safe) * (h * ca.sin(L) - k * ca.cos(L)) * delta3
243    + ca.sqrt(MU * ca.fmax(p, eps)) * (q / ca.fmax(p, eps)) ** 2
244)
245wdot = -T_THRUST * (1.0 + 0.01 * tau) / ISP
246
247# Scaled dynamics
248phase.dynamics(
249    {
250        p_s: (pdot / P_SCALE) * T_SCALE,
251        f: fdot * T_SCALE,
252        g: gdot * T_SCALE,
253        h: hdot * T_SCALE,
254        k: kdot * T_SCALE,
255        L_s: (Ldot / L_SCALE) * T_SCALE,
256        w: wdot * T_SCALE,
257    }
258)
259
260
261# ============================================================================
262# Constraints
263# ============================================================================
264
265# Path constraint: Unit thrust vector magnitude
266thrust_magnitude_squared = u1**2 + u2**2 + u3**2
267phase.path_constraints(
268    thrust_magnitude_squared >= 1.0 - EQ_TOL,
269    thrust_magnitude_squared <= 1.0 + EQ_TOL,
270)
271
272# Event constraints: Final orbital element targets
273phase.event_constraints(
274    p_s.final == PTF / P_SCALE,
275    ca.sqrt(f.final**2 + g.final**2) == EVENT_FINAL_9,
276    ca.sqrt(h.final**2 + k.final**2) == EVENT_FINAL_10,
277    f.final * h.final + g.final * k.final == EVENT_FINAL_11,
278)
279
280# Additional final constraint
281gtf_htf_minus_ktf_ftf = g.final * h.final - k.final * f.final
282phase.event_constraints(
283    gtf_htf_minus_ktf_ftf >= EVENT_FINAL_12_LOWER,
284    gtf_htf_minus_ktf_ftf <= EVENT_FINAL_12_UPPER,
285)
286
287
288# ============================================================================
289# Objective
290# ============================================================================
291
292# Maximize final mass (minimize negative mass)
293problem.minimize(-w.final)
294
295
296# ============================================================================
297# Mesh Configuration and Initial Guess
298# ============================================================================
299
300# Mesh configuration
301phase.mesh(
302    [8, 8, 8, 8, 8, 8, 8, 8],
303    [-1.0, -6 / 7, -4 / 7, -2 / 7, 0, 2 / 7, 4 / 7, 6 / 7, 1.0],
304)
305
306# Generate initial guess from external module
307states_guess, controls_guess, final_time_guess = generate_initial_guess()
308
309phase.guess(
310    states=states_guess,
311    controls=controls_guess,
312    terminal_time=final_time_guess,
313)
314
315problem.parameter_guess(tau=-25.0)
316
317# ============================================================================
318# Solve
319# ============================================================================
320
321solution = mtor.solve_adaptive(
322    problem,
323    error_tolerance=1e-4,
324    max_iterations=30,
325    min_polynomial_degree=3,
326    max_polynomial_degree=8,
327    ode_solver_tolerance=1e-4,
328    nlp_options={
329        "ipopt.print_level": 0,
330        "ipopt.max_iter": 500,
331        "ipopt.tol": 1e-4,
332        "ipopt.constr_viol_tol": 1e-4,
333        "ipopt.acceptable_tol": 1e-3,
334        "ipopt.mu_strategy": "adaptive",
335        "ipopt.linear_solver": "mumps",
336    },
337)
338
339
340# ============================================================================
341# Results
342# ============================================================================
343
344if solution.status["success"]:
345    final_mass = solution[(1, "w")][-1]
346    final_time_scaled = solution.phases[1]["times"]["final"]
347    final_time = final_time_scaled * T_SCALE
348
349    print(f"Final mass: {final_mass:.6f}")
350    print(f"Final time: {final_time:.1f} seconds ({final_time / 3600:.2f} hours)")
351
352    p_final_scaled = solution[(1, "p_scaled")][-1]
353    p_final_physical = p_final_scaled * P_SCALE
354    print(f"Final p (scaled): {p_final_scaled:.6f}")
355    print(f"Final p (physical): {p_final_physical:.1f} ft")
356
357    solution.plot()
358else:
359    print(f"Failed: {solution.status['message']}")