from godot import cosmos
# optionally avoid verbose logging messages
import godot.core.util as util
util.suppressLogger()
# create the universe
uniConfig = cosmos.util.load_yaml('universe.yml')
uni = cosmos.Universe(uniConfig)
Generate and Check the Reference Trajectory¶
Now we will load a reference Trajectory file and a spacecraft trajectory file
GODOT cannot currently store the reference trajectory in an IPF, this should be added
# create the trajectory
traConfig = cosmos.util.load_yaml('trajectory_ref.yml')
refTra = cosmos.Trajectory(uni, traConfig)
refTra.compute(False)
When we generate the reference trajectory we normally use the following function to close the orbit. Here we use it to check that the longitude of the ascending nodes is the same at the first and 176th node crossing:
from godot.model import common
import godot.core.tempo as tempo
import godot.core.astro as astro
class LonCon(common.ScalarTimeEvaluable):
def __init__(self, uni, tra):
common.ScalarTimeEvaluable.__init__(self)
self.__uni = uni
self.__tra = tra
def eval(self, epoch):
t1 = self.__tra.point('refASC00001') if type(epoch) is tempo.Epoch else self.__tra.xpoint('refASC00001')
x1 = self.__uni.frames.vector3('Earth', 'SCREF_center', 'ITRF', t1)
lon1 = astro.sphericalFromCart(x1)[1]
t2 = self.__tra.point('refASC00176') if type(epoch) is tempo.Epoch else self.__tra.xpoint('refASC00176')
x2 = self.__uni.frames.vector3('Earth', 'SCREF_center', 'ITRF', t2)
lon2 = astro.sphericalFromCart(x2)[1]
return lon2 - lon1
err = LonCon(uni,refTra)
print ( "The discrepancy in radians is: ", err.eval(tempo.Epoch()))
The discrepancy is close enough.
Controlling The Spacecraft Trajectory¶
Next we can load the spacecraft trajectory
# create the trajectory
traConfig = cosmos.util.load_yaml('trajectory.yml')
tra = cosmos.Trajectory(uni, traConfig)
Now we load all the other modules we need for this computation
import godot.core.autodif as ad
import godot.core.autodif.bridge as br
from godot.model import common, eventgen
import godot.core.tempo as tempo
import godot.core.astro as astro
import godot.core.events as events
The first thing to be done, is to create a ScalarTimeEvaluable to compute the ground track deviation: this is taken from the Napeos problem formulation
class GroundTrackDeviation ( common.ScalarTimeEvaluable ):
def __init__(self, spacecraft:str , reference:str ):
super().__init__()
self.__spacecraft = spacecraft
self.__reference = reference
self.__radiusEarth = uni.constants.getRadius('Earth')
self.__flatteningEarth = uni.constants.getFlattening('Earth')
def __subsatellitePoint(self, e, uni, point):
"""The subsatellite point in ICRF
"""
inertialPos = uni.frames.vector3('Earth', point, 'ITRF', e)
geodeticPos = astro.geodeticFromCartesian(inertialPos, self.__radiusEarth, self.__flatteningEarth)
# FEEDBACK. Now to obtain the sub satellite point we make height = 0.0.
# What happens with the partials when we do this?
geodeticPos[2] = 0.0
res = astro.cartesianFromGeodetic(geodeticPos, self.__radiusEarth, self.__flatteningEarth)
return res
def eval(self, e ):
referenceInertialState = uni.frames.vector6('Earth', self.__reference, 'ITRF', e)
referenceInertialPos = referenceInertialState[0:3]
referenceInertialVel = referenceInertialState[3:6]
referenceSubsat = self.__subsatellitePoint(e, uni, self.__reference)
vert = referenceInertialPos/br.norm(referenceInertialPos)
norm = referenceInertialVel - br.dot(referenceInertialVel, vert)
norm = norm/br.norm(norm)
orto = br.cross(norm, vert)
def funcNorm(t):
pos = uni.frames.vector3('Earth', self.__spacecraft, 'ITRF', t)
return br.dot(norm, pos - referenceSubsat)
grid = [e.value()-100.0,e.value(),e.value()+100]
if (br.isAutodifType(e)):
es = events.generateXEventSet(funcNorm, 1e-6, grid, 1e-9)
else:
es = events.generateEventSet(funcNorm, 1e-6, grid, 1e-9)
subsat = self.__subsatellitePoint(es[0].value(), uni, self.__spacecraft)
return br.dot(orto, subsat - referenceSubsat)
We add an instance of this class to the universe time evaluables:
if uni.evaluables.contains( "groundTrackDeviation" ):
uni.evaluables.remove( "groundTrackDeviation" )
uni.evaluables.add( "groundTrackDeviation", GroundTrackDeviation("SCORB_center", "SCREF_center") )
We want to evaluate the current state of the orbit, so we first compute the trajectory (without partials)
tra.compute(False)
Now we can plot the ground track deviation
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [25, 10]
def plot_vs_time( name, label ):
plt.grid()
plt.xlabel('Time [MJD]')
plt.ylabel(label)
func = uni.evaluables.get(name)
span = tempo.EpochRange(refTra.range().start() + 101.0,refTra.range().end() - 101.0 )
x = [ e.mjd() for e in span.createGrid(360.0) ]
y = [ func.eval( e ) for e in span.createGrid(360.0) ]
plt.plot(x,y)
plot_vs_time("groundTrackDeviation",'GT Dev [km]' )
plt.show()
This is no good for ground track orbit control because it varies strongly within a single orbit. We will sample the function at the ascending nodes and compute a polynomial fit. We will use the following class to convert a Time Evaluable to a Lagrange Polynomial Time Evaluable.
this functionality is being added to godot so the smoothing can be done in c++
We compute the ascending nodes times in a list:
def nodecrossingpoints(span):
# Function with zeros at node crossings
def z(e):
return uni.frames.vector3('Earth', "SCORB_center", 'ITRF', e)[2]
# compute the event interval set
eis = events.generateEventIntervalSet(z, 1e-6, span.createGrid(3600.0), 1e-9 )
# Generate the list of points for the interpolation
ascNode = [ ei.start().value() for ei in eis ]
return ascNode
span = tempo.EpochRange(tra.point('mano1_end') + 300, refTra.range().end() - 101.0 )
ascNode = nodecrossingpoints(span)
points = ascNode[15::30]
points
Now we use the points to compute a smoothed version of the ground track deviation, and its derivative
# import local lagrange
import lagrange
smoothFunc = lagrange.LagrangePolynomialScalarEvaluable( points, uni.evaluables.get("groundTrackDeviation").eval )
if uni.evaluables.contains( "smoothGroundTrackDeviation" ):
uni.evaluables.remove( "smoothGroundTrackDeviation" )
uni.evaluables.add( "smoothGroundTrackDeviation", smoothFunc )
derivSmoothFunc = smoothFunc.deriv()
if uni.evaluables.contains( "smoothGroundTrackDeviationDot" ):
uni.evaluables.remove( "smoothGroundTrackDeviationDot" )
uni.evaluables.add( "smoothGroundTrackDeviationDot", derivSmoothFunc )
Now we plot to show the smoothed function we will use:
plot_vs_time("groundTrackDeviation",'GT Dev [km]' )
plot_vs_time("smoothGroundTrackDeviation",'GT Dev [km]' )
xp = [ e.mjd() for e in points ]
yp = [ smoothFunc.eval( e ) for e in points ]
plt.plot(xp,yp,'o')
xn = [ e.mjd() for e in ascNode ]
yn = [ uni.evaluables.get("groundTrackDeviation").eval( e ) for e in ascNode ]
plt.plot(xn,yn,'.')
plt.show()
We will construct some constraint TimeEvaluable using this smoothed curve. In order to initialise the orbit control process we will introduce a constraint which fixes a point of this curve at a certain value.
Now, we want to optimise a manoeuvre delta v to meet this constraint. So we specify the parameter to track and check it is freed. We will also track the drag scale factor, Cd for analysis purposes.
uni.parameters.track(["mano1_end_dv","drag_cd"])
uni.parameters
Now, we have to re-propagate the trajectories with partial derivatives
tra.compute(True)
refTra.compute(True)
Now we perform the orbit control.
We control a single parameter $p$ which is mano1_end_dv, the control manoeuvre delta-v.
If the constraint is $g_0$, $g_1$ before and after the update then
$$ g_1 = g_0 + \frac { \partial{g} } {\partial{p}} dp $$Therefore the parameter update is:
$$ dp = \frac { g_1 - g_0 } { \frac { \partial{g} } {\partial{p}} } $$def control_deltav(target:float, constraint:common.ScalarTimeEvaluable, xepoch: tempo.XEpoch ):
# Evaluate the constraint at any point
con = common.resolve( constraint.eval( xepoch) )
# Extract the partial you want
dg_dp = con.at("mano1_end_dv")
# Get the initial value of the contraint
g0 = con.value()
# Set the target value
g1 = target
# update the value of the parameter in the parameter book
p = uni.parameters.get("mano1_end_dv").getPhysicalValue()
uni.parameters.get("mano1_end_dv").setPhysicalValue( p + ( g1 - g0 ) / dg_dp )
We use the smoothGroundTrackDeviation as the constraint and set the target value of the constraint to be -30 m, so $g_1= -0.030$ at the end of the arc
there is a bug in the lagrange class, if we interpolate at the points nan is returned, so we so I seconds before.
control_deltav(+0.030, uni.evaluables.get("smoothGroundTrackDeviation"),tempo.XEpoch(points[-1]-1) )
# print out the book to check it is updated
uni.parameters
Now we need to update the trajectory
tra.compute(True)
We check that the procedure was successful by plotting the ground track deviation. At the same time we will plot error bars related to a 10% Cd error.
plot_vs_time("groundTrackDeviation",'GT Dev [km]' )
plot_vs_time("smoothGroundTrackDeviation",'GT Dev [km]' )
xp = [ e.mjd() for e in points ]
yp = [ smoothFunc.eval( e ) for e in points ]
plt.plot(xp,yp,'o')
xn = [ e.mjd() for e in ascNode ]
# get errors from drag
Yn = [ uni.evaluables.get("groundTrackDeviation").eval( tempo.XEpoch( e ) ) for e in ascNode ]
yn = [ y.value() for y in Yn ]
ynErr = [ common.resolve(y).at("drag_cd")[0]* 0.1 for y in Yn ]
plt.errorbar(xn, yn, yerr=ynErr, elinewidth=3)
plt.ylim((-0.15,+0.15))
plt.show()
Now this is correct, the constraint was placed at the end of the cycle.
Next we define a constraint TimeEvaluable we can use to control the minimum of the smooth curve. This takes 2 functions and a time grid. Events are generated from the second function and the time grid and the first of these events is used to evaluate the first function. If no events are found the end point of the search grid is used to evaluate the first function. This is used to define a constraint to be the smoothed ground track deviation, evaluated at the minimum of this function.
this will be implemented in godot optimisation, users will be able to specify a constraint using any functions in the evaluables.
class ConstraintAtEvent(common.ScalarTimeEvaluable):
def __init__(self, uni, func1, func2, grid):
common.ScalarTimeEvaluable.__init__(self)
self.__uni = uni
self.__func1 = func1
self.__func2 = func2
self.__grid = grid
def eval(self, epoch):
g = eventgen.EventGenerator(self.__func2, 1e-6,1e-6)
if (br.isAutodifType(epoch)):
es = g.computeXEvents(self.__grid)
end = tempo.XEpoch(self.__grid[-1])
else:
es = g.computeEvents(self.__grid)
end = self.__grid[-1]
if len(es) == 0:
e = end
else:
func1_values = [ self.__func1.eval(e.value().value()) for e in es ]
ind = func1_values.index(min(func1_values))
e = es[ind].value()
return self.__func1.eval(e)
So we create an instance of our ConstraintAtEvent class, using the smoothGroundTrackDeviation evaluable and its derivative, smoothGroundTrackDeviationDot, to set the set the constraint at the extreme of that function. We perform the grid search for extrema with a step of 1 hour.
grid=span.createGrid(3600.0)
gtdev = ConstraintAtEvent(uni,uni.evaluables.get("smoothGroundTrackDeviation") , uni.evaluables.get("smoothGroundTrackDeviationDot"), grid)
We now need to iterate.
- recompute the ascending nodes and the the smooth lagrange evaluables
- perform the parameter update
- check the result by re-propagating the orbit
# recompute the cross points and create a new smoothed function
ascNode = nodecrossingpoints(span)
points = ascNode[15::30]
smoothFunc = lagrange.LagrangePolynomialScalarEvaluable( points, uni.evaluables.get("groundTrackDeviation").eval )
if uni.evaluables.contains( "smoothGroundTrackDeviation" ):
uni.evaluables.remove( "smoothGroundTrackDeviation" )
uni.evaluables.add( "smoothGroundTrackDeviation", smoothFunc )
derivSmoothFunc = smoothFunc.deriv()
if uni.evaluables.contains( "smoothGroundTrackDeviationDot" ):
uni.evaluables.remove( "smoothGroundTrackDeviationDot" )
uni.evaluables.add( "smoothGroundTrackDeviationDot", derivSmoothFunc )
# perform the parameter update
control_deltav(-0.030, gtdev, tempo.XEpoch())
# compute the trajectory
tra.compute(True)
We can visualise the ground track deviation again:
plot_vs_time("groundTrackDeviation",'GT Dev [km]' )
plot_vs_time("smoothGroundTrackDeviation",'GT Dev [km]' )
xp = [ e.mjd() for e in points ]
yp = [ smoothFunc.eval( e ) for e in points ]
plt.plot(xp,yp,'o')
xn = [ e.mjd() for e in ascNode ]
# get errors from drag
Yn = [ uni.evaluables.get("groundTrackDeviation").eval( tempo.XEpoch( e ) ) for e in ascNode ]
yn = [ y.value() for y in Yn ]
ynErr = [ common.resolve(y).at("drag_cd")[0]* 0.1 for y in Yn ]
plt.errorbar(xn, yn, yerr=ynErr, elinewidth=3)
plt.ylim((-0.15,+0.15))
plt.show()
And then we are done. The parameter value used to generate the manoeuvre commands is:
print( "Delta-v is {0:6.4f} mm/s " .format(1e6*uni.parameters.get("mano1_end_dv").getPhysicalValue()) )
This approach is not limited to ground track control.
Any function can be used in place of GroundTrackDeviation, smoothed if necessary with a suitable choice of sample points. The ConstraintAtEvent class provides functionality to be implemented in GODOT.