[Solve_ivp积分会在初始条件触发终端为True的事件时卡住

Question

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

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

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

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[1,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
最好同时更新t0和Y0，因此设置

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