Solve_ivp integration gets stuck if initial condition triggers event with terminal = True
I work in spacecraft trajectory design and I am implementing an algorithm for finding periodic orbits that needs to detect whenever a trajectory crosses a specific plane. I am using scipy.integrate's solve_ivp and its events functionality.
Ideally, I would like to count when events are triggered only stop the integration after a certain number of events are triggered, but since I couldn't find any similar option in solve_ivp's documentation, I am simply using a for cycle to chain these integration arcs and count the number of events. However, my problem is that solve_ivp gets stuck when the initial condition triggers the event and event.terminal = True.
You can find an example code snipet below:
import numpy as np
from scipy.integrate import solve_ivp
def eq_motion(t,Y):
Ydot = np.zeros((6,1))
r = (Y[0,0]**2 + Y[1,0]**2 + Y[2,0]**2)**0.5
Ydot[0] = Y[3,0]
Ydot[1] = Y[4,0]
Ydot[2] = Y[5,0]
Ydot[3] = 2*Y[4,0] + 3*Y[0,0] - Y[0,0]/r**3
Ydot[4] = -2*Y[3,0] - Y[1,0]/r**3
Ydot[5] = -Y[2,0] - Y[2,0]/r**3
return Ydot
def event_xz_plane(t,Y):
return Y[1]
event_xz_plane.terminal = True
event_xz_plane.direction = 0
# initial conditions
Y0 = np.array([-0.006489114696252, 0, 0.107462057974196 + 0.0006, 0, 4.165257687202607,0])
t_span = np.array([0, 1.892845635529040*2])
n_plane_intersections = 4
for i in np.arange(n_plane_intersections):
integrate = solve_ivp(fun=eq_motion, t_span=t_span, y0=Y0, method = 'LSODA',
vectorized = True, rtol = 1e-13, atol = 1e-14,
events = (event_xz_plane),dense_output=True)
print('Number of integrations before event: {}'.format(integrate.y.shape[1]-1))
print(integrate.y_events[0])
print('\n')
# attempting to chain the different integration arcs
Y0 = integrate.y[:,-1]
The initial conditions chosen are close to those of a periodic orbit and produce the following trajectory, where you can see there are 4 intersections with the xz-plane, including the initial condition, which is on that plane.
However, the above script prints out the following:
Number of integrations before event: 1
[[-0.00648911 0. 0.10806206 0. 4.16525769 0. ]]
Number of integrations before event: 1
[[-0.00648911 0. 0.10806206 0. 4.16525769 0. ]]
Number of integrations before event: 1
[[-0.00648911 0. 0.10806206 0. 4.16525769 0. ]]
Number of integrations before event: 1
[[-0.00648911 0. 0.10806206 0. 4.16525769 0. ]]
where you can see that it basically gets stuck in the initial condition. This happens, of course, because the initial condition is situated at a terminal condition, however, there seems to be no option to "pass" the first iteration of the integration.
I tried writing the event function as
def event_xz_plane(t,Y):
return Y[1] + 1e3*(t==0)
so that it would ignore the initial state, but then it simply gets "stuck" on the next plane intersection, and either way it doesn't seem like a good solution.
The only other solution that I thought of was to set terminate=False on the event function and check the number of events triggered and the state vector at these instances. However, this is not ideal for my implementation since it depends on the t_span time set, which is an unknown in the algorithm. That is, if I need to get and n number of plane intersections, I would have to iteratively adjust t_span until I reach that number, instead of just allowing the integration to run until n plane intersections are triggered.
Is there anyway to overcome this?
First check the system of differential equations. How do the linear terms make sense? The third component of the acceleration has its central force part proportional to the third coordinate Y[2]. It is quite possible that with this substantial correction your original approach using a non-terminal event will work.
If you decide to implement the ODE function in the vectorized version, then be fully compliant with the requirement. You can not suppose that the input vector contains only one state, make the output Ydot have the same format as the input Y.
def eq_motion(t,Y):
Ydot = np.zeros(Y.shape)
r = (Y[0]**2 + Y[1]**2 + Y[2]**2)**0.5
Ydot[0] = Y[3]
Ydot[1] = Y[4]
Ydot[2] = Y[5]
Ydot[3] = - Y[0]/r**3 + 2*Y[4] + 3*Y[0]
Ydot[4] = - Y[1]/r**3 - 2*Y[3]
Ydot[5] = - Y[2]/r**3 - Y[2]
return Ydot
It is better to also update t0 along with Y0, so set
t0, tfinal = 0, 10
Use event direction to prohibit the solver to detect the last event again, ideally you would from Y0 and eq_motion(t0,Y0) detect the current sign of the event in positive time direction, here positive, and set the direction opposite to it
event_xz_plane.direction = -1
n_plane_intersections = 4
for i in np.arange(n_plane_intersections):
integrate = solve_ivp(fun=eq_motion, t_span=[t0,tfinal], y0=Y0, method = 'LSODA',
vectorized = True, rtol = 1e-13, atol = 1e-14,
events = (event_xz_plane),dense_output=True)
print(integrate.message)
print('Number of integrations before event: {}'.format(integrate.y.shape[1]-1))
print(integrate.t_events[0],integrate.y_events[0])
print('\n')
# attempting to chain the different integration arcs
Y0 = integrate.y[:,-1]
t0 = integrate.t[-1]
event_xz_plane.direction *= -1
This gives the output
A termination event occurred.
Number of integrations before event: 288
[0.96633212] [[ 3.54891748e-01 -6.67868538e-17 -8.77207053e-01 5.03973235e-02
-7.78052855e-01 -1.76219528e-02]]
A termination event occurred.
Number of integrations before event: 263
[2.03187275] [[ 7.62458860e-03 1.66901228e-15 1.66823579e-01 -2.71527186e-01
3.27865867e+00 1.07443749e-01]]
A termination event occurred.
Number of integrations before event: 300
[3.22281063] [[ 3.86773384e-01 -1.92554306e-16 -8.24376230e-01 -6.86117012e-02
-8.96669320e-01 2.07695448e-01]]
A termination event occurred.
Number of integrations before event: 258
[4.10715034] [[-2.07546473e-02 -9.56266316e-16 1.42957911e-01 9.37459643e-02
3.56344204e+00 7.24687825e-02]]
To compare with, the Dormand-Prince 853 method is much faster (and Radau much slower) and gives the same results
A termination event occurred.
Number of integrations before event: 67
[0.96633212] [[ 3.54891748e-01 -2.46764414e-16 -8.77207053e-01 5.03973235e-02
-7.78052855e-01 -1.76219528e-02]]
A termination event occurred.
Number of integrations before event: 52
[2.03187275] [[ 7.62458861e-03 -5.10008702e-16 1.66823579e-01 -2.71527186e-01
3.27865867e+00 1.07443749e-01]]
A termination event occurred.
Number of integrations before event: 58
[3.22281063] [[ 3.86773384e-01 1.21430643e-17 -8.24376230e-01 -6.86117012e-02
-8.96669320e-01 2.07695448e-01]]
A termination event occurred.
Number of integrations before event: 52
[4.10715034] [[-2.07546473e-02 9.19403442e-17 1.42957911e-01 9.37459642e-02
3.56344204e+00 7.24687825e-02]]