Оптимизация параметров для системы ODE с использованием scipy.differential_evolution завершается ошибкой из-за того, что resolve_ivp выбрасывает runtime_warnings

Я пытаюсь решить модель на Python и подогнать неизвестные параметры модели к экспериментальным данным. Модель состоит из двух ODE, и я решаю ее с помощью scipy.integrate.solve_ivp. Параметры модели неизвестны, поэтому я хочу подогнать их с помощью множества методов. Мой первый выбор - diff_evolution, так как он может обеспечить действительно хорошие результаты для гораздо более сложных моделей (раньше использовал его как часть другого пакета). Однако проблема заключается в том, что когда я даю свою модель дифференциальной эволюции (я хочу, чтобы она находила глобальный минимум для наименьших квадратов между расчетными и экспериментальными точками), она находит параметры, при которых модель становится нестабильной (LSODA не может интегрировать ее из-за внутреннего шаг равен 0). Я пытался поймать предупреждения среды выполнения, которые выдает LSODA в этих обстоятельствах, но это не помогло. Как лучше всего решить проблему?

Моя модель и мой код ниже:

def aggregation_model(time, z, k1, k2, k3, k_1, k_2):
    GPVI_0 = 55.5
    GPVI, GPVI_Clust = z

    dGPVIdt = - k1 * GPVI_Clust * GPVI + k_1 * (GPVI_0 - GPVI) - 2 * k2 * GPVI * GPVI
    dGPVI_Clustdt = - (k_2 * GPVI_Clust + k3) * GPVI_Clust + k2 * GPVI * GPVI

    return [dGPVIdt, dGPVI_Clustdt]


def res_squares(parameters):

    time = np.linspace(0, 300, 1000)

    timepoints = [0, 25, 50, 75, 100, 125, 150, 300]
    val = [0, 2, 5, 10, 15, 20, 20, 18]

    model_calc = solve_ivp(aggregation_model, [0, 300], [55.5, 0], args=parameters, max_step=100000,
                           dense_output=True, method='LSODA', rtol=1e-12, atol=1e-6)

    solution = model_calc.sol(time)
    transposed = list(map(list, zip(*solution.T)))

    size = [0, ]
    for i in range(1, len(transposed[0])):
        size.append((55.5 - transposed[0][i]) / transposed[1][i])

    diff = []
    for i in range(len(timepoints)):
        indexval = min(range(len(time)), key=lambda j: abs(time[j] - timepoints[i]))
        diff.append((val[i] - size[indexval])**2)

    # print(np.sqrt(np.sum(diff)))

    output = np.sum(diff)

    return output

if __name__ == '__main__':

    from scipy.integrate import solve_ivp
    from scipy.optimize import differential_evolution
    import numpy as np


    parameterBounds = []
    parameterBounds.append([-1000000, 1000000])  # parameter bounds for a
    parameterBounds.append([-1000000, 1000000])  # parameter bounds for a
    parameterBounds.append([-1000000, 1000000])  # parameter bounds for a
    parameterBounds.append([-1000000, 1000000])  # parameter bounds for a
    parameterBounds.append([-1000000, 1000000])  # parameter bounds for a

    result = differential_evolution(res_squares, parameterBounds, seed=3, disp=True, init='random')
    print(result.x)

    sol = solve_ivp(aggregation_model, [0, 300], [55.5, 0], args=parvalues,
                    dense_output=True, method='LSODA', rtol=1e-6, atol=1e-12)