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)
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)
In [ ]: