I've just done this recently. I't a classical "(dynamic) system identification problem" consiting of a nonlinear least-squares problem with an ode that needs to be integrated repeatedly in the iteration loop. Optimally, you would be able to put bounds on the parameters guiding the solution of the nonlinear optimization process. Unfortunately, scipy doesn't provide such a routine yet for vector systems (for the future a python interface to e.g. Omuses/HQP could provede a solution) so that for the moment we're stuck with unbounded optimization --- this may/will work with not too many parameters in vector p and a good conditioned system. An example is attached. Success. Regards, Henk Jansen provide Nils Wagner wrote:
Hi,
Suppose it is desired to fit a set of data y_i to a known model [given in form of an IVP (initial value problem)]
y' = f(y,p,t) , y(0) = y_0(p), y' = dy/dt
where p is a vector of parameters for the model that need to be found. y denotes the state-variables. The initial conditions y(0) may also depend on parameters.
How can I solve this problem using scipy's optimization and ode tools ?
A small example would be appreciated.
Thanks in advance
Nils _______________________________________________ SciPy-user mailing list SciPy-user@scipy.net http://www.scipy.net/mailman/listinfo/scipy-user
from scipy.optimize import minpack class SystemIdentification: """ SystemIdentification uses the LeastSquares algorithm to find the set of unknown parameters of a dynamic system. Least Squares algorithm: leastsq(func, x0, args=(), Dfun = None, full_output = 0, col_deriv = 0, ftol = 1.49012e-8, xtol = 1.49012e-8, gtol = 0.0, maxfev = 0, epsfcn = 0.0, factor = 100, diag = None) \n\n\t """ __doc__ += minpack.leastsq.__doc__ def __init__(_, residuals = None, p_start = None, args = None, Dfun = None, full_output = 0, col_deriv = 0, ftol = 1.49012e-8, xtol = 1.49012e-8, gtol = 0.0, maxfev = 0, epsfcn = 0.0, factor = 100, diag = None ): _.__residuals = residuals _.__p_start = p_start _.__args = args _.__Dfun = Dfun _.__full_output = full_output _.__col_deriv = col_deriv _.__ftol = ftol _.__xtol = xtol _.__gtol = gtol _.__maxfev = maxfev _.__epsfcn = epsfcn _.__factor = factor _.__diag = diag def initialize (_, residuals, p_start, args = (), Dfun = None ): _.__residuals = residuals _.__p_start = p_start _.__args = args _.__Dfun = Dfun def run(_): if not (_.__residuals and _.__p_start and _.__args): print "*** Error: bad initialization ***" return None return minpack.leastsq(_.__residuals, _.__p_start, _.__args, Dfun = _.__Dfun, full_output = _.__full_output, col_deriv = _.__col_deriv, ftol = _.__ftol, xtol = _.__xtol, gtol = _.__gtol, maxfev = _.__maxfev, epsfcn = _.__epsfcn, factor = _.__factor, diag = _.__diag) if __name__=="__main__": print "Test: System Identification" from test.ode_system import TestOdeSystem full_output = 1 col_deriv = 1 si = SystemIdentification(full_output = full_output, col_deriv = col_deriv) #print si.__doc__ # =================== print "Test 1: Static model (tutorial example)" from scipy import Numeric, RandomArray from Numeric import * from RandomArray import random x = arange (0, 6e-2, 6e-2/30) A, k, theta = 10, 1./3e-2, pi/6 y_true = A*sin(2*pi*k*x+theta) y_meas = y_true + 2.*random(len(x)) #y_meas = y_true def residuals(p, y, x): A,k,theta = p print p err = y-A*sin(2*pi*k*x+theta) return err #p0 = [8, 1/2.3e-2, pi/3] # causes unstable solution process p0 = [12, 35, 0.6] # causes stable solution process print "p0 = ", array(p0) print si.initialize(residuals, p0, args=(y_meas, x)) sol = si.run() print "Solution = ", sol[0] print "Exact = ", array([A, k, theta]) print "Number of function evals = ", sol[1]["nfev"], "(= # iterations)" print print "Rerun again and you'll find the solution unstable (with certain initial conditions)!\n" print "\nNow with Jacobian:\n" def jacobian(p, y, x): A,k,theta = p # err = y-A*sin(2*pi*k*x+theta) z = 2*pi*k*x+theta dfdA = -sin(z) dfdk = -A*cos(z)*2*pi*x dfdtheta = -A*cos(z) jac = zeros((3,len(x)),"d") jac[0] = dfdA # row 1 jac[1] = dfdk # row 2 jac[2] = dfdtheta # row 3 #print jac return jac si.initialize(residuals, p0, args=(y_meas, x), Dfun=jacobian) sol = si.run() print "Solution = ", sol[0] print "Exact = ", array([A, k, theta]) print "Number of function evals = ", sol[1]["nfev"], "(= # iterations)" print "Number of jacobian evals = ", sol[1]["njev"] print print "Rerun again and you'll find THIS solution (using the Jacobian) is STABLE!\n" # print sol # =================== print "Test 2: Dynamic model (paper example)" from scipy import integrate full_output = 1 col_deriv = 0 si = SystemIdentification(full_output = full_output, col_deriv = col_deriv) n = 10 dt = 0.1 #dt = 0.01 k = range(n+1) t = array(k)*dt t_long = concatenate((t,[t[-1]+dt])) # to avoid boundary overflow in InterpolatingFunction class BenchmarkSys: def __init__(_,u, a, x0 = 0. ): _.__u = u _.a = a _.x = array([x0],typecode="f") def f(_,x,t): """ dx/dt = f(x,t) """ return _.a*sin(x)+_.__u bsys = BenchmarkSys(1,-1,0) y_true = integrate.odeint(bsys.f, bsys.x[0], t) y_true = y_true[:,0] scale = 0.25*max(y_true) y_meas = y_true + scale*(random(len(t))-0.5) #y_meas = y_true # print y_meas class OdeSystem: def __init__(_): pass def f(_,x,t): return None def Df_x(_,x,t): return None class UnknownSys(OdeSystem): """ Unknown system. jac - Jacobian system (sensitivity equations; optional) """ def __init__(_,u,a,x0=0., jac = None # jacobian (or sensitivity equations) ): _.u = u _.a = a _.x = array([x0],typecode="f") _.jac = jac _.jac.setSystem(_) OdeSystem.__init__(_) def f(_,x,t): """ dx/dt = f(x,t) """ return _.a*sin(x)+_.u def Df_x(_,x,t): """ Optionally used for time-integration. """ return _.a*cos(x) from Scientific.Functions.Interpolation import InterpolatingFunction class SensitivityEquations(OdeSystem): """ Optional: System of sensitivity equations. Provides the Jacobian for the least-squares problem. """ def __init__(_,x0=0.0): _.x = array([x0],typecode="f") _.sys = None _.z = None # interpolating function OdeSystem.__init__(_) def f(_,x,t): """ dx/dt = f(x,t) """ return sin(_.z(t))+_.sys.a*cos(_.z(t))*x def Df_x(_,x,t): """ Optionally used for time-integration. """ return _.sys.a*cos(_.z(t)) def setSystem(_,sys): _.sys = sys def setInterpFunc(_,t): _.z = InterpolatingFunction((t_long,), sys.x) jac = SensitivityEquations() sys = UnknownSys(1,2,jac=jac) def residuals(p, y, t, sys, integrator): sys.a = p # Optional: Dfun - integration Jacobian # x_sim = integrator(sys.f, sys.x, t) x_sim = integrator(sys.f, sys.x[0], t_long, Dfun=sys.Df_x) sys.x = x_sim[:,0] err = y - sys.x[:-1] return err def jacobian(p, y, t, sys, integrator): """Assumes same signature as residuals()""" # Optional: least-squares Jacobian sys.jac.setInterpFunc(t) x_p_sim = -integrator(sys.jac.f, sys.jac.x[0], t, Dfun=sys.jac.Df_x) sys.jac.x = x_p_sim[:,0] return sys.jac.x with_jacobian = 1 #with_jacobian = 0 if with_jacobian: si.initialize(residuals, sys.a, args=(y_meas, t, sys, integrate.odeint), Dfun=jacobian) else: si.initialize(residuals, sys.a, args=(y_meas, t, sys, integrate.odeint)) sol = si.run() print "Solution = ", sol[0] print "Exact = ", bsys.a print "Number of function evals = ", sol[1]["nfev"], "(= # iterations)" if with_jacobian: print "Number of jacobian evals = ", sol[1]["njev"] #print sol