如何加速SciPy的odeint？

Question

我调用了一个函数，该函数在每次通过for循环时都使用odeint (遗憾的是，我不能中断该循环中的任何内容)。但是，事情运行得比我希望的要慢得多。代码如下：

def get_STM(t_i, t_f, X_ref_i, dxdt, Amat):
    """Evaluate the state transition matrix rate of change for a given A matrix.
    """

    STM_i = np.eye(X_ref_i.size).flatten()
    args = (dxdt, Amat)
    X_aug_i = np.hstack((X_ref_i, STM_i))
    t = [t_i, t_f]

    # Propogate reference trajectory & STM together!    
    X_aug_f = odeint(dxdt_interface, X_aug_i, t, args=args)
    X_f = X_aug_f[-1, :X_ref_i.size]
    STM_f = X_aug_f[-1, X_ref_i.size:].reshape(X_ref_i.size, X_ref_i.size)

    return X_f, STM_f

def dxdt_interface(X,t,dxdt,Amat):
    """
    Provides an interface between odeint and dxdt
    Parameters :
    ------------
    X : (42-by-1 np array) augmented state (with Phi)
    t : time
    dxdt : (function handle) time derivative of the (6-by-1) state vector
    Amat : (function handle) state-space matrix
    Returns:
    --------
    (42-by-1 np.array) time derivative of the components of the augmented state 
    """
    # State derivative
    Xdot = np.zeros_like(X)
    X_stacked = np.hstack((X[:6], t))
    Xdot_state = dxdt(*(X_stacked))
    Xdot[:6] = Xdot_state[:6].T

    # STM
    Phi = X[6:].reshape((Xdot_state.size, Xdot_state.size))

    # State-Space matrix
    A = Amat(*(X_stacked))
    Xdot[6:] = (A .dot (Phi)).reshape((A.size))

    return Xdot

问题是，我每次调用get_STM大约8640次，这导致了232217次dxdt_interface调用，大约占我总计算时间的70%，每次get_STM调用5ms (其中99.9%是由于odeint)。

我对SciPy的集成技术是个新手，我根本想不出如何基于odeint的documentation来加速这个过程。我研究过用Numba实现dxdt_interface，但我不能让它工作，因为dxdt和Amat是象征性的。

有没有什么技术可以加速我错过的odeint？

编辑:下面包含了Amat和dxdt函数。注意，这些函数不是在我的主要for循环中调用的，它们为传递给我的get_STM函数(我称为import sympy as sym)的符号化函数创建句柄。

def get_A(use_j3=False):
    """ Returns the jacobian of the state time rate of change
    Parameters
    ----------
    R : Earth's equatorial radius (m)
    theta_dot : Earth's rotation rate (rad/s)
    mu : Earth's standard gravitationnal parameter (m^3/s^2)
    j2 : second zonal harmonic coefficient
    j3 : third zonal harmonic coefficient
    Returns
    ----------    
    A : (function handle) jacobian of the state time rate of change
    """
    theta_dot = EARTH['rotation rate']
    R = EARTH['radius']
    mu = EARTH['mu']
    j2 = EARTH['J2']
    if use_j3:
        j3 = EARTH['J3']
    else:
        j3 = 0

    # Symbolic derivations
    x, y, z, mus, j2s, j3s, Rs, t = sym.symbols('x y z mus j2s j3s Rs t', real=True)
    theta_dots = sym.symbols('theta_dots', real=True)
    xdot,ydot,zdot = sym.symbols('xdot ydot zdot ', real=True)

    X = sym.Matrix([x,y,z,xdot,ydot,zdot])

    A_mat = sym.lambdify( (x,y,z,xdot,ydot,zdot,t), dxdt_s().jacobian(X).subs([
        (theta_dots, theta_dot),(Rs, R),(j2s,j2),(j3s,j3),(mus,mu)]), modules='numpy')

    return A_mat

def Dxdt(use_j3=False):
    """ Returns the time derivative of the state vector
    Parameters
    ----------
    R : Earth's equatorial radius (m)
    theta_dot : Earth's rotation rate (rad/s)
    mu : Earth's standard gravitationnal parameter (m^3/s^2)
    j2 : second zonal harmonic coefficient
    j3 : third zonal harmonic coefficient
    Returns
    ----------    
    dxdt : (function handle) time derivative of the state vector
    """

    theta_dot = EARTH['rotation rate']
    R = EARTH['radius']
    mu = EARTH['mu']
    j2 = EARTH['J2']
    if use_j3:
        j3 = EARTH['J3']
    else:
        j3 = 0

    # Symbolic derivations
    x, y, z, mus, j2s, j3s, Rs, t = sym.symbols('x y z mus j2s j3s Rs t', real=True)
    theta_dots = sym.symbols('theta_dots', real=True)
    xdot,ydot,zdot = sym.symbols('xdot ydot zdot ', real=True)

    dxdt = sym.lambdify( (x,y,z,xdot,ydot,zdot,t), dxdt_s().subs([
        (theta_dots, theta_dot),(Rs, R),(j2s,j2),(j3s,j3),(mus,mu)]), modules='numpy')

    return dxdt

Answer 1

在dxdt和Amat作为黑盒的情况下，您不能做很多事情来加速这一过程。一种可能性是简化对它们的调用。hstack可能是大材小用了。

In [355]: def dxdt_quiet(*args):
    x=args
    return x
   .....: 
In [356]: t=1.23
In [357]: dxdt_quiet(*xs)
Out[357]: (0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 1.23)
In [358]: dxdt_quiet(*tuple(x[:6])+(t,))
Out[358]: (0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 1.23)

元组方法要快得多：

In [359]: timeit dxdt_quiet(*tuple(x[:6])+(t,))
100000 loops, best of 3: 5.1 µs per loop
In [360]: %%timeit
xs=np.hstack((x[:6],1.234))
dxdt_quiet(*xs)
   .....: 
10000 loops, best of 3: 25.4 µs per loop

我会做更多这样的测试来优化dxdt_interface调用。