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

H Jansen h.jansen at fel.tno.nl
Fri Aug 23 10:00:43 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.
>
> 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://projects.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://projects.scipy.org/pipermail/scipy-user/attachments/20020823/fb1d298d/attachment.vcf 


More information about the SciPy-user mailing list