Finding parameters of ODEs using lsqcurvefit

I have written the following function, and then trying to call this function using lsqcurvefit to fit the experimental data to my ordinary differential equations. it seems that I have a problem in calling lsqcurvefit. Please take a look at my codes and tell me what my problem is. Thank you.

here is my function:

function S = Kinetics(B, t)
% KINETICS codes the system of differential equations describing 
%  COD and FCL behavior in the washing system:
%       dO/dt = k0;
%       dC/dt = k1*O*C;
%   with:
%       Variables:  x(1) = O,  x(2) = C
%       Parameters: k0 = B(1),  k1 = B(2)
x0 = [305,25];
[T,Sv] = ode45(@DifEq, t,x0);
      function dS = DifEq(t, x)
          xdot = zeros(2,1);
          xdot(1) = B(1);
          xdot(2) = -B(2) .* x(2) .* x(1);
          dS = xdot;
      end
S = Sv;
end

     And this is how I have used "lsqcurvefit ":

Time = [0 2 4 6 8 10 12];
COD=[307.18 394.39 441.93 516.62 565.13 636.74 653.68];
Fcl = [21.06666667 15.4 9.633333333 4.666666667 0.753333333 0.403333333 0.206666667];
B0 = [COD(1),Fcl(1)];

[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@Kinetics, B0, Time(:), COD(:), Fcl(:));
parameters = B

[Jan]

 it seems that I have a problem

Then please explain the problem you observe.

[Ryan Abnavi]

well, my problem is that the following errors comes up: Error using lsqcurvefit (line 251) Function value and YDATA sizes are not equal. and I don't know how to fix it

[Alex Sha]

Please refer to the results below, stable and unique：

Root of Mean Square Error (RMSE): 11.0726433185545
Sum of Squared Residual: 1716.44802083901
Correlation Coef. (R): 0.992799844135045
R-Square: 0.98565153051457

Parameter          	Best Estimate
--------------------	-------------
k0	28.482650276005
k1	-0.000561833672907692

Answer 1

[Star Strider]

You are actually fitting two vectors, ‘COD’ and ‘Fcl’. (That was not immediately obvious.) They have to be in one matrix, that you then pass to lsqcurvefit as the dependent variables you want to fit. You also have to return both columns of ‘Sv’ to fit them.

This works, however you may need to experiment with different initial values to get parameters that fit well (or just run it a few times with the same initial parameters to see if different runs produce different results, as they often will):

function S = Kinetics(B, t)
% KINETICS codes the system of differential equations describing 
%  COD and FCL behavior in the washing system:
%       dO/dt = k0;
%       dC/dt = k1*O*C;
%   with:
%       Variables:  x(1) = O,  x(2) = C
%       Parameters: k0 = B(1),  k1 = B(2)
x0 = [305,25];
[T,Sv] = ode45(@DifEq, t,x0);
      function dS = DifEq(t, x)
          xdot = zeros(2,1);
          xdot(1) = B(1);
          xdot(2) = -B(2) .* x(2) .* x(1);
          dS = xdot;
      end
S = Sv;
end

%      And this is how I have used "lsqcurvefit ":

Time = [0 2 4 6 8 10 12];
COD=[307.18 394.39 441.93 516.62 565.13 636.74 653.68];
Fcl = [21.06666667 15.4 9.633333333 4.666666667 0.753333333 0.403333333 0.206666667];
COD_Fcl = [COD(:) Fcl(:)];
B0 = [COD(1),Fcl(1)];

[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@Kinetics, B0, Time(:), COD_Fcl);
parameters = B

Tv = linspace(min(Time), max(Time), 25);
Sv = Kinetics(B,Tv);

figure
plot3(Time, COD, Fcl, 'p')
hold on
plot3(Tv, Sv(:,1), Sv(:,2), '-r')
hold off
grid

This takes a very long time to run. You might get better results using your own fitness function based on your ODE system, and using the Global Optimization Toolbox genetic algorithm ga function to estimate the parameters.

If you want to let the optimization algorithm determine the initial conditions for the ODEs as well as the parameters, change ‘Kinetics’ to:

function S = Kinetics(Bv, t)
% KINETICS codes the system of differential equations describing 
%  COD and FCL behavior in the washing system:
%       dO/dt = k0;
%       dC/dt = k1*O*C;
%   with:
%       Variables:  x(1) = O,  x(2) = C
%       Parameters: k0 = B(1),  k1 = B(2)
x0 = Bv(3:4);
B = Bv(1:2);
[T,Sv] = ode45(@DifEq, t,x0);
      function dS = DifEq(t, x)
          xdot = zeros(2,1);
          xdot(1) = B(1);
          xdot(2) = -B(2) .* x(2) .* x(1);
          dS = xdot;
      end
S = Sv;
end

and ‘B0’ to:

B0 = rand(4,1);

or some vector with 4 elements.

[Ryan Abnavi]

Thank you so much, as you said it takes a very long time to run, but I do not know how to fix it, I mean I didn't exactly understand what you said, can you please explain it more?

[Star Strider]

As always, my pleasure.

I am not certain what you want me to explain. Your ODE system produces a 2-column matrix, and you are fitting 2 vectors to it. You need to fit the 2 columns to your data, so your data (your dependent variable), needs to be in a matching 2-column matrix. I created ‘COD_Fcl’ to do that.

I had problems fitting your model to your data. The problem may be that your model does not explain your data well, or that there was a problem coding the model. I am not familiar with your experiment, or what you are modeling, so I cannot be certain.

A gradient-descent approach to estimating the parameters, such as used in lsqcurvefit, may not be the best approach here. (Gradient-descent algorithms can be trapped in a local minimum, and not discover the global minimum that provides the best parameter estimates.) Optimization procedures that do not use a gradient-descent approach, such as the ga (genetic algorithm) function, may provide better parameter estimates. It searches the entire parameter space (if you design your initial population correctly), and may — probably over a day or more with your model — eventually discover a much better set of parameters. (No promises.)

Also, I always let the optimization algorithm determine the initial conditions for the ODE system, so I added the option to your parameter vector. The first 2 elements are your ‘B’ parameters, and the last 2 are the initial conditions. I usually get a much better fit this way. I included that in the later version of ‘Kinetics’.

I will of course provide more details if I still have not adequately explained my code and my approach.