MATLAB Answers

Parameter Estimation for a System of Differential Equations

87 views (last 30 days)
Daphne Deidre Wong
Daphne Deidre Wong on 22 Sep 2020 at 4:06
Commented: Star Strider on 6 Oct 2020 at 11:06
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 Comments

Show 4 older comments
Alan Stevens
Alan Stevens on 24 Sep 2020 at 7:00
According to the model the sum of the concentrations is 1 for all time, not just at t = 0. To see this, add the four rates. They sum to zero, which means the sum of the concentrations is constant. (In fact it's not difficult to solve the equatons analytically). So, yes, I agree, we're not seeing the whole system.
Daphne Deidre Wong
Daphne Deidre Wong on 24 Sep 2020 at 15:27
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.

Sign in to comment.

Accepted Answer

Star Strider
Star Strider on 22 Sep 2020 at 14:23
Edited: Star Strider on 22 Sep 2020 at 14:37
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
.

  22 Comments

Star Strider
Star Strider on 5 Oct 2020 at 11:00
If you want to post your code, I will of course look at it and do my best to answer whatever questions you have about it.
I do not understand about ‘time constraint’ however. Please explain that in a bit more detail.
Daphne Deidre Wong
Daphne Deidre Wong on 6 Oct 2020 at 4:17
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
Star Strider
Star Strider on 6 Oct 2020 at 11:06
I fixed the datetime problem 5 days ago in my previous Comment!
Change ‘DT_GA_Time’ to this:
DT_GA_Time = datetime([0 0 0 0 0 GA_Time], 'Format','HH:mm:ss.SSS');
and the problem ino longer exists.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by