XEvents

XEvents and optimisation problems¶

Note this doesn't work all the way at the moment.

Using autodif in many algorithms allows easy computation of partial derivatives, root finding is an example. If Event times are computed using autodif, the result is called an XEvent, which is an event containing an XEpoch including gradients.

In this example we are going to use XEvents to compute gradients for constraints.

When we specify a constraint, we define a point at which that constraint applies, this can be a fixed Epoch, or an Event, e.g. a periapsis passage or an eclipse entry.

Generally speaking Events are defined by an event function $g\left(t, p\right)$ generating an event set $E = \{t_g\} $ by root finding.

For each event in $E$, by definition $g\left(t_g, p\right)=0$.
If an optimisation results in an updated value of the parameters, $p'$, then $g\left(t_g, p' \right) \ne 0 $
Solving for the roots $E' = \{t'_g\}$, we again have $g\left(t'_g, p' \right) = 0 $

The gradient of $g\left(t_g, p\right)$, taking into account the dependency of the function on the event time, is expected to be zero:

$$ \frac{ \partial g}{\partial p}\Bigr|_{g=0} = \frac{ \partial g}{\partial p} + \frac{ \partial g}{\partial t_g}\frac{ \partial t_g}{\partial p} = 0 $$

Suppose we compute another function $f\left(t, p\right)$ at the XEpoch of a solved event from $E$, using $g$. The gradients of $f$ then take into account the gradient of the event time to the parameters.

$$ \frac{ \partial f}{\partial p}\Bigr|_{g=0} = \frac{ \partial f}{\partial p} + \frac{ \partial f}{\partial t_g}\frac{ \partial t_g}{\partial p} $$

This is the gradient of $f$, accounting for the condition that it applies at a root of $g$, so the dependency of $f$ on the roots of $g$ is included.

If we omit the extra term, the gradient will be wrong and the constraint will not be satisfied after an update. This can result in convergence problems, or a failure of the optimisation to converge.

Now we will demonstrate this process.

In [16]:

from godot.core import tempo, events
from godot.model import common
from godot import cosmos
import ruamel.yaml as yaml

# optionally avoid verbose logging messages
import godot.core.util as util
util.suppressLogger()

Load the universe and trajectory objects:

In [17]:

uni_file = 'universe.yml'
with open(uni_file) as unifile:
    uni_config = yaml.safe_load(unifile)

uni = cosmos.Universe(uni_config)

tra_file = 'trajectory.yml'
with open(tra_file) as trafile:
    tra_config = yaml.safe_load(trafile)

tra = cosmos.Trajectory(uni, tra_config)

Compute the trajectory, tracking the gradients of the spacecraft pericenter and apocenter radius.

In [18]:

uni.parameters.track(['launch_GeoSat_center_rpe', 'launch_GeoSat_center_rap'])
tra.compute(True)

Define an event function, the distance of the spacecraft from Earth, minus 20k km.

In [19]:

def g(e):
    d = uni.frames.distance('Earth', 'GeoSat_center', e)
    return d - 20000.0

Compute the events in the range of the trajectory and print their values relative to the start of the arc, resolving the gradients:

In [20]:

tol = 1e-9
eps = 1e-3
grid = tra.range().createGrid(3600.0)
grid[0]+=1
grid[-1]-=1
eventset = events.generateXEventSet(g, eps, grid, tol)

print("The event times relative to the start of the propagation, including gradients:")
for evt in eventset:
    print( common.resolve( evt.value() - tra.range().start() ) )

The event times relative to the start of the propagation, including gradients:
         Value | launch_GeoSat_center_rpe | launch_GeoSat_center_rap
  3.266717e+03 |             9.918227e-02 |            -5.521718e-03

Next we compute the value of the event function at the solved event times. We do this using the event gradients in the event times, and then removing the event gradient in the event times.

In each case we resolve the gradient partial derivatives. This means applying the chain rule to compute the partial derivatives with respect to the original parameters, the pericenter and apocenter radii.

In [21]:

print("The gradients of g using the gradients of the event times:\n")
for evt in eventset:
    xepochWithPartials = evt.value()
    print(common.resolve(g( xepochWithPartials )))

print("The gradients of g NOT using the gradients of the event times:\n")
for evt in eventset:
    xepochWithoutPartials = tempo.XEpoch(evt.value().value()) # create an xepoch with no partials
    print(common.resolve(g(xepochWithoutPartials)))

The gradients of g using the gradients of the event times:

         Value | launch_GeoSat_center_rpe | launch_GeoSat_center_rap
  6.839400e-10 |            -2.220446e-16 |            -3.469447e-18

The gradients of g NOT using the gradients of the event times:

         Value | launch_GeoSat_center_rpe | launch_GeoSat_center_rap
  6.839400e-10 |            -4.245562e-01 |             2.363608e-02

Notice that the gradient of the event function with respect to the free parameters at the solved event time is (close to) zero when the event gradient are used, in the computation and non zero when the event gradient is not used.

Now we introduce another function, $f$ the radial component of the velocity

In [22]:

import godot.core.autodif.bridge as br
def f( epo ):
    r = uni.frames.vector6('Earth', 'GeoSat_center', 'ICRF', epo)
    return br.dot(r[:3],r[3:])/br.norm(r[:3])

We compute the function $f$ and resolved gradient with and without the event gradient.

In [23]:

print("The gradients of f using the gradients of the event times:\n")
for evt in eventset:
    xepochWithPartials = evt.value()
    print(common.resolve(f(xepochWithPartials)))
print("The gradients of f NOT using the gradients of the event times:\n")
for evt in eventset:
    xepochWithoutPartials = tempo.XEpoch(evt.value().value())
    print(common.resolve(f(xepochWithoutPartials)))

The gradients of f using the gradients of the event times:

         Value | launch_GeoSat_center_rpe | launch_GeoSat_center_rap
  4.280566e+00 |            -1.828480e-04 |             1.101240e-05

The gradients of f NOT using the gradients of the event times:

         Value | launch_GeoSat_center_rpe | launch_GeoSat_center_rap
  4.280566e+00 |            -1.449101e-04 |             8.900303e-06

Notice the gradients are different. This is because the first result computes the gradient of the function $f$ at the event time using $g(t)=0$.

Now we are going to compute the update required in launch_GeoSat_center_rpe to set $f$ equal to a target value. We compute $g$ and $f$ at the a priori and the updated event times. First without using the event time partials:

In [26]:

useEventPartials = False
par = 'launch_GeoSat_center_rpe'

uni_config = yaml.safe_load(open('universe.yml'))
tra_config = yaml.safe_load(open('trajectory.yml'))
uni = cosmos.Universe(uni_config)
tra = cosmos.Trajectory(uni, tra_config)

def f( e ):
    r = uni.frames.vector6('Earth', 'GeoSat_center', 'ICRF', e)
    return br.dot(r[:3],r[3:])/br.norm(r[:3])

f_t = -5.17

def g(e):
    d = uni.frames.distance('Earth', 'GeoSat_center', e)
    return d - 20000.0

uni.parameters.track([par])
tra.compute(True)

# compute the a priori values of g and f, and their gradients
grid = tra.range().createGrid(3600.0)
#grid[-1]-=1
eventset = events.generateXEventSet(g, 1e-3, grid[1:-1], 1e-9)
print(eventset)

xepoch = eventset[0].value()
if useEventPartials:
    xepoch = eventset[0].value()
else:
    xepoch = tempo.XEpoch(eventset[0].value().value())    

print('A priori, g, f')
print(common.resolve(g(xepoch)))
print(common.resolve(f(xepoch)))

# estimate the linear correction to p to set f to the target value
f0 = common.resolve(f(xepoch))
dfdp = f0.at(par)
dp = (f_t - f0)/dfdp
p = uni.parameters.get(par)
uni.parameters.get(par).setPhysicalValue( p.getPhysicalValue() + dp.value() )

# recompute
tra.compute(True)
grid = tra.range().createGrid(3600.0)
grid[-1]-=1
eventset = events.generateXEventSet(g, 1e-3, grid[1:-1], 1e-9)
if useEventPartials:
    xepoch = eventset[0].value()
else:
    xepoch = tempo.XEpoch(eventset[0].value().value())

print('Updated g, f')
print(common.resolve(g(xepoch)))
print(common.resolve(f(xepoch)))

And then using the event time partials:

In [25]:

useEventPartials = True
par = 'launch_GeoSat_center_rpe'

uni_config = yaml.safe_load(open('universe.yml'))
tra_config = yaml.safe_load(open('trajectory.yml'))
uni = cosmos.Universe(uni_config)
tra = cosmos.Trajectory(uni, tra_config)

def f( e ):
    r = uni.frames.vector6('Earth', 'GeoSat_center', 'ICRF', e)
    return br.dot(r[:3],r[3:])/br.norm(r[:3])

f_t = -5.17

def g(e):
    d = uni.frames.distance('Earth', 'GeoSat_center', e)
    return d - 20000.0

uni.parameters.track([par])
tra.compute(True)

# compute the a priori values of g and f, and their gradients
grid = tra.range().createGrid(3600.0)
grid[-1]-=1
eventset = events.generateXEventSet(g, 1e-3, grid[1:-1], 1e-9)

xepoch = eventset[0].value()
if useEventPartials:
    xepoch = eventset[0].value()
else:
    xepoch = tempo.XEpoch(eventset[0].value().value())    

print('A priori, g, f')
print(common.resolve(g(xepoch)))
print(common.resolve(f(xepoch)))

# estimate the linear correction to p to set f to the target value
f0 = common.resolve(f(xepoch))
dfdp = f0.at(par)
dp = (f_t - f0)/dfdp
p = uni.parameters.get(par)
uni.parameters.get(par).setPhysicalValue( p.getPhysicalValue() + dp.value() )

# recompute
tra.compute(True)
grid = tra.range().createGrid(3600.0)
grid[-1]-=1
eventset = events.generateXEventSet(g, 1e-3, grid[1:-1], 1e-9)
if useEventPartials:
    xepoch = eventset[0].value()
else:
    xepoch = tempo.XEpoch(eventset[0].value().value())

print('Updated g, f')
print(common.resolve(g(xepoch)))
print(common.resolve(f(xepoch)))

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-25-56671b984e08> in <module>
     25 eventset = events.generateXEventSet(g, 1e-3, grid[1:-1], 1e-9)
     26 
---> 27 xepoch = eventset[0].value()
     28 if useEventPartials:
     29     xepoch = eventset[0].value()

RuntimeError: 
Exception raised in GODOT: 
What : Pre-condition check failed: index < size. Help message: Element access outside of maximum bound
Func : transformIndex
Src  : /tmp/pip-req-build-xmk07fz9/godotpy/common/Slice.h
Line : 12

In [ ]:

import sys
sys.path

Out[ ]:

['/home/mas/emcd/mcd53/mcd/_build/lib/python',
 '/usr/lib64/python36.zip',
 '/usr/lib64/python3.6',
 '/usr/lib64/python3.6/lib-dynload',
 '',
 '/local/scratch/gvarga/godottut/env/lib64/python3.6/site-packages',
 '/local/scratch/gvarga/godottut/env/lib/python3.6/site-packages',
 '/local/scratch/gvarga/godottut/env/lib64/python3.6/site-packages/IPython/extensions',
 '/home/gvarga/.ipython']

The event time is updated after the iteration, and the constraint function $f$ is closer to to the target value in the second case, when its gradient includes contributions to account for changes in the event time due to changes in the parameters.

This is how constraint and cost functions (optimisation functions) are computed during an optimisation. The optimisation function is computed at the xepoch obtained from a zero of the event function which defines its location.

Optimisation functions can be added after the trajectory has been computed, this makes the configuration of the optimisation problem simpler, since only aspects related to the motion need to be specified and not aspects related to the optimisation problem.