1
$\begingroup$

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)}.

$\endgroup$
  • $\begingroup$ Welcome to SciComp.SE. It is not clear what your question is. What is the problem that you have? $\endgroup$ – nicoguaro Jul 1 '15 at 21:20
  • $\begingroup$ Hi. Thanks for your reply. When I run the code. the dimension of the OO and FF vector is not the same as the time vector. My time vector has dimension 111. I need my OO and FF vector to be the same. $\endgroup$ – User1 Jul 1 '15 at 22:33
  • $\begingroup$ The length of the time subvector is not given, and the double loops you use will, for the pathological case, produce a vector of length l=length(time)*length(t1). The length of FF and OO can be of any length between 0 and my l defined before depending on the level of agreement. Maybe you could draw a picture of what you're trying to accomplish? $\endgroup$ – Bill Barth Jul 1 '15 at 23:52
  • $\begingroup$ I use ode45 to solve the system of DE over the time interval t=[0 110] However MATLAB, when computing numerically the solution, generates a discrete time interval whose elements are not equally spaced. I have a code where MATLAB picks 111 elements of the vector t (closest to integer time values). So the length of "time" is 111. When I run the above code. The length of OO and FF is 221 (they must be of same length as time) For the picture, after getting OO and FF, I want to compute the error between the simulated model and data. The lengths of the data vectors are also 111. $\endgroup$ – User1 Jul 2 '15 at 1:23
  • $\begingroup$ Finally, I get a vector whose entries are errors for each parameter set. I want to find the minimum of this vector and the only output that i want to see from the code is the choice of parameter set that minimizes the error. It will be my baseline choice for the model. $\endgroup$ – User1 Jul 2 '15 at 1:24
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Hi Doug, I decided to give your code a try. In my code, I had to "manually" zoom into the 5-dimensional grid (5 parameter space). I am also not satisfied with the results. parameters are very much correlated. I ran your code but I think I have a tolerance problem. The code stops after 1 second and says that the initial point is a local minimum. values of yReference are between 0 and 1 and so errors should be of the order 0.001 or less. I tried changing Tolerance but I couldn't manage. How should I proceed? @User1 $\endgroup$ – User1 Jul 30 '15 at 16:51
  • $\begingroup$ i used options=optimset('TolX',a,'TolFun',b) and I am letting a and b be of very small orders e.g. 10^-10 and it still takes one second. I also defined a lower and upper bound for each of the parameters to increase the grid search. Nothing is helping @Doug $\endgroup$ – User1 Jul 30 '15 at 17:27
  • $\begingroup$ If you are getting "The code stops after 1 second and says that the initial point is a local minimum" from lsqnonlin(), this is commonly due to a discontinuous objective function. A common case is where the residual function samples some data using nearest-neighbor interpolation, e.g. Residual = @(t) y - interp1(tdat,ydat,t,'nearest'); In this case, changing from 'nearest' to 'linear' will fix the problem, allowing lsqnonlin() to "see" small changes in t. $\endgroup$ – GeoMatt22 Oct 6 '15 at 1:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.