MATLAB: How to estimate coefficients for a system of ODE’s

coefficientsestimateestimatesodesystem 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?

Best Answer

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.

Related Solutions

MATLAB: Solving nonlinear equations using fmincon

The error is straightforward:

Error using fmincon (line 635)
Supplied objective function must return a scalar value.

This is the reason fsolve works, since it is a root-finding algorithm.

In order for this to work with fmincon (and certain other optimisation functions), you would need to do something like this:

[x,fval,exitflag,output] = fmincon(@(x)norm(Obj_Fun(x)),x0,A,b,Aeq,beq,lb,ub)

That will llikely work with fmincon. Whether it produces the desired result is less certain.

MATLAB: Finding parameters of ODEs using lsqcurvefit

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.

Best Answer

Related Solutions

MATLAB: Solving nonlinear equations using fmincon

MATLAB: Finding parameters of ODEs using lsqcurvefit

Related Question