如何在solve_ivp中将终端属性用于事件功能

Question

我正在使用Solve_ivp解决大量的一阶ODE。我正在尝试使用事件函数终止求解器，但由于某些原因，它仍会运行到最大时间（1e15）

Solve_Ivp呼叫：

solution = scipy.integrate.solve_ivp(
    lambda t,y: self.rhs(t, y, hasHS),(0, 1e15),simvector,method='BDF', events = lambda t, y: self.eventfun(t, y, hasHS))

这是我的事件功能：

def eventfunWrapper(func):
    @functools.wraps(func)
    def wrapper(*arg,**kwargs):
        func.terminal = True
        func.direction = 0
        return func(*arg,**kwargs)
    return wrapper
## Event function
@eventfunWrapper
def eventfun(self, t, stateVector,hasHS):
    dy = self.rhs(t, stateVector,hasHS)
    x = norm(dy) - 1e-8
    return x

运行代码后，您可以看到找到了多个t_event，但求解器仍在继续运行：

 message: 'The solver successfully reached the end of the integration interval.'
     nfev: 188
     njev: 13
      nlu: 43
      sol: None
   status: 0
  success: True
        t: array([0.00000000e+00, 4.13385988e-05, 8.26771976e-05, 4.96063186e-04,
       9.09449174e-04, 1.54691091e-03, 2.18437264e-03, 2.82183437e-03,
       3.96472459e-03, 5.10761480e-03, 6.25050501e-03, 7.39339522e-03,
       9.55710076e-03, 1.17208063e-02, 1.38845118e-02, 1.60482174e-02,
       2.28891941e-02, 2.97301708e-02, 3.65711476e-02, 5.70008840e-02,
       7.74306204e-02, 1.25358621e-01, 1.73286622e-01, 2.21214623e-01,
       3.79170155e-01, 5.37125688e-01, 7.97053765e-01, 1.05698184e+00,
       1.31690992e+00, 1.79126726e+00, 2.26562460e+00, 2.73998194e+00,
       3.49608272e+00, 4.25218349e+00, 5.00828427e+00, 5.76438504e+00,
       6.89976824e+00, 8.03515144e+00, 9.17053463e+00, 1.50812551e+01,
       2.09919755e+01, 8.00991800e+01, 1.39206384e+02, 7.30278429e+02,
       1.32135047e+03, 7.23207092e+03, 1.31427914e+04, 7.22499958e+04,
       1.31357200e+05, 7.22429245e+05, 1.31350129e+06, 7.22422174e+06,
       1.31349422e+07, 7.22421467e+07, 1.31349351e+08, 7.22421396e+08,
       1.31349344e+09, 7.22421389e+09, 1.31349343e+10, 7.22421388e+10,
       1.31349343e+11, 7.22421388e+11, 1.31349343e+12, 7.22421388e+12,
       1.31349343e+13, 7.22421388e+13, 1.31349343e+14, 7.22421388e+14,
       1.00000000e+15])
 t_events: [array([ 1305948.92714663,  3775274.62982964, 11770897.57749387,
       60058078.65234353, 89829759.09852111])]
        y: array([[  0.  ,   0.  ,   0.  , ...,   0.  ,   0.  ,   0.  ],
       [  0.  ,   0.  ,   0.  , ...,   0.  ,   0.  ,   0.  ],
       [273.15, 273.15, 273.15, ..., 273.15, 273.15, 273.15],
       ...,
       [  0.  ,   0.  ,   0.  , ...,   0.  ,   0.  ,   0.  ],
       [  0.  ,   0.  ,   0.  , ...,   0.  ,   0.  ,   0.  ],
       [  0.  ,   0.  ,   0.  , ...,   0.  ,   0.  ,   0.  ]])

我的猜测是，我没有正确分配终端属性，但不知道如何解决它。

Answer 1

我认为问题出在您的装饰器上。永远不会设置属性。基于此answer，以下代码段应进行设置。

def eventAttr():
    def decorator(func):
        func.direction = 0
        func.terminal = True
        return func
    return decorator


@eventAttr()
def event(self, t, stateVector, hasHS):
    dy = self.rhs(t, stateVector, hasHS)
    x = abs(dy) - 1e-8
    return x


print(event.direction)
print(event.terminal)

按照相同的答案，如果您需要包装，以下代码段将适用

from functools import wraps


def eventAttr():
    def decorator(func):
        @wraps(func)
        def wrapper(*args, **kwargs):
            return func(*args, **kwargs)
        wrapper.direction = 0
        wrapper.terminal = True
        return wrapper
    return decorator


@eventAttr()
def event(self, t, stateVector, hasHS):
    dy = self.rhs(t, stateVector, hasHS)
    x = abs(dy) - 1e-8
    return x


print(event.direction)
print(event.terminal)

全部归功于@Martijn Pieters♦