[SciPy-user] Parameter estimation / Data fitting in scipy
H Jansen
h.jansen at fel.tno.nl
Fri Aug 23 11:00:43 EDT 2002
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 at scipy.net
> http://www.scipy.net/mailman/listinfo/scipy-user
-------------- next part --------------
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
-------------- next part --------------
A non-text attachment was scrubbed...
Name: sysident.ps
Type: application/postscript
Size: 115192 bytes
Desc: not available
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20020823/fb1d298d/attachment.ps>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: h.jansen.vcf
Type: text/x-vcard
Size: 471 bytes
Desc: Card for H Jansen
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20020823/fb1d298d/attachment.vcf>
More information about the SciPy-User
mailing list