# [SciPy-user] Parameter estimation / Data fitting in scipy

H Jansen h.jansen at fel.tno.nl
Fri Aug 23 10:04:28 CDT 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.
>
>
>                       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: h.jansen.vcf
Type: text/x-vcard
Size: 471 bytes
Desc: Card for H Jansen
Url : http://projects.scipy.org/pipermail/scipy-user/attachments/20020823/5f0820c9/attachment.vcf
```