Как корректно реализовать обработку нескольких одновременных условий при решении ОДУ scipy_ivp?

Question

Решается задача описания работы двухстороннего пневматического цилиндра с нахождением положения поршней и величины давления в полостях.

При описании работы пневмоцилиндра формируются следующие ограничения:

Если расчётная координата уходит за край системы то обнуляется координата и скорость пневмоцилиндра
```
if X < 0.0:
# Если координата уходит за край системы
X = 0.0  # Координата x = 0
U_X = 0.0  # Скорость  u_x = 0
```
Если рабочий поршень столкнулся с поршнем тормозным, то его координата не может быть больше чем разность максимального хода системы и текущей координаты поршня тормозного
```
if (L_max - Z) < X <= x_max:  # Максимальный ход системы - L_max, поршня - x_max
    X = L_max - z  # z - координата тормозного поршня
```
Поршням ехать за край системы нельзя, но отскакивать от края можно
```
if X > x_max and u_x > 0:
X = x_max
u_x  = 0.0
```
Такие же условия, только для поршня тормозного.

Реализуется решение системы дифференциальных уравнений, описывающих изменение давления в полостях нагнетания перед поршнями и приведение их в движение (с обоюдным взаимодействием)

Система решается с помощью решателя RK45 из стандартного набора решателя solve_ivp, однако решение чрезвычайно длительное, плюс, как видится, неустойчивое в точках разрыва, генерируемых условиями "if-else". Возникло решение провести рефакторинг кода с применением вызываемых событий "events", как точек остановки и последующего "monkey-patching" результатов работы решателя "на лету" с последующим перезапуском решения с новыми начальными условиями. Были сформированы следующие события:

    def work_hit_zero(t, y):
        return y[0]
    work_hit_zero.direction = -1
    work_hit_zero.terminal = True

    def work_hit_brake(t, y):
        return y[0] + y[1] - udp.L_max
    work_hit_brake.terminal = True

    def work_hit_max(t, y):
        return y[0] - udp.x_max
    work_hit_max.terminal = True
    work_hit_max.direction = -1

    def brake_hit_zero(t, y):
        return y[1]
    brake_hit_zero.terminal = True

    def brake_hit_work(t, y):
        return y[0] + y[1] - udp.L_max
    brake_hit_work.terminal = True

    def brake_hit_max(t, y):
        return y[1] - udp.z_max
    brake_hit_max.terminal = True

после формирования системы ОДУ:

   def func(t,u):

   """ODE definition"""

   return d_x, d_z, d_u1, d_u2, d_p1, d_p2, d_p3, d_p4

запускается решатель:

    u = ["Начальные условия"]

    event = [work_hit_zero, work_hit_brake, work_hit_max,
             brake_hit_zero, brake_hit_max, brake_hit_work]

    while True:
        sol = integrate.solve_ivp(func, (t, tend), u, events=event)
        ts.append(sol.t)
        ys.append(sol.y)
        if sol.status == 1:
            t = sol.t[-1]  # старт с новой точки интегрирования
            u = sol.y[:, -1].copy()  # обновление начальной точки
            
            if abs(sol.t_events[5](t,y)) < 1e-12:
                u[1] = z_max
                u[3] = 0.0
                print('Тормозной поршень находится в крайнем положении')

            if abs(events[0](t, y)) < 1e-12:
                u[0] = 0.0  # Координата 0
                u[2] = 0  # Скорость 0.0
                print('Столкновение рабочего поршня со стенкой (началом координат)')

            if abs(events[1](t, y)) < 1e-12:
                u[0] = L_max - u[1]
                print('Столкновение рабочего поршня с тормозным поршнем')
            if abs(events[2](t, y)) < 1e-12:
                u[0] = x_max
                u[2] = 0.0
                print('Рабочий поршень достиг максимальной координаты перемещения')
            if abs(events[3](t, y)) < 1e-12:
                u[0] = 0  # Координата 0
                u[2] = 0  # Скорость 0
                print('Столкновение рабочего поршня со стенкой (началом координат)')
            if abs(events[4](t, y)) < 1e-12:
                u[1] = z_max
                u[3] = 0.0
                print('Тормозной поршень находится в крайнем положении')
            if abs(events[5](t, y)) < 1e-12:
                u[0] = L_max - u[1]
                print('Столкновение рабочего поршня с тормозным поршнем')

        elif sol.status == -1:
            print('Fail!')
            break
        else:
            break

При такой постановке задачи после срабатывания одного из events решатель начинает циклически перезапускать решение, либо после модификации, падает с ошибкой

    if abs(events[0](t, y)) < 1e-12:
    NameError: name 'events' is not defined

Замена events[0] на sol.events[0]:

    KeyError: 'events'

Замена sol.events[0] на sol.y_events[0]:

        if abs(sol.y_events[0](t, y)) < 1e-12:
           NameError: name 'y' is not defined

Замена sol.y_events[0](t, y) на sol.y_events[0](t, u):

        if abs(sol.y_events[0](t, u)) < 1e-12:
           TypeError: 'numpy.ndarray' object is not callable

Что я делаю не так?

Решение такого вида подсмотрел тут.

БЛОГ НА HUSL

Как корректно реализовать обработку нескольких одновременных условий при решении ОДУ scipy_ivp?

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