Curve Fitting problem MATLAB

Question

I have a system of differential equations with some unknown parameters and I need to find the optimal parameter set that fits best with the data I have.

For each parameter set choice (using for loops), my code solves the system of DE numerically, then extracts from the time vector a time subvector that agrees with data time values. The next step requires to extract a new output subvector for each of the variables that match with the time subvector (in order to compute the error).

Here is how my code is looking:

OO=[];
y2=(y(:,2))';
t1=t';
FF=[];
y1=(y(:,1))';

global "parameters and other global variables"

for parameter_1
    for parameter_2
        [t,y]=ode45('system',[time interval], [initial conditions]);
        time; % "time" is a code that extracts the time subvector

        for i=1:length(time)
            for j=1:length(t1)
                if time(i)==t1(j);
                    FF(end+1) = y1(j);
                end
            end
        end

        for k=1:length(time)
            for m=1:length(t1)
                if time(k)==t1(m);
                    OO(end+1) = y2(m);
                end
            end
        end

    end
end

At this point the code is not over, but I would be expecting it to return the OO vector and the FF vector for the last parameter set choice of the loop. However when the code runs (without any errors) it produces vectors OO and FF that are longer than vector time: length(time) < {length(OO),length(FF)}.

Answer 1

This is a more general nonlinear optimization problem than would usually be described as simply "curve fitting". I suggest you approach it as such. You will need to define a few things:

1. The parameters you are optimizing.
2. The quantity to be minimized.
3. Helper functions to simplify things.

Matlab has tools that can help and it will help you to break the problem into a couple parts as follows:

Main script

clear

% load the curve and time values you want the ODE to fit
[ tReference, yReference ] = loadKnownValues();

% load the initial conditions for the ODE:
initialConditions = loadInitialConditions();

% initial guess for parameter values (a column vector)
parametersGuess = [ 1; 6.83; ... ];

% define an objective function for the optimization
% I will assume you want to minimize the sum of the errors squared
% lsqnonlin is in the Optimization toolbox
%
% You could also use fminsearch. In that case you would want to use the
% objective function F = @(param) sum( objectiveFunction(param).^2 );
objectiveFunction = @(parameters) yReference ...
    - solveODE(tReference,initialConditions,parameters);

optimizedParameters = lsqnonlin(objectiveFunction,parametersGuess);

This main script relies on a few other functions that do specific things. Obviously you will need to load in the curve you're trying to fit and the initial conditions for the ODE (loadKnownValues() and loadInitialConditions()). You will also need a function that solves the ODE given a set of initial conditions and parameters and returns the values at the desired times. Those first two should be easy, the last one can be done as follows:

solveODE.m

function yVals = solveODE(tReference,initialConditions,parameters)
% yVals = solveODE(tReference,initialConditions,parameters)
% Solves the ODE for your problem and with the given initialConditions
% and parameters and returns the values at the specified times in
% tReference. The first time in tReference must correspond to the
% initial time.

% function that returns the right hand side of the ODE,
% use the parameters in the input vector "parameters"
rightHandSide = @(t,y) ... ;

[~, yVals] = ode45(rightHandSide, tReference, initialConditions);

Obviously I can't comment on the convergence or performance of the optimization without knowing a lot more about the system, but this should get you started.