我使用SciPy的solve_ivp来模拟一个由5个方程组成的模型：

Upsilon = 0.305
Alpha = 0.00000055
r = 0.176
b = 0.0000000000005

def model(t, y0, Beta, Eta, Rmin, p1, p2, p3, A, a, Xi, Epsilon, Lambda, Theta, Mu, Delta, Gamma):

    CD, CT, CM, CE, T, TOTAL = y0

    def k_t(t):
        return Rmin + (p1 / (1 + (p2 * t) ** p3))

    def F_t(T):
        return T / (A + T)

    def f_t(CF, T):
        return (CF / T) / (Upsilon + ((a + CF) / T))

    dCD_dt = -(Beta + Eta) * CD
    dCT_dt = (Eta * CD) + (k_t(t) * F_t(T) * CT) - ((Xi + Epsilon + Lambda) * CT) + (Theta * T * CM) - (Alpha * T * CT)
    dCM_dt = (Epsilon * CT) - (Theta * T * CM) - (Mu * CM)
    dCE_dt = (Lambda * CT) - (Delta * CE)
    dT_dt = ((r * T) * (1 - (b * T))) - (Gamma * f_t((CD + CT), T) * T)

    TOTAL = dCD_dt + dCT_dt + dCM_dt + dCE_dt

    print(t)
    return [dCD_dt, dCT_dt, dCM_dt, dCE_dt, dT_dt, TOTAL]

我正在执行敏感性分析，因此我在输入参数变化的情况下运行了大约37，000次模型。当我运行模型时，我使用以下代码来执行：

#p_values is a numpy array

final = 100
t_range = (0, final)
t_eval = np.arange(0, final, 0.1)

# run model for each parameter combination
for i in range(len(p_values)):
    CD = p_values[i][15]
    CT, CM, CE = 0, 0, 0
    T = 10 ** 7
    TOTAL = CD
    y0 = [CD, CT, CM, CE, T, TOTAL]

    solution = solve_ivp(model, t_range, y0, t_eval=t_eval, args=(p_values[i][0], p_values[i][1], p_values[i][2], p_values[i][3], p_values[i][4], p_values[i][5], p_values[i][6], p_values[i][7], p_values[i][8], p_values[i][9], p_values[i][10], p_values[i][11], p_values[i][12], p_values[i][13], p_values[i][14])).y

但在运行的某个时刻它会...有时它发生得很早（在前一千个解决方案内），有时它可以在卡住之前达到10，000个解决方案或更多。
为了尝试解决这个问题，我开始为模型的每个循环打印t，我注意到一些非常奇怪的事情-当它卡住时，似乎t保持停滞，或者至少增长非常缓慢。这对我来说没有意义，因为我假设它应该总是在每一步增加0.1。每当这种情况发生时，t值都很高（高于95），但我不知道这是否是巧合。如果不是，我想这可能是模型中函数k（t）产生的数字太小的问题，但我不明白为什么这应该是一个问题。
如果有比我更有经验的人有一个建议来尝试和修复它，那显然会非常感激。作为替代，我也可以在运行时间太长时跳过模型的运行并进入下一个运行，但我不知道如何实现这样的东西而不会使运行非常缓慢。