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 не работали так, как предполагалось согласно документации и "как в матлабе". Сейчас необходимо привести её к читабельному виду и повысить производительность. Буду очень благодарен за помощь!


Ответы (0 шт):