Как использовать lmfit python для получения наилучших коэффициентов коэффициентов обыкновенных дифференциальных уравнений (ODE)?

Question

Я пытаюсь минимизировать функцию потерь по сравнению с проблемой обыкновенного дифференциального уравнения (ODE), используя python.

Эта функция потерь предназначена для получения наилучших коэффициентов коэффициента для определяемой пользователем ODE.

Несмотря на мои усилия и примеры, найденные в inte rnet, все результаты заканчивались ODE, чьи решения возвращали значения коэффициентов, которые находились за пределами их соответствующих границ.

Я знаю, что это может произойти по нескольким причинам: среди прочего, выбор минимизатора, неверное определение ODE, размер шага интеграции.

Тем не менее, моя главная проблема - найти пример, где Минимизатор может поддерживать интеграцию «Breaks» в интеграции ODE.

Например, если меня интересует простая модель ODE, например, ODE Prey-Predator (Lotka-Volterra), необходимо установить одно правило: никакие популяции не могут получить отрицательные значения.

Следовательно, если минимизатор ODE по какой-либо заданной причине возвращает отрицательное значение для любой из двух моделируемых совокупностей (Prey и Predator), минимизатор должен прекратить свою интеграцию.

Тем не менее, так как он Похоже, что этот тип правил не поддерживается общими минимизациями, по крайней мере, для объектов Minimizer протестированного мною lmfit.

Даже с ограничивающими правилами, установленными для каждого объекта Parameter, минимизаторы продолжают интерполировать данные за пределами заданных порогов (то есть: отрицательные популяции).

---

Например:

Давайте предположим, что я хочу найти наилучшие коэффициенты коэффициента ODE Prey-Predator по отношению к некоторым эмпирические данные.

В этом случае, если я использую библиотеку lmfit python для своей задачи, я мог бы сделать что-то вроде следующего:

---

from lmfit import Parameters, report_fit, Minimizer

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import ode

end = '\n'*2 + '-'*50 + '\n'*2

def ode_f(t, xs, ps):
    """Lotka-Volterra predator-prey model."""
    if isinstance(ps, Parameters):

        Prey_Natality = ps[r'Prey_Natality'].value
        Predation_rate = ps['Predation_rate'].value
        Predator_Natality = ps['Predator_Natality'].value
        Predator_Mortality = ps['Predator_Mortality'].value
        Prey_natural_mortality = ps['Prey_natural_mortality'].value


    else:
        (Prey_Natality, Predation_rate, Predator_Natality, 
         Predator_Mortality, Prey_natural_mortality) = ps



    Prey, Predator = xs

    dPrey_dt = ( + Prey_Natality*Prey  
                 - Predation_rate*(Prey)*Predator  
                 - Prey_natural_mortality*Prey)

    dPred_dt = ( + Predation_rate* (Prey ) *Predator  
                 + Predator_Natality*Predator  
                 - Predator_Mortality*Predator)

    return [dPrey_dt, dPred_dt]



def solout(y):

    if np.any( y <=0):
        return -1 # stop integration
    else:
        return 0

def ode_solver(t, 
               x0, 
               ps, 
               mode='vode', 
               nsteps=500, 
               method='bdf'):
    """
    Solution to the ODE x'(t) = f(t,x,k) with initial condition x(0) = x0
    """




    r = ode(ode_f).set_integrator(mode, 
                        nsteps=nsteps, 
                        method=method)

    t0 = t.min()
    tmax = t.max()
    dt = np.diff(t)[0]


    r.set_initial_value(x0, t0).set_f_params(ps)

    y = []
    times = []
    integration_time = r.t

    while r.successful() and integration_time < tmax:

        r.integrate(integration_time + dt)


        integration_time = r.t

        yi = r.y

        y.append(yi)

        times.append(integration_time)

        if solout(yi) == -1:
            print('Time stoped at: {0}'.format(integration_time))

            break


    return (np.array(y)).astype(float), times




def ODE_solver_residual_evaluator(ps, ts, data):
    x0 = [ps['Prey_Pop'].value, ps['Predator_Pop'].value]

    model, times = ode_solver(ts, x0, ps)


    # if data.shape[0] <= model.shape[0]:

    #     data.resize((model.shape[0],data.shape[1]),refcheck=False)
    # else:
    #     model.resize((data.shape[0],data.shape[1]),refcheck=False)


    return ( model[:len(data )] - data[:len(model)] )


def residual_dim_reducer(residual_array):

        return np.square(residual_array).sum()

if '__main__' == __name__:
######## 

    dt = 10**(-4)
    t = np.arange(0, 100, dt)

    Prey_initial_pop = 1200
    Predator_initial_pop = 50


    x0 = np.array([Prey_initial_pop,
                   Predator_initial_pop])


    Prey_Natality = 2.6

    Predation_rate = 0.12

    Predator_Natality = 0.401

    Predator_Mortality = 0.0025

    Prey_natural_mortality = 0.001

    true_params = np.array((Prey_Natality, 
                            Predation_rate, 
                            Predator_Natality, 
                            Predator_Mortality,
                            Prey_natural_mortality))

    data, times = ode_solver(t, x0, true_params)

    data += np.random.lognormal(size=data.shape)*0.5


    Param_Names = ['Prey_Natality', 
                   'Predation', 
                   'Predator_Natality', 
                   'Predator_Mortality',
                   'Prey_natural_mortality']

    Populations = ['Prey population',
                   'Predator population']

    for i in range(data.shape[1]):
        plt.plot(times, np.real(data[:,i]), 'o', label=r'original {0}'.format(Populations[i]))

    plt.legend()    
    plt.show()

    import os

    to_save = os.path.join(os.getcwd(), 'original data.png')
    plt.savefig(to_save)



    print('Creating the minizer object', end=end)



    ################3



    # set parameters incluing bounds
    params = Parameters()
    params.add('Prey_Pop', value= 100, min=1, max=600)
    params.add('Predator_Pop', value=10, min=1, max=400)
    params.add('Prey_Natality', value=0.03, min=0.000001, max=3.5)
    params.add('Predation_rate', value=0.02, min=0.00003, max=3.5)
    params.add('Predator_Natality', value=0.0004, min=0.00001, max=3.2)
    params.add('Predator_Mortality', value=0.003, min=0.000001, max=3.2)
    params.add('Prey_natural_mortality', value=0.001, min=0.0000001, max=0.2)



    fitter = Minimizer(ODE_solver_residual_evaluator, 
                       params, 
                       fcn_args=(t, data), 
                       iter_cb=None, 
                       scale_covar=True,
                       nan_policy='propagate', 
                       reduce_fcn=residual_dim_reducer,
                       calc_covar=True, 

                       )

    fitted_ODE = fitter.minimize(method='Nelder-Mead')


    Optimum_params = fitted_ODE.params.valuesdict()


    x0_optimum = np.array([Optimum_params.pop('Prey_Pop'),
                           Optimum_params.pop('Predator_Pop')])

    Y_fitted, times_fitted = ode_solver(t, x0_optimum, Optimum_params.values())


    Param_Names = list(params.keys())


    print(end, 'Param_Names: {0}'.format(Param_Names))


    ##################

    data = data[:len(Y_fitted)]
    Y_fitted = Y_fitted[:len(data)]
    from sklearn.metrics import (explained_variance_score,
                                 r2_score,
                                 mean_squared_error)


    explained_variance_score = explained_variance_score(Y_fitted, data)
    R2 = r2_score(Y_fitted, data)
    RMSE = np.sqrt(mean_squared_error(Y_fitted, data))

    print('Explained variance: {0} \n R2: {1} \n RMSE: {2}'.format(explained_variance_score,
                                                                   R2, RMSE))

    print(end, 'Fitting by Minizer is complete', end=end)

    # display fitted statistics
    report_fit(fitted_ODE)

    print(end)



    ######################3


    fig2 = plt.figure()
    n_ref_markers = 12

    markers_on = np.linspace(0, data[:,i].size-1, n_ref_markers).astype(int).tolist()

    for i in range(Y_fitted.shape[1]):

        plt.plot(times_fitted[:len(times)], 
                 Y_fitted[:len(times),i], 
                 '-', 

                 linewidth=1.1, 
                 label=r'fitted {0} '.format(Param_Names[i]))



        plt.plot(np.arange(len(data)), data[:,i], 'o', 
                 markevery=markers_on,
                 label=r'original {0}'.format(Param_Names[i]))

    fig2.legend()    
    fig2.show()

---

Возвращаемые значения:

1. Две цифры: A) исходные данные B) исходные данные и моделированные данные ODE в соответствии с наилучшим коэффициентом показатели, которые вернул lmfit
2. Таблица возвратов, содержащая подходящие измерения и статистические описания параметров.

---

Возвращенный результат:

--------------------------------------------------

 Fitting by Minizer is complete

--------------------------------------------------
 Models fitting error:

 Explained variance: -90.1682809468072 
 R2: -3125.4358694840084 
 RMSE: 785.9581933129715

--------------------------------------------------

[[Fit Statistics]]
    # fitting method   = Nelder-Mead
    # function evals   = 66437
    # data points      = 364
    # variables        = 7
    chi-square         = 2.2485e+08
    reduced chi-square = 629842.649
    Akaike info crit   = 4867.50583
    Bayesian info crit = 4894.78590
##  Warning: uncertainties could not be estimated:
[[Variables]]
    Prey_Pop:                117.479436 +/- 16.3998313 (13.96%) (init = 100)
    Predator_Pop:            397.552948 +/-        nan (nan%) (init = 10)
    Prey_Natality:           1.05567443 +/-        nan (nan%) (init = 0.03)
    Predation_rate:          3.46543190 +/- 0.05124925 (1.48%) (init = 0.02)
    Predator_Natality:       0.48528830 +/-        nan (nan%) (init = 0.0004)
    Predator_Mortality:      2.76733581 +/- 0.06777831 (2.45%) (init = 0.003)
    Prey_natural_mortality:  0.03928761 +/- 0.00503378 (12.81%) (init = 0.001)
[[Correlations]] (unreported correlations are < 0.100)
    C(Predation_rate, Prey_natural_mortality)     = -0.596
    C(Prey_Pop, Predator_Mortality)               =  0.179
    C(Predation_rate, Predator_Mortality)         =  0.141
    C(Prey_Pop, Prey_natural_mortality)           =  0.127
    C(Predator_Mortality, Prey_natural_mortality) = -0.112

---

---

Эта проблема возникла из примера didacti c, извлеченного из здесь :

Любая помощь по этому вопросу будет принята.

С уважением,