SciPy solve_ivp: почему некоторые события (events) срабатывают, а некоторые - нет?
Описание задачи
Пишу программу по баллистическому расчёту ракеты. Каждый участок траектории описывается собственной системой дифференциальных уравнений. Для решения каждой системы использую scipy.solve_ivp (версия 1.4.1). Интервал времени для решения каждой системы задаю заведомо большим, чтобы наверняка достичь необходимых параметров движения.
Проблема
Чтобы не производить лишние расчёты и затем не "чистить" массивы рассчитанных данных (т.е. не тратить время на поиск и удаление бесполезных данных, типа отрицательной высоты), необходимо прерывать решение систем ДУ.
И здесь начались странности, связанные с events.
В документации на solve_ivp https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html написано, что на вход в качестве events можно подавать как callable, так и список callables (list of callables). Т.к. движение ракеты довольно сложное, то зачастую приходится одновременно использовать несколько "событий" для прерывания решения.
Но заметил (используя debug в виде print()), что для __move_const_theta() всё работает, как надо. А для остальных функций ___move... работает только событие __event_hit_ground(). Объяснения как раз этому феномену я дать не могу. Причём не имеет значения, сколько событий передаётся в events - работает только __event_hit_ground().
Повторюсь: для __move_const_theta() всё работает, как должно.
Код
Основная функция класса ракеты
def start(self, max_step=0.2):
"""
Launches the rocket.
:param max_step: max step for the ODE's solver
"""
# Fly
self.engine.on(False)
t_span = np.array([0.0, self.optim_params.k_t_start * self.traj.t_start_end])
self.__move_const_theta(t_span, max_step / 3)
# Check engine
if self.state.mu > self.engine.mu_break[self.engine.mode]:
self.engine.off()
self.state.Theta = self.optim_params.Theta_start
if self.engine.is_run is True:
t_span = np.array([self.state.t, self.state.t + 500.0])
events = [self.__event_t_start, self.__event_mu_start]
self.__move_const_theta(t_span, max_step / 2, events)
self.engine.off()
t_span = np.array([self.state.t, self.state.t + 500.0])
events = [self.__event_max_height, self.__event_march_height, self.__event_hit_ground]
self.__move_passive(t_span, max_step / 2, events)
self.state.Theta = 0.0
if self.engine.on(True) is True:
t_span = np.array([self.state.t, self.state.t + 100.0])
events = [self.__event_mu_max]
self.__move_const_theta(t_span, max_step, events)
self.engine.off()
t_span = np.array([self.state.t, self.state.t + 1000.0])
events = [self.event_hit_ground]
self.__move_passive(t_span, max_step / 4, events)
Функции расчёта высокого уровня
def __move_const_theta(self, t_span: np.ndarray, max_step: float, events=None):
"""
:type events: list
"""
y0 = self.state.to_ndarray_xy()
if events is None:
sol = solve_ivp(self.__ode_const_theta, t_span, y0, max_step=max_step)
else:
sol = solve_ivp(self.__ode_const_theta, t_span, y0, max_step=max_step, events=events)
self.update(sol.t, sol.y)
def __move_passive(self, t_span: np.ndarray, max_step: float, events=None):
"""
:type events: list
"""
y0 = self.state.to_ndarray_xy()
if events is None:
sol = solve_ivp(self.__ode_passive, t_span, y0, max_step=max_step)
else:
sol = solve_ivp(self.__ode_passive, t_span, y0, max_step=max_step, events=events)
self.update(sol.t, sol.y)
def __move_maneuver(self, t_span: np.ndarray, max_step: float, events=None):
"""
:type events: list
"""
y0 = self.state.to_ndarray_xy()
if events is None:
if self.engine.is_run is True:
sol = solve_ivp(self.__ode_active_man, t_span, y0, max_step=max_step)
else:
sol = solve_ivp(self.__ode_passive_man, t_span, y0, max_step=max_step)
else:
if self.engine.is_run is True:
sol = solve_ivp(self.__ode_active_man, t_span, y0, max_step=max_step, events=events)
else:
sol = solve_ivp(self.__ode_passive_man, t_span, y0, max_step=max_step, events=events)
self.update(sol.t, sol.y)
Думаю, низкоуровневые функции типа __ode... не представляют интереса - в них всего лишь записана система дифуров.
Функции событий
def __event_t_start(self, t, y):
return t - self.traj.t_start_end
__event_t_start.direction = 1
__event_t_start.terminal = True
def __event_theta_reached(self, t, y):
print(y[4] - self.maneuver.Theta_end)
return y[4] - self.maneuver.Theta_end
__event_theta_reached.direction = 0
__event_theta_reached.terminate = True
def __event_max_height(self, t, y):
return y[2] - self.constraints.max_height
__event_max_height.direction = 1
__event_max_height.terminate = True
def __event_march_height(self, t, y):
return self.optim_params.h_march - y[2]
__event_march_height.direction = 1
__event_march_height.terminate = True
def __event_hit_ground(self, t, y):
return y[2] - 500.0 # accuracy
__event_hit_ground.direction = -1
__event_hit_ground.terminal = True
def __event_max_dist(self, t, y):
return y[1] - self.constraints.max_dist
__event_max_dist.direction = 1
__event_max_dist.terminal = True
def __event_mu_start(self, t, y):
return y[3] - self.engine.mu_break[0]
__event_mu_start.direction = 1
__event_mu_start.terminal = True
def __event_mu_max(self, t, y):
return y[3] - self.engine.mu_break[-1]
__event_mu_max.direction = 1
__event_mu_max.terminal = True
P.S.
Остальной код проверен и проблем не вызывает. Раньше данная программа задачу решала, но была напичкана костылями именно по причине того, что events не работали так, как предполагалось согласно документации и "как в матлабе". Сейчас необходимо привести её к читабельному виду и повысить производительность. Буду очень благодарен за помощь!