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

Question

Danny Helwegen 2019-1-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/441363-how-to-estimate-coefficients-for-a-system-of-ode-s

评论： Nisrina Pargustan 2020-11-21

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?

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Nisrina Pargustan 2020-11-21

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 2019-1-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/441363-how-to-estimate-coefficients-for-a-system-of-ode-s#answer_357920

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.

9 个评论
显示 7更早的评论隐藏 7更早的评论

Danny Helwegen 2019-1-25

在 MATLAB Online 中打开

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);

Danny Helwegen 2019-1-25

Than the following plot is obtained:

Star Strider 2019-1-25

在 MATLAB Online 中打开

Igor_Moura_GA.m

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.

请先登录，再进行评论。

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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

9 个评论
显示 7更早的评论隐藏 7更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

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

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

9 个评论 显示 7更早的评论隐藏 7更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

9 个评论
显示 7更早的评论隐藏 7更早的评论