Notebook source code:
examples/kepler/kepler.ipynb
Run the notebook yourself on binder
Kepler (minimum time orbit transfer) - Python code¶
Minimum time control of the Kepler equation (CNES / TAS / Inria / CNRS collaboration):
\[t_f \to \min,\]
\[\ddot{q} = -\mu\frac{q}{|q|^3}+\frac{u}{m}\,,\quad t \in [0,t_f],\]
\[\dot{m} = -\beta|u|,\quad |u| \leq T_{\mathrm{max}}.\]
Fixed initial and final Keplerian orbits (free final longitude).
Initializations¶
[1]:
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import time
from nutopy import nle
from nutopy import tools
from nutopy import ocp
ctmax = (3600**2) / 1e6 # Conversion from Newtons
mass0 = 1500. # Initial mass of the spacecraft
beta = 1.42e-02 # Engine specific impulsion
mu = 5165.8620912 # Earth gravitation constant
t0 = 0. # Initial time (final time is free)
x0 = np.array([ 11.625, 0.75, 0., 6.12e-02, 0., 3.14159265358979 ]) # Initial state (fixed initial longitude)
xf_fixed = np.array([ 42.165, 0., 0., 0., 0. ]) # Final state (free final longitude)
# tmax = 60 Newtons
#tmax = ctmax * 60.; tf = 15.2055; p0 = -np.array([ .361266, 22.2412, 7.87736, 0., 0., -5.90802 ]); N = 1000
# tmax = 6 Newtons
tmax = ctmax * 6.; tf = 1.32e2; p0 = -np.array([ -4.743728539366440e+00, -7.171314869854240e+01, -2.750468309804530e+00, 4.505679923365745e+01, -3.026794475592510e+00, 2.248091067047670e+00 ]); N = 1000
# tmax = 0.7 Newtons
#tmax = ctmax * 0.7; tf = 1.210000000000000e+03; p0 = -np.array([ -2.215319700438820e+01, -4.347109477345140e+01, 9.613188807286992e-01, 3.181800985503019e+02, -2.307236094862410e+00, -5.797863110671591e-01 ]); N = 5000
# tmax = 0.14 Newtons
#tmax = ctmax * 0.14; tf = 6.08e3; p0 = -np.array([ -1.234155379067110e+02, -6.207170881591489e+02, 5.742554220129187e-01, 1.629324243017332e+03, -2.373935935351530e+00, -2.854066853269850e-01 ]) ; N = 5000
p0 = p0 / np.linalg.norm(p0) # Normalization |p0|=1 for free final time
y = np.hstack((p0, tf)) # initial guess, y = (p0, tf)
Hamiltonian (fortran wrapper)¶
The first and second derivatives of the code are generated by Tapenade. Important: arguments dx and dx0 (and dp and dp0) have been exchanged; the Tapenade generated signature was
SUBROUTINE HFUN_D_D(t, x, xd0, xd, p, pd0, pd, tmax, mass0, beta, mu, h, hd, hdd, n)
while the corrected one is
SUBROUTINE HFUN_D_D(t, x, xd, xd0, p, pd, pd0, tmax, mass0, beta, mu, h, hd, hdd, n)
This ensures that the signature matches what is expected by @tensorize: variations, up to any order, come as x, dx, d2x, d3x, ...
[2]:
!python -m numpy.f2py -c hfun.f90 -m hfun > /dev/null 2>&1
!python -m numpy.f2py -c hfun_d.f90 -m hfun_d > /dev/null 2>&1
!python -m numpy.f2py -c hfun_d_d.f90 -m hfun_d_d > /dev/null 2>&1
from hfun import hfun
from hfun_d import hfun_d
from hfun_d_d import hfun_d_d
hfun = tools.tensorize(hfun_d, hfun_d_d, tvars=(2, 3), full=True)(hfun)
h = ocp.Hamiltonian(hfun)
f = ocp.Flow(h)
Shooting function¶
[3]:
def dshoot(t0, dt0, x0, dx0, p0, dp0, tf, dtf, next=False):
(xf, dxf), (pf, dpf) = f((t0, dt0), (x0, dx0), (p0, dp0), (tf, dtf), tmax, mass0, beta, mu)
s = np.zeros(7) # code duplication and full=True
s[0:5] = xf[0:5] - xf_fixed
s[5] = pf[5] # free final longitude
s[6] = p0[0]**2 + p0[1]**2 + p0[2]**2 + p0[3]**2 + p0[4]**2 + p0[5]**2 - 1.
ds = np.zeros(7)
ds[0:5] = dxf[0:5]
ds[5] = dpf[5]
ds[6] = 2*p0[0]*dp0[0] + 2*p0[1]*dp0[1] + 2*p0[2]*dp0[2] + 2*p0[3]*dp0[3] + 2*p0[4]*dp0[4] + 2*p0[5]*dp0[5]
if not next: return s, ds
else: return s, ds, ((tf, dtf), (xf, dxf), (pf, dpf), None)
@tools.vectorize(vvars=(3,))
@tools.vectorize(vvars=(4,), next=True)
@tools.tensorize(dshoot, full=True)
def shoot(t0, x0, p0, tf, next=False):
"""s = shoot(t0, x0, p0, tf)
Shooting function associated with h
"""
xf, pf = f(t0, x0, p0, tf, tmax, mass0, beta, mu)
s = np.zeros(7)
s[0:5] = xf[0:5] - xf_fixed
s[5] = pf[5] # free final longitude
s[6] = p0[0]**2 + p0[1]**2 + p0[2]**2 + p0[3]**2 + p0[4]**2 + p0[5]**2 - 1.
if not next: return s
else: return s, (tf, xf, pf, None)
Solve¶
[4]:
dfoo = lambda y, dy: shoot(t0, x0, (y[:-1], dy[:-1]), (y[-1], dy[-1]))
foo = lambda y: shoot(t0, x0, y[:-1], y[-1])
foo = tools.tensorize(dfoo, full=True)(foo)
nleopt = nle.Options(SolverMethod='hybrj', Display='on', TolX=1e-8)
et = time.time(); sol = nle.solve(foo, y, df=foo, options=nleopt); y_sol = sol.x; et = time.time() - et
print('Elapsed time:', et)
print('y_sol =', y_sol)
print('foo =', foo(y_sol))
Calls |f(x)| |x|
1 6.288197653086673e+00 1.320037878244408e+02
2 7.100017914373875e-01 1.398218271582472e+02
3 1.971014432778390e+00 1.411556611371509e+02
4 2.079137312804174e+00 1.419381872915138e+02
5 1.698641948837419e-01 1.411913913720076e+02
6 5.530356739475702e-02 1.412436418925176e+02
7 5.160467577766739e-03 1.412343505428431e+02
8 3.609540041386610e-03 1.412315324616509e+02
9 2.304328588779287e-03 1.412296626257120e+02
10 3.861188438084359e-05 1.412304710130267e+02
11 3.444274156379021e-06 1.412304911883799e+02
12 1.223115342687875e-06 1.412304904958646e+02
13 3.838113948682708e-08 1.412304906089938e+02
14 1.469368913720655e-08 1.412304906245045e+02
15 4.686170678922632e-09 1.412304906240111e+02
16 1.669128031059662e-11 1.412304906238177e+02
Results of the nle solver method:
xsol = [ 5.88696798e-02 8.12460412e-01 1.21287152e-02 -5.79655506e-01
1.39917457e-02 -9.95283135e-03 1.41226950e+02]
f(xsol) = [-1.65982783e-11 -2.31574742e-13 -1.34379389e-12 6.69900333e-14
-4.87884632e-13 -1.01267627e-13 9.91873250e-13]
nfev = 16
njev = 2
status = 1
success = True
Successfully completed: relative error between two consecutive iterates is at most TolX.
Elapsed time: 14.979362487792969
y_sol = [ 5.88696798e-02 8.12460412e-01 1.21287152e-02 -5.79655506e-01
1.39917457e-02 -9.95283135e-03 1.41226950e+02]
foo = [-1.65982783e-11 -2.31574742e-13 -1.34379389e-12 6.69900333e-14
-4.87884632e-13 -1.01267627e-13 9.91873250e-13]
Plots¶
[5]:
p0 = y_sol[:-1]
tf = y_sol[-1]
tspan = list(np.linspace(0, tf, N+1))
xf, pf = np.array( f(t0, x0, p0, tspan, tmax, mass0, beta, mu) )
P = xf[:, 0]
ex = xf[:, 1]
ey = xf[:, 2]
hx = xf[:, 3]
hy = xf[:, 4]
L = xf[:, 5]
cL = np.cos(L)
sL = np.sin(L)
W = 1+ex*cL+ey*sL
Z = hx*sL-hy*cL
C = 1+hx**2+hy**2
q = np.zeros((N+1, 3))
q[:, 0] = P*( (1+hx**2-hy**2)*cL + 2*hx*hy*sL ) / (C*W)
q[:, 1] = P*( (1-hx**2+hy**2)*sL + 2*hx*hy*cL ) / (C*W)
q[:, 2] = 2*P*Z / (C*W)
plt.rcParams['legend.fontsize'] = 20
plt.rcParams['figure.figsize'] = (20, 20)
fig1 = plt.figure()
ax = fig1.gca(projection='3d')
ax.set_xlim3d(-60, 60) # ax.axis('equal') not supported
ax.set_ylim3d(-60, 60)
ax.set_zlim3d(-40, 40)
u, v = np.mgrid[ 0:2*np.pi:20j, 0:np.pi:10j ]
r = 6.378 # Earth radius (in Mm)
x1 = r*np.cos(u)*np.sin(v)
x2 = r*np.sin(u)*np.sin(v)
x3 = r*np.cos(v)
ax.plot_wireframe(x1, x2, x3, color='b')
ax.plot(q[:, 0], q[:, 1], q[:, 2], label='Orbit transfer')
#ax.quiver(q[:, 0], q[:, 1], q[:, 2], u, v, w, length=0.1, normalize=True)
ax.legend()
[5]:
<matplotlib.legend.Legend at 0x7f0ab0adca10>
