# Circular Restricted Three Body Problem Derrivation

SymPy is usefull for: 1) Derriving the full equations of motion 2) Optimizing code through common subexpression elimination 3) Generating complicated $\LaTeX$ expressions. SymPy is open-source.



In [1]:

from sympy import *




In [2]:

# State
x, y, z, vx, vy, vz, m = symbols('x y z vx vy vz m')
r                      = Matrix([x, y, z])
v                      = Matrix([vx, vy, vz])
s                      = Matrix([r, v, [m]])

# Costate
lx, ly, lz, lvx, lvy, lvz, lm = symbols('lx ly lz lvx lvy lvz lm')
lr                            = Matrix([lx, ly, lz])
lv                            = Matrix([lvx, lvy, lvz])
l                             = Matrix([lr, lv, [lm]])

# Full State
fs = Matrix([s, l])

# Controls
u, ax, ay, az = symbols('u ax ay az')
a             = Matrix([ax, ay, az])
cont          = Matrix([[u], a])

# Parameters
Isp, g0, mu, T, eps = symbols('Isp g0 mu T eps')
c                   = Isp*g0




In [3]:

# Intermediates
r1 = sqrt((x + mu)**2 + y**2 + z**2)
r2 = sqrt((x + mu -1)**2 + y**2 + z**2)




In [4]:

# Position Dependence
g = Matrix([
x - (1 - mu)*(x + mu)/r1**3 - mu*(x + mu - 1)/r2**3,
y - (1 - mu)*y/r1**3 - mu*y/r2**3,
-(1 - mu)*z/r1**3 - mu*z/r2**3
])
# Velocity Dependence
h = Matrix([2*vy, -2*vx, 0])




In [5]:

# State Equations of Motion
ds = Matrix([v, g + h + u*T*a/m, [-u*T/c]])




In [ ]:

# Homotopic Cost Lagrangian
L = T*(u - eps*u*(1 - u))/c




In [ ]:

# Hamiltonian
Hamiltonian = l.dot(ds) + L




In [ ]:

# Location Dependence Derrivative
G = g.jacobian(r)
# Velocity Dependence Derrivative
H = h.jacobian(v)




In [ ]:

# Costate Equations of Motion
dl = -Matrix([Hamiltonian.diff(i) for i in s])




In [ ]:

# Fullstate Equations of Motion
dfs = Matrix([ds, dl])




In [ ]:

# Switching Function
S = -lv.norm()*c/m - lm + 1




In [ ]:

# Optimal Thrust Direction
aopt = -lv.normalized()




In [ ]:

# Derrivative of Fullstate
ddfs = dfs.jacobian(fs)



# Implementation

Initialisation parameters sourced from Low-Thrust Minimum Fuel Optimization in the Circular Restricted Three-Body Problem (Zhang et. al.)

The nonlinear programming problem's decision vector is $$[t_f, \lambda_x(t_i), \lambda_y(t_i), \lambda_z(t_i), \lambda_{v_x}(t_i), \lambda_{v_y}(t_i), \lambda_{v_z}(t_i), \lambda_m(t_i)]$$



In [1]:

import sys
sys.path.append('../')
from trajectory.space import CRTBP
from numpy import *
from PyGMO import *




In [2]:

# Parameters
mu  = 1.21506683e-2 # Gravitational parameter
T   = 10. # Maximum thrust [N]
Isp = 2000. # Specific impulse
eps = 1. # Homotopy paramter (energy-minimisation)




In [3]:

# Initial state
si = array([-0.019488511458668, -0.016033479812051, 0,
8.918881923678198, -4.081793688818725, 0,
1000], float)
# Target state
st = array([0.823385182067467, 0, -0.022277556273235,
0, 0.134184170262437, 0,
1000], float)
# Decision variable guess
di = array([8.6,                            # Time of flight [days]
5.616017, 32.875896, -0.094522, # Position costates
-0.101606, 0.044791, -0.000150, # Velocity costates
0.133266], float)               # Mass costate




In [4]:

# Intialise the problem
prob = CRTBP(mu, T, Isp, eps, si, st)




In [5]:

algo = algorithm.scipy_slsqp(max_iter=20, screen_output=True)
algo = algorithm.mbh(algo, screen_output=True)




In [ ]:

pop = population(prob, 1)




/usr/lib/python2.7/dist-packages/scipy/integrate/_ode.py:1035: UserWarning: dop853: larger nmax is needed
self.messages.get(idid, 'Unexpected idid=%s' % idid))




In [ ]: