Parameter Estimation for a System of Differential Equations

Question

Daphne Deidre Wong 2020-9-22

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/597652-parameter-estimation-for-a-system-of-differential-equations

回答： Jeremy Huard 2024-11-19，12:24

Hello. I am new to matlab and I have a set of kinetic equations that need to be solve. The reaction has 1 reactant in which it will decompose into 3 other product over time. I have the experimental results that shows the degradation of reactant and the concentration of products over a period of time, as below. The differential equations I need to solve are also as below.

Reaction Kinetics:

d[A]/dt = -k1[A]-k2[A]

d[B]/dt = k1[A]-k3[B]-k4[B]

d[C]/dt = k2[A]+k3[B]

d[D]/dt = k4[B]

Experimental result

Time (min) / Conc (g/L) [A] [B] [C] [D]

0 1.000 0 0 0

20 0.8998 0.0539 0.0039 0.0338

30 0.8566 0.0563 0.00499 0.0496

60 0.7797 0.0812 0.00715 0.0968

90 0.7068 0.0918 0.00937 0.1412

120 0.6397 0.0989 0.01028 0.1867

I have browse through various code and the one solved StarStrider for Igor Moura shows the result that I wish to achieve where plotted is the graph of the concentrations of reactant and products over time as well as solving the reaction rate constants "ki". However I have modified the code according to my equations and the graph came out weird and the results are not fitted well into with the experimental results. Can someone help me with this? Attach below is the code I have tried out.

function Igor_Moura
% 2016 12 03
% NOTES: 
% 
%   1. The ‘theta’ (parameter) argument has to be first in your 
%       ‘kinetics’ funciton,
%   2. You need to return ALL the values from ‘DifEq’ since you are fitting
%       all the values
function C=kinetics(k,t)
c0=[1;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(4,1);
    dcdt(1)=-k(1).*c(1)-k(2).*c(1);
    dcdt(2)= k(1).*c(1)-k(3).*c(2)-k(4).*c(2);
    dcdt(3)= k(2).*c(1)+k(3).*c(2);
    dcdt(4)= k(4).*c(2);
    dC=dcdt;
    end
C=Cv;
end
t=[20
    30
    60
    90
    120];
c=[0.91257 0.02086 0.00939 0.00725
    0.88232 0.02531 0.01664  0.00897
   0.83324 0.02882 0.03927 0.01195
  0.76289 0.03137 0.06834 0.01514  
  0.70234 0.03380 0.10542 0.01679 ]; 
k0=[1;1;1;1];
[k,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,k0,t,c);
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(k)
    fprintf(1, '\t\k(%d) = %8.5f\n', k1, k(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics(k, tv);
figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit);
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'Location','N')
end

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

Daphne Deidre Wong 2020-9-24

Sorry for the late reply. Indeed the experimental result do not show the whole system as the experiment was of a decomposition reaction that forms multiple product and side product. The products were analyze using mass spectrometry and only significant products concentration was taken into consideration.

Star Strider 2020-9-24

Daphne Deidre Wong — Thank you!

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2020-9-22

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/597652-parameter-estimation-for-a-system-of-differential-equations#answer_498550

编辑：Star Strider 2020-9-22

在 MATLAB Online 中打开

Using this version of ‘kinetics’ (that also estimates the initial conditions for ‘DifEq’),

function C=kinetics(k,t)
c0=k(5:8);
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(4,1);
    dcdt(1)=-k(1).*c(1)-k(2).*c(1);
    dcdt(2)= k(1).*c(1)-k(3).*c(2)-k(4).*c(2);
    dcdt(3)= k(2).*c(1)+k(3).*c(2);
    dcdt(4)= k(4).*c(2);
    dC=dcdt;
    end
C=Cv;
end

I got a good fit with these parameters (note that the first 4 are ‘k(1)’ through ‘k4’ and the last 4 are the initial conditions):

	Rate Constants:
		Theta(1) =  0.00180
		Theta(2) =  0.00049
		Theta(3) =  0.02622
		Theta(4) =  0.01094
		Theta(5) =  0.89995
		Theta(6) =  0.00011
		Theta(7) =  0.00213
		Theta(8) =  0.00436

producing this plot:

I used the ga (genetic algorithm) function to find them. I am cleaning up the code, and will post the entire code file in a few minutes.

EDIT — (22 Sep 2020 at 14:36)

The complete (cleaned) code:

function Daphne_Deidre_Wong_GA
% 2020 09 22
% NOTES: 
% 
%   1. The ‘k’ (parameter) argument has to be first in your 
%       ‘kinetics’ function,
%   2. You need to return ALL the values from ‘DifEq’ since you are fitting
%       all the values
function C=kinetics(k,t)
c0=k(5:8);
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(4,1);
    dcdt(1)=-k(1).*c(1)-k(2).*c(1);
    dcdt(2)= k(1).*c(1)-k(3).*c(2)-k(4).*c(2);
    dcdt(3)= k(2).*c(1)+k(3).*c(2);
    dcdt(4)= k(4).*c(2);
    dC=dcdt;
    end
C=Cv;
end
format long
t=[ 20
    30
    60
    90
   120];
c=[0.91257 0.02086 0.00939 0.00725
   0.88232 0.02531 0.01664 0.00897
   0.83324 0.02882 0.03927 0.01195
   0.76289 0.03137 0.06834 0.01514  
   0.70234 0.03380 0.10542 0.01679]; 
ftns = @(theta) norm(c-kinetics(theta,t));
PopSz = 500;
Parms = 8;
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-5, 'MaxGenerations',2E3, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
[theta,fval,exitflag,output] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],opts)
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
GA_Time = etime(t1,t0)
DT_GA_Time = datetime(num2str(GA_Time,'%.6f'),'InputFormat','ss.SSSSSS', 'Format','HH:mm:ss.SSS')
fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);
figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
    CV(k1,:) = hd(k1).Color;
    hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
    hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'Location','E')
format short eng
end

.

30 个评论
显示 28更早的评论隐藏 28更早的评论

Star Strider 2020-9-25

Hi, Daphne!

The only suggestion I have with respect to copyright is parts c. and d. from the Content (section 2) of MATLAB Central Terms of Use.

If you have questions beyond that, use the Contact Us telephone handset icon in the upper right corner of this page. I suppose this would be a Technical Support problem. Another option is to get legal assistance from your university. Copyright law (and most other aspects of the law) are very far from my areas of expertise.

With respect to adding another reaction, since I have no idea what the system is or what the parameters or compartments are, I cannot determine that. My modeling of such systems is mostly with respect to endocrine and metabilic systems, where the pathways are known and the only unknowns are some of the concentrations (in a partially ‘uncontrollable’ and ‘unobservable’ system), and most of the parameters. I may be able to help if I understand more about the system you are modeling. This assumes that I have sufficient knowledge of it to understand how to include the additional reaction.

Daphne Deidre Wong 2020-10-6

在 MATLAB Online 中打开

Currently the code that i have modify is shown below. However, the problem I have with it is the datetime error.

function Daphne_Deidre_Wong_GA_V1
% 2020 09 22
% NOTES: 
% 
%   1. The ‘k’ (parameter) argument has to be first in your 
%       ‘kinetics’ function,
%   2. You need to return ALL the values from ‘DifEq’ since you are fitting
%       all the values
function C=kinetics(k,t)
c0=k(5:9);
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(5,1);
    dcdt(1)=k(1).*c(2)-k(2).*c(1)-k(3).*c(1)-k(4).*c(1);
    dcdt(2)= k(2).*c(1)-k(1).*c(2)-k(5).*c(2)-k(6).*c(2);
    dcdt(3)= k(3).*c(1)+k(5).*c(2);
    dcdt(4)= k(6).*c(2);
    dcdt(5)= k(4).*c(1);
 
    
    dC=dcdt;
    end
C=Cv;
end
format long
t=[ 20
    30
    60
    90
   120];
c=[0.89978 0.05385 0.00387 0.03384 0.00499
    0.85655 0.05634 0.00499 0.0496 0.00665
    0.77968 0.08119 0.00715 0.09682 0.00929
    0.70684 0.09177 0.00937 0.14118 0.0114
    0.63967 0.09893 0.01028 0.18666 0.01277]; 
ftns = @(theta) norm(c-kinetics(theta,t));
PopSz = 500;
Parms = 9;
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-5, 'MaxGenerations',2E3, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
[theta,fval,exitflag,output] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],opts)
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
GA_Time = etime(t1,t0)
DT_GA_Time = datetime(num2str(GA_Time,'%.6f'),'InputFormat','ss.SSSSSS', 'Format','HH:mm:ss.SSS');
fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);
figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
    CV(k1,:) = hd(k1).Color;
    hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
    hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'C_5(t)', 'Location','E')
format short eng
end

Daphne Deidre Wong 2020-10-27

在 MATLAB Online 中打开

Hello Star Strider,

I am sorry to bother you again. I was thinking if you can help me take a look at the code below as the curve have strayed quite a lot from the data and I do not know how to fix it. I am so sorry for the inconveniences, but I really need this help.

function Gluc_200
% 2020 09 22
% NOTES: 
% 
%   1. The ‘k’ (parameter) argument has to be first in your 
%       ‘kinetics’ function,
%   2. You need to return ALL the values from ‘DifEq’ since you are fitting
%       all the values
function C=kinetics(k,t)
c0=k(5:8);
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
   dcdt=zeros(4,1);
    dcdt(1)=-k(1).*c(1)-k(2).*c(1)-k(3).*c(1);
    dcdt(2)= k(1).*c(1)-k(4).*c(2);
    dcdt(3)= k(2).*c(1)+k(4).*c(2);
    dcdt(4)= k(3).*c(1);
    dC=dcdt;
    end
C=Cv;
end
format long
t=[ 20
    40
    60
    90
   100];
c=[0.59511 0.12077 0.01160 0.01548
    0.39150 0.08768 0.01284 0.01784
    0.28869 0.05067 0.01251 0.01633
    0.19520 0.03158 0.01158 0.01317
    0.15019 0.02249 0.00883 0.01141]; 
ftns = @(theta) norm(c-kinetics(theta,t));
PopSz = 500;
Parms = 8;
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-5, 'MaxGenerations',2E3, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
[theta,fval,exitflag,output] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],opts)
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
GA_Time = etime(t1,t0)
DT_GA_Time = datetime(num2str(GA_Time,'%.6f'),'InputFormat','ss.SSSSSS', 'Format','HH:mm:ss.SSS')
fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);
figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
    CV(k1,:) = hd(k1).Color;
    hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
    hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'Location','E')
format short eng
end

Kaihua Liu 2023-6-23

sorry to bother you again, I successfully made graphs using your program, but how should I separate the graphs (Make a picture for each substance)?

Star Strider 2023-6-23

Use subplot, tiledlayout, or stackedplot. (I believe that exhausts the options.)

A Vote would be appreciated!

请先登录，再进行评论。

Answer 2

Jeremy Huard 2024-11-19，12:24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/597652-parameter-estimation-for-a-system-of-differential-equations#answer_1546858

You could also consider using SimBiology which provides parameter estimation capabilities for ODE systems.

Here is a short video tutorial.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Parameter Estimation for a System of Differential Equations

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

采纳的回答

30 个评论
显示 28更早的评论隐藏 28更早的评论

更多回答（1 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

Parameter Estimation for a System of Differential Equations

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

采纳的回答

30 个评论 显示 28更早的评论隐藏 28更早的评论

更多回答（1 个）

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

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

30 个评论
显示 28更早的评论隐藏 28更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论