python-3.x - 如果初始条件触发终端 = True 的事件，Solve_ivp 集成将卡住

Question

我从事航天器轨迹设计工作，我正在实现一种算法来寻找周期性轨道，该算法需要在轨迹穿过特定平面时进行检测。我正在使用 scipy.integrate 的 solve_ivp 及其 events 功能。

理想情况下，我想计算事件触发的时间，只在触发一定数量的事件后停止集成，但由于我在 solve_ivp's documentation 中找不到任何类似的选项，我只是使用 for 循环来链接这些积分弧并计算事件数。但是，我的问题是当初始条件触发事件并且 event.terminal = True 时 solve_ivp 卡住了。

您可以在下面找到示例代码片段:

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]

选择的初始条件接近于周期轨道的条件并产生以下 trajectory ，您可以在其中看到与 xz 平面有 4 个交点，包括初始条件，它位于该平面上。

但是，上面的脚本打印出以下内容:

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.        ]]

你可以看到它基本上停留在初始状态。当然，发生这种情况是因为初始条件位于终止条件，但是，似乎没有选项可以“通过”集成的第一次迭代。

我试着把事件函数写成

def event_xz_plane(t,Y):
    return Y[1] + 1e3*(t==0) 

这样它会忽略初始状态，但随后它只会“卡在”下一个平面交叉点上，无论哪种方式，它似乎都不是一个好的解决方案。

我想到的唯一其他解决方案是在事件函数上设置 terminate=False 并检查触发事件的数量和这些实例的状态向量。但是，这对我的实现来说并不理想，因为它取决于 t_span 时间集，这在算法中是未知的。也就是说，如果我需要获取 n 个平面交叉点，我将不得不反复调整 t_span 直到达到该数字，而不是只允许集成运行到 n 平面相交被触发。

有没有办法克服这个问题？

Answer 1

首先检查微分方程组。线性项如何有意义？加速度的第三个分量的中心力部分与第三个坐标 Y[2] 成正比。通过这种实质性的更正，您使用非终端事件的原始方法很可能会奏效。

如果您决定在矢量化版本中实现 ODE 函数，则完全符合要求。你不能假设输入向量只包含一个状态，使输出 Ydot 与输入 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

最好连同Y0一起更新t0，所以设置

t0, tfinal = 0, 10

使用事件方向来禁止求解器再次检测最后一个事件，理想情况下，您可以从 Y0 和 eq_motion(t0,Y0) 检测事件的当前符号正时间方向，这里为正，设置与之相反的方向

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

这给出了输出

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]]

相比之下，Dormand-Prince 853 方法要快得多(而 Radau 要慢得多)并且给出相同的结果

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]]