Make propagation functions higher level - poliastro/poliastro GitHub Wiki
Problem
Target: poliastro 0.6
The propagation functions (especially cowell
) are too low level, some unexpected use cases have appeared and their limitations are blocking the development of low-thrust studies.
Explanation
In our implementation of the Cowell method, the raw position and velocity vectors are available in each integration step:
https://github.com/poliastro/poliastro/blob/93d58d/src/poliastro/twobody/propagation.py#L15-L100
However, there are some cases where we need to access other orbital parameters, say inside the acceleration function. For instance, this is an attempt to change the sign of the velocity when the argument of latitude reaches a certain value:
https://github.com/Juanlu001/pfc-uc3m/blob/82f9846/edelbaum.py#L54-L75
Notice how the commented-out debugging code builds the corresponding State
object to access more interesting properties.
Solution
We should evaluate the impact on code and performance of having a higher level object available at each time step, and also provide a way to properly retrieve the state of an object in a sequence of times. The integrator works naturally with times and state vectors, so we could use a decorator to adapt the signatures of the objective function, the callbacks and the perturbation accelerations.
The chosen integration method, Dormand & Prince 8(5,3), "uses a 5th order error estimator and bootstraps on a 3rd order estimator to obtain a dense output of order 7". "Dense output is an interesting concept -- essentially, a dense output of a method provides a solution that can be used as a good linear interpolant for the solution everywhere between output points."
http://math.stackexchange.com/a/941958/24849
Luckily, we can use the set_solout
method to add a callback to retrieve all the intermediate points, instead of sampling:
http://docs.scipy.org/doc/scipy-0.18.0/reference/generated/scipy.integrate.ode.set_solout.html
Keep in mind that this solution is only valid with the cowell
propagator, as kepler
only yields the final state without intermediate calculations. In fact, even if we wanted the Kepler or the naïve method to compute intermediate states too we would be forced to decide a variable step size, probably regularizing through some Sundmann transformation.
On the other hand, this integrator has no prior knowledge of the nature of the problem, so it is going to work for all cases - in a suboptimal way.
Idea: There's the possibility that, regarding propagators and types of orbits, there's no silver bullet.
More ideas
- Implement the naïve propagator used in plotting as another analytical propagator, besides
kepler
- Propagate to a certain value of the (true) anomaly, instead of handling just the time
- Have a high level object that can be queried at any time to retrieve any parameter?
- Export this object as a
DataFrame
with certain values for easier reporting - Use it for plotting
- Export this object as a
Performance issues
The problem is that just instantiating an Orbit
object has a huge performance cost that we cannot afford. Some possible solutions:
- Forget about having high level objects when integrating.
- Use mutable objects instead of creating them in every step.
- Try something like
numba.jitclass
(probably it won't work) or Cython (probably it will be too much work).
⇢ The problem was that astropy.units.quantity_input
is extremely slow and it was called twice. Perhaps we should remove it from *State
classes (hence making them unsafe but noticeably faster). Therefore, we should remove it from all constructors and leave it only in factory methods that are clearly safe.
References
- Ongoing pull request to add IVP solvers with event detection to SciPy http://nbviewer.jupyter.org/gist/nmayorov/ca99381bb20a87762758f7dad4fc2220
- Loosely based on this MATLAB documentation http://www.mathworks.com/help/matlab/math/ode-event-location.html
- In particular, notice how the solver returns an object that can be evaluated in the desired points. "We see that at some regions the solver takes quite large steps, but an accurate solution is available at all points nevertheless." https://github.com/nmayorov/scipy/blob/466a283/scipy/integrate/_py/ivp.py#L149-L176
- Some theory behind event detection http://scicomp.stackexchange.com/a/18907/782
- Alternative approaches:
- SUNDIALS Python wrapper https://github.com/bmcage/odes
- PyDSTool http://www.ni.gsu.edu/~rclewley/PyDSTool/FrontPage.html
- Assimulo http://www.jmodelica.org/assimulo
- Procrastination: Hybrid Systems https://en.wikipedia.org/wiki/Hybrid_system
- More about integration methods https://books.google.es/books?id=F93u7VcSRyYC&printsec=frontcover