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
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
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
------------------- Thank you for your prompt reply. Please can you send me the missing module
python hansen.py
Test: System Identification Traceback (most recent call last): File "hansen.py", line 96, in ? from test.ode_system import TestOdeSystem ImportError: No module named ode_system Cheers, Nils
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
Sorry. Although it's used, it doesn't add much functionality though. Regards, Henk Jansen My VDI Freemail wrote:
-------------------
Thank you for your prompt reply. Please can you send me the missing module
python hansen.py
Test: System Identification Traceback (most recent call last): File "hansen.py", line 96, in ? from test.ode_system import TestOdeSystem ImportError: No module named ode_system
Cheers,
Nils
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
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.net http://www.scipy.net/mailman/listinfo/scipy-user
from scipy.integrate import odeint class OdeSystem: def __init__(_): pass def f(_): pass def y(_): pass from math import * class TestOdeSystem (OdeSystem): def __init__(_): _._x = 0. _._u = 0. _._theta = 0. # OdeSystem.__init__(_) def f(_): """dxdt = f(x,u; theta)""" return _._theta*sin(_._x)+u def y(_): """y=x""" return _.__x if __name__=="__main__": print "Test: OdeSystem"
SciPy testing currently runs in an infinite loop for me. Not sure why. import scipy scipy.test() just keeps running repeating how it's building tests for the different files. Is this happening to anyone else? -Travis O.
On Tue, 17 Sep 2002, Travis Oliphant wrote:
SciPy testing currently runs in an infinite loop for me. Not sure why.
import scipy scipy.test()
just keeps running repeating how it's building tests for the different files.
Is this happening to anyone else?
Same here. It appears that weave tests cause this loop (I haven't check the details though) as scipy.linalg.test(1) scipy.stats.test(1) scipy.optimize.test(1) etc run fine but scipy.weave.test(1) runs into infinite loop. Pearu
Travis> SciPy testing currently runs in an infinite loop for me. Not Travis> sure why. ... Travis> Is this happening to anyone else? I've been wondering why the build on the Sun took so looooong. All this time I thought it was the slow machine or lack of memory. Based on Pearu's response about weave problems, I wonder if something I checked in hosed things. (I did change the directory structure so that multiple users on the same machine can build weave modules.) Skip
-----Original Message----- From: scipy-user-admin@scipy.net [mailto:scipy-user-admin@scipy.net] On Behalf Of Skip Montanaro Sent: Tuesday, September 17, 2002 1:07 PM To: scipy-user@scipy.net Subject: Re: [SciPy-user] Scipy Testing in an infinite loop
Travis> SciPy testing currently runs in an infinite loop for me. Not Travis> sure why. ... Travis> Is this happening to anyone else?
I've been wondering why the build on the Sun took so looooong. All
time I thought it was the slow machine or lack of memory.
Based on Pearu's response about weave problems, I wonder if something I checked in hosed things. (I did change the directory structure so
Not sure what the issue is. I did make a change last night to weave concerning the location of testing.py, but it should cause failures, not infinite loops do to an import error (I haven't moved the actual file yet.) I remember seeing this a looong time ago. I think it had something to do with the test harness. Seems like it "just started working" after tinkering though, and I never found the real problem. I'll spend some time on that this afternoon and try to chase it down. Also, weave has had many changes in the last few days so that could also be of issue. On testing.py's location. After further consideration, I think it should live in its own package as Travis O. suggests. Pearu is right that much of the scipy_distutils stuff could be folded back into distutils and so that package might actually disappear. I think its total disappearance is unlikely, but scipy_testing still makes prtty good sense. This is going to touch *a lot* of code. I will be making the change today. eric this that
multiple users on the same machine can build weave modules.)
Skip _______________________________________________ SciPy-user mailing list SciPy-user@scipy.net http://www.scipy.net/mailman/listinfo/scipy-user
participants (7)
-
eric jones -
H Jansen -
My VDI Freemail -
Nils Wagner -
Pearu Peterson -
Skip Montanaro -
Travis Oliphant