fitting parameters of ODEs while using octave/matlab ODE solver

Question

I am using OdePkg in Octave to solve a system of stiff ODEs, e.g. by ode5r:

function yprime = myODEs(t,Y,param)
    yprime = [
        - param(1) * Y(1);                      # ODE for Y(1)
        param(1) * Y(1) - param(2) Y(2) * Y(3); # ODE for Y(2)
        param(2) Y(2) * Y(3)                    # ODE for Y(3)
                                                # etc.
];

time_span = [1, 24]         # time span of interest
Y0        = [1.0, 1.1, 1.3] # initial values for dependent variables Y
param     = [7.2, 8.6, 9.5] # parameters, to be optimized

[t, Y] = ode5r(@myODEs, time_span, Y0, ..., param);

The solver stores the dependent variables Y in a matrix with respect to time t (vector):

t     Y(1)  Y(2)  Y(3)
0.0   1.0   1.1   1.3
0.1   ...   ...   ...
0.5   ...   ...   ...
0.9   ...   ...   ...
...   ...   ...   ...
4.0   ...   ...   ...
...   ...   ...   ...
24.0  ...   ...   ...

I want to fit the parameters in param, so that the resulting variables Y best fit my reference values, e.g.:

t         Y(1)  Y(2)  Y(3)
0.5       1.1   N/A   N/A
1.0       1.9   N/A   N/A
4.0       2.3   2.7   2.1
5.0       N/A   2.6   2.2
24.0      0.9   1.5   2.0

Which Octave/Matlab (other languages are welcome) routine can perform a multi-parameter (least square / spline) fit? How is it possible to combine parameter sets for different initial values Y0 in the fit? I would be glad if you could provide me with some hints and possibilities.

Best regards, Simon

Answer 1

Do you mean that each function y(t) needs to be fitted? In that case, a lease squares or spline fitting for each set of Yi versus time will work just fine. Can't tell which one would be best without seeing your data.

You'll have to come up with another independent variable if you mean that you want to fit a curve across all values of Yi for a given time point and then watch that curve evolve over time.

UPDATE: Least square fitting is what it is - I don't have a particular routine to recommend. SciPy has one, I'm sure. I'm sorry that I don't have a better recommendation. I'm only learning Python now.

I don't know what you mean by "fitness indicator". Least squares fitting calculates coefficients that minimize the mean square of the error between the fit and the data at each point.

Just one way to combine several data sets into a single fit: merge them and re-run the calculation.

Answer 2

This should be relatively straightforward with scipy. scipy.optimize.leastsq() takes a function that should return an array of residuals for a given parameter vector. It will minimize the sum of the squares of the residuals. To handle multiple datasets with different initial values, you just run the ODE once for each dataset, compute the residuals for each dataset/run pair, and then concatenate the residual vectors together. Here is a rough sketch:

import numpy
from scipy import integrate, optimize

# The initial guess.
p0 = numpy.array([7.2, 8.6, 9.5])

# The collected datasets.
# A list of (t, y0, y) tuples.
# The y's are all (len(y0), len(t))-shaped arrays. The output of
# integrate.odeint is also in this format.
datasets = [...]

def odes(y, t, params):
    dydt = [
        -params[0] * y[0],
        params[0]*y[0] - params[1]*y[1]*y[2],
        params[1]*y[1]*y[2],
    ]
    return np.array(dydt)

def residuals(params, datasets):
    res = []
    for t, y0, y in datasets:
        res.append(integrate.odeint(odes, y0, t, args=(params,)) - y)

    # Stack them horizontally and flatten the array into the expected vector.
    # You're on your own for handling missing data. Look into the numpy.ma
    # module.
    all_residuals = numpy.hstack(res).ravel()
    return all_residuals

opt_params, err = optimize.leastsq(residuals, p0, args=(datasets,))