Подгонка данных к системе ODE с использованием Python через Scipy & Numpy

У меня возникли проблемы с переводом моего кода MATLAB на Python через Scipy & Numpy. Я застрял в том, как найти оптимальные значения параметров (k0 и k1 )для моей системы ОДУ, чтобы они соответствовали моим десяти наблюдаемым точкам данных. В настоящее время у меня есть начальное предположение для k0 и k1. В MATLAB я могу использовать что-то под названием «fminsearch», которая представляет собой функцию, которая принимает систему ОДУ, наблюдаемые точки данных и начальные значения системы ОДУ. Затем он рассчитает новую пару параметров k0 и k1, которые будут соответствовать наблюдаемым данным. Я включил свой код, чтобы узнать, можете ли вы помочь мне реализовать какой-то «fminsearch», чтобы найти оптимальные значения параметров k0 и k1, которые будут соответствовать моим данным. Я хочу добавить любой код для этого в мой файл lsqtest.py.

У меня есть три файла.py -ode.py, lsq.py и lsqtest.py

. ode.py:

def f(y, t, k): 
return (-k[0]*y[0],
      k[0]*y[0]-k[1]*y[1],
      k[1]*y[1])

лск.py:

import pylab as py
import numpy as np
from scipy import integrate
from scipy import optimize
import ode
def lsq(teta,y0,data):
    #INPUT teta, the unknowns k0,k1
    # data, observed 
    # y0 initial values needed by the ODE
    #OUTPUT lsq value

    t = np.linspace(0,9,10)
    y_obs = data #data points
    k = [0,0]
    k[0] = teta[0]
    k[1] = teta[1]

    #call the ODE solver to get the states:
    r = integrate.odeint(ode.f,y0,t,args=(k,))


    #the ODE system in ode.py
    #at each row (time point), y_cal has
    #the values of the components [A,B,C]
    y_cal = r[:,1] #separate the measured B
    #compute the expression to be minimized:
    return sum((y_obs-y_cal)**2)

lsqtest.py:

import pylab as py
import numpy as np
from scipy import integrate
from scipy import optimize
import lsq


if __name__ == '__main__':

    teta = [0.2,0.3] #guess for parameter values k0 and k1
    y0 = [1,0,0] #initial conditions for system
    y = [0.000,0.416,0.489,0.595,0.506,0.493,0.458,0.394,0.335,0.309] #observed data points
    data = y
    resid = lsq.lsq(teta,y0,data)
    print resid
5
задан tacaswell 19 December 2012 в 15:03
поделиться