How to estimate coefficients for a system of ODE's?

Hi guys, I need some help with how to process the code in matlab. I have set up a system of ODE's, got the initial values and have values for A at specific times and now I need to estimate the values of k1, k2 and k3. I tried several things, but I still don't have a clue on how to process this into matlab, this is what i got:

ODE's

dAdt = -k1*A*B-k2*A*C-k3*A*D
dBdt = -k1*A*B
dCdt =  k1*A*B-k2*A*C
dDdt =  k1*A*B+k2*A*C-k3*A*D
dEdt =  k3*A*D
dFdt =  k2*A*C+k3*A*D

Inital values

A = 20.09; B = 6.96; C = 0; D = 0; E = 0; F = 0

Specific values for A at different times:

numdata = xlsread('dataset.xlsx')
t = numdata(:,1);
A = numdata(:,2);

eg. numdata =

    4.5000   15.4000
    8.6700   14.2200
   12.6700   13.3500

Is there someone that can help me?

[Nisrina Pargustan]

hai danny, i got the exact same problem as you. so i tried to run your code but i got the error message. can you please send me your matlab code for this problem?

Answer 1

[Star Strider]

If you are doing parameter estimation of a system of ordinary differential equations, these will likely provide some guidance:

Monod kinetics and curve fitting (link)

Solving Coupled Differential Equations (link)

Parameter Estimation for a System of Differential Equations (link)

This is usually a fairly straightforward problem.

[Danny Helwegen]

Hey thanks for your answer, I used the links and set up my code. Only, when i execute my code i get the following error:

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the default value of the function tolerance.

<stopping criteria details>


Optimization stopped because the relative sum of squares (r) is changing
by less than options.FunctionTolerance = 1.000000e-06.

Optimization Metric                                  Options
relative change r =   2.64e-08             FunctionTolerance =   1e-06 (default)

Any idea where the mistake can be? This is my written code:

function A = Kinetics(B, t)
% System of ODE's:
%       dAdt=-k1*A*B - k2*A*C - k3*A*D
%       dBdt=-k1*A*B
%       dCdt= k1*A*B - k2*A*C
%       dDdt= k1*A*B + k2*A*C - k3*A*D
%       dEdt= k3*A*D
%       dFdt= k2*A*C + k3*A*D
%       
%   with:
%       Variables:  x(1) = A,  x(2) = B, x(3) = C, x(4) = D, x(5) = E, 
%       x(6) = F
%       Parameters: k1 = B(1),  k2 = B(2),  k3 = B(3)
x0 = [20.09; 6.96; 0; 0; 0; 0];
[T,Av] = ode45(@DifEq, t, x0);
      function dA = DifEq(t, x)
          xdot = zeros(6,1);
          xdot(1) = (-B(1) .* x(1) .* x(2) - B(2) .* x(1) .* x(3) - B(3) .* x(1) .* x(4));
          xdot(2) = (-B(1) .* x(1) .* x(2));
          xdot(3) = (B(1) .* x(1) .* x(2) - B(2) .* x(1) .* x(3));
          xdot(4) = (B(1) .* x(1) .* x(2) + B(2) .* x(1) .* x(3) - B(3) .* x(1) .* x(4));
          xdot(5) = (B(3) .* x(1) .* x(4));
          xdot(6) = (B(2) .* x(1) .* x(3) + B(3) .* x(1) .* x(4));
          dA = xdot;
      end
  A = Av(:,1);
end

And I call it via:

B0 = rand(3,1)*100;
numdata = xlsread('ASSIGNMENT5_dataset.xlsx');
time = numdata(:,1);
A_data = numdata(:,2);
[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@Kinetics,B0,time,A_data);

[Torsten]

This is not an error message, but an information for you.

[Star Strider]

@ Danny Helwegen —

When you plot the solution (using the estimated parameters lsqcurvefit returns) and your data, how accurate is the fit?

@Torsten —

Thank you.

[Danny Helwegen]

Oh okay I get it now. I have only one question left, the function now depens on a random generated B0 that I first need to multiply with 100 too let the function work properly, is there a way to make it indipendent of this B0 because the output depends on the values of B0?

[Star Strider]

You can set ‘B0’ to be anything you want. It is best if you set it to values you believe are close to the correct parameters (at least with respect to magnitude).

You can have the optimization function (here lsqcurvefit) estimate the initial values for your system of differential equations as well. I believe I did that in the ‘Monod Kinetics’ code.

You can also constrain the parameters with lsqcurvefit, so for example if they can only be positive, set the lower bound equal to a vector of zeros.

If you are having problems finding the best parameter set for your system of differential equations, and if you have the Global Optimization Toolbox, I suggest using a genetic algorithm (the ga funciton). I have some example code I can attach if you want to do that.

[Danny Helwegen]

Okay thanks, if you would attach the example code that would be really nice. And when i try to plot the solution it gives me an error that the T is not defined, so this is the solution plot:

plot(T,Av)

But i get the error: Undefined function or variable 'T'.

[Torsten]

B0 = rand(3,1)*100;
numdata = xlsread('ASSIGNMENT5_dataset.xlsx');
time = numdata(:,1);
A_data = numdata(:,2);
[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@Kinetics,B0,time,A_data);
Av = Kinetics(B,time);
plot(time,A_data,time,Av)

[Danny Helwegen]

Than the following plot is obtained:

[Star Strider]

Here it is (a version of previous code, this time using ga), attached.

Since all the parameters are positive (by definition of the model), I constrained them to be greater than or equal to zero. Change that if your parameters are allowed to be negative.

Note that you will have to re-state your objective function in terms of the fitness function (in this code, ‘ftns’):

ftns = @(B) norm(A_data - Kinetics(B, time));

You mentioned that your parameters are on the order of 100, so change the options structure to:

opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+2,PopSz,Parms), 'MaxGenerations',2E3, 'PlotFcn','gaplotbestf');

That should make it converge faster, however it will likely still take a few minutes. Also, ‘theta’ in my code is ‘B’ in yours.