Parameter Estimation for a System of Differential Equations

125 views (last 30 days)
Igor Moura
Igor Moura on 1 Dec 2016
Commented: Kishen Mahadevan on 21 Nov 2024
Hello. Ok, so I'm new to matlab and I've got a question regarding parameter estimation for a kinetic model. I have 4 different reactants and their concentrations are c1, c2, c3 and c4. I also have 4 differential equations, each one related to a concentration (c1, c2, c3 and c4, respectively -see below-) and experimental data for all these concentrations on 12 different times plus the initial condition. The k's are the rate coefficients. I want to solve this system of ODE's using ode45 and then use the output to compute the experimental data minus the observed data and use these results to estimate the values of k's using lsqnonlin, but apparently I can't solve these ODE's without numerical values for k -which is what I want to know-. Any help on how to set up the command to solve this?
function dcdt=batch(t,c,k)
dcdt=zeros(4,1);
dcdt(1)=-k(1)*c(1)-k(2)*c(1);
dcdt(2)= k(1)*c(1)+k(4)*c(3)-k(3)*c(2)-k(5)*c(2);
dcdt(3)= k(2)*c(1)+k(3)*c(2)-k(4)*c(3)+k(6)*c(4);
dcdt(4)= k(5)*c(2)-k(6)*c(4);
end
Data:
t c1 c2 c3 c4
0 1 0 0 0
0.1 0.902 0.06997 0.02463 0.00218
0.2 0.8072 0.1353 0.0482 0.008192
0.4 0.6757 0.2123 0.0864 0.0289
0.6 0.5569 0.2789 0.1063 0.06233
0.8 0.4297 0.3292 0.1476 0.09756
1 0.3774 0.3457 0.1485 0.1255
1.5 0.2149 0.3486 0.1821 0.2526
2 0.141 0.3254 0.194 0.3401
3 0.04921 0.2445 0.1742 0.5277
4 0.0178 0.1728 0.1732 0.6323
5 0.006431 0.1091 0.1137 0.7702
6 0.002595 0.08301 0.08224 0.835
Thanks in advance!
  7 Comments
Show 5 older comments
Budda
Budda on 2 Jan 2021
Hi Star Srider, How would I change the code if only say the data for first variable ( equation) available in the Igor_Moura.m?
Thank you in advance
Star Strider
Star Strider on 2 Jan 2021
Budda —
Make appropriate changes to the ‘kinetics’ function (specifically to the ‘C’ variable) to output only the variable you need. In that code, there are 4 columns, so address only the column you want returned in ‘C’.

Sign in to comment.

Sign in to answer this question.

Accepted Answer

Star Strider
Star Strider on 1 Dec 2016
See if the techniques in Monod kinetics and curve fitting will do what you want.
  52 Comments
Show 50 older comments
Star Strider
Star Strider on 3 Aug 2024
Ran in:
This runs, though it tales too long to run here.
I created:
SetCell{1} = {t1, c1, temp1, c0_1};
SetCell{2} = {t2, c2, temp2, c0_2};
SetCell{3} = {t3, c3, temp3, c0_3};
and then a loop to loop through them. Note that I included the various ‘temp’ variables in my cell array.
I also corrected several errors, all of which were introduced here and are not part of my original code.
If the first seet of parameters produces a decent fit to the datta, preserve them and use them as the initial estimates for the subsequent sets. It may make this more efficient.
A Vote for my Answer here would be appreciated!
PE_LScurvefit_Multiple_data_sets_non_isothermal
Warning: Length of lower bounds is < length(x); filling in missing lower bounds with -Inf.
Warning: Length of upper bounds is < length(x); filling in missing upper bounds with +Inf.
Local minimum possible. lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance. Rate Constants: Theta( 1) = 0.05025 Theta( 2) = 0.08310 Theta( 3) = 0.66899 Theta( 4) = 0.66844 Theta( 5) = 0.00006 Theta( 6) = 0.19195 Theta( 7) = 1367.83624 Theta( 8) = 2619.94112 Theta( 9) = 2964.46774 Theta(10) = 13236.98878 Theta(11) = -3738.19745 Theta(12) = 529.43631
Warning: Length of lower bounds is < length(x); filling in missing lower bounds with -Inf.
Warning: Length of upper bounds is < length(x); filling in missing upper bounds with +Inf.
Local minimum possible. lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance. Rate Constants: Theta( 1) = 0.00908 Theta( 2) = 0.71083 Theta( 3) = 0.43675 Theta( 4) = 0.75300 Theta( 5) = 0.00062 Theta( 6) = 0.61537 Theta( 7) = -58.95522 Theta( 8) = 3751.88913 Theta( 9) = 2646.80729 Theta(10) = 17763.40859 Theta(11) = -2821.47869 Theta(12) = 462.55345
Warning: Length of lower bounds is < length(x); filling in missing lower bounds with -Inf.
Warning: Length of upper bounds is < length(x); filling in missing upper bounds with +Inf.
Local minimum possible. lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance. Rate Constants: Theta( 1) = 0.00000 Theta( 2) = 0.00232 Theta( 3) = 0.11307 Theta( 4) = 0.78547 Theta( 5) = 0.86578 Theta( 6) = 0.90149 Theta( 7) = 8.45868 Theta( 8) = 9.95011 Theta( 9) = -440.41576 Theta(10) = 103.26330 Theta(11) = 59.29220 Theta(12) = 0.24303
Elapsed time is 3.514244 seconds.
function PE_LScurvefit_Multiple_data_sets_non_isothermal
tic
%Attempt to update from 1 dataset (t,c) and isothermal operation to
% three data sets (t1,c1), (t2,c2), (t3,c3), each set at different
% temperatures, meaning now theta is theta_temp, so theta_temp(1),
% theta_temp(2) and theta_temp(3) have different values.
%So, the parameter to be estimated are k(1 to 6) and E(1 to 6).
%The plotting part is not updated because the code is not running
%T14 oxidation 80ºC 6bar
t1=[0
16
30
45
61
96
121
180
240
362
492
640];
%T14 oxidation 80ºC 6bar
c1=[105.84 0.00 0.00 0.00
83.03 17.36 2.29 2.47
74.82 19.85 1.63 3.61
69.28 21.98 1.33 4.81
62.72 25.19 1.26 6.18
51.98 28.93 0.82 9.16
43.11 29.59 0.60 9.97
28.58 28.07 0.51 12.28
17.87 24.32 0.53 12.98
7.84 16.22 0.53 11.40
2.68 8.22 0.83 6.92
0.66 2.77 0.28 0.85];
%T24 oxidation 90ºC 6bar
t2=[0
15.25
30
45
60
90
120
180
240
300
360
412
];
%T24 oxidation 90ºC 6bar
c2=[110.16 0.00 0.00 0.00
79.13 17.04 1.22 2.82
68.77 21.16 0.88 4.74
56.83 25.07 0.74 6.46
45.17 27.84 0.71 8.03
35.49 28.68 0.87 11.40
20.73 25.15 0.49 10.94
9.77 18.07 0.54 11.40
2.78 10.48 0.00 6.45
0.98 5.35 0.00 3.71
0.42 2.13 0.00 1.69
0.35 0.76 0.00 0.89
];
% T3 oxidation 70ºC 6bar DIFFERENT SIZED
t3 = [0; 50; 100; 150; 200; 250; 300; 350; 400];
c3 = [120.5 0.0 0.0 0.0
110.2 5.6 0.1 4.3
100.8 3.2 0.2 2.1
90.5 2.1 0.3 1.5
80.2 1.5 0.4 1.0
70.9 1.1 0.5 0.7
60.6 0.8 0.6 0.5
50.3 0.6 0.7 0.4
40.0 0.4 0.8 0.3];
temp1 = 80; % Temperature for data set 1
temp2 = 90; % Temperature for data set 2
temp3 = 70; % Temperature for data set 3
% Initial conditions
c0_1 = [105.84; 0; 0; 0];
c0_2 = [110.16; 0; 0; 0];
c0_3 = [120.5; 0; 0; 0];
SetCell{1} = {t1, c1, temp1, c0_1};
SetCell{2} = {t2, c2, temp2, c0_2};
SetCell{3} = {t3, c3, temp3, c0_3};
% Combine the additional data sets and initial conditions
% t = {t1, t2, t3};
% c = {c1, c2, c3};
% c0 = {c0_1, c0_2, c0_3};
% Kinetics function to include the temperature variable
function C = kinetics(theta, t, c0,temp)
[~, Cv] = ode15s(@(t,c) DifEq (t,c,temp), t, c0);
function dC = DifEq(t, c, temp)
dcdt = zeros(4, 1);
% Introduce temperature dependence on theta parameter - theta = k *
% Exp (-E/(R*temp))
k = theta(1:6);
E = theta(7:12);
R = 8.314; % Gas constant
% whos
theta_temp = k .* exp(-E ./ (R * temp));
dcdt(1) = -theta_temp(1) * c(1) - theta_temp(2) * c(1);
dcdt(2) = theta_temp(1) * c(1) - theta_temp(3) * c(2) - theta_temp(4) * c(2);
dcdt(3) = theta_temp(2) * c(1) + theta_temp(4) * c(2) - theta_temp(6) * c(3);
dcdt(4) = theta_temp(3) * c(2) - theta_temp(5) * c(4);
dC = dcdt;
end
C = Cv;
end
% Initialize theta values
theta0 = rand(12, 1);
% Boundaries for theta
lb = [0; 0; 0; 0; 0; 0];
ub = [1; 1; 1; 1; 1; 1];
% Combine the data sets
% t = {t1, t2};
% c = {c1, c2};
% c0 = {c0_1, c0_2};
% Non-linear solver
for k = 1:numel(SetCell)
t = SetCell{k}{1};
c = SetCell{k}{2};
temp = SetCell{k}{3};
c0 = SetCell{k}{4};
[theta, Rsdnrm, Rsd, ExFlg, OptmInfo, Lmda, Jmat] = lsqcurvefit(@(theta, t) kinetics(theta, t, c0, temp), theta0, t, c, lb, ub);
fprintf(1, '\tRate Constants:\n')
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit1 = kinetics(theta, tv, c0_1, temp);
% Cfit2 = kinetics(theta, tv, c0_2);
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, Cfit1);
for k1 = 1:size(c, 2)
hlp(k1).Color = CV(k1, :);
end
% hd = plot(t, c2, 'p');
% for k1 = 1:size(c2, 2)
% CV(k1, :) = hd(k1).Color;
% hd(k1).MarkerFaceColor = CV(k1, :);
% end
% hlp = plot(tv, Cfit2);
% for k1 = 1:size(c2, 2)
% hlp(k1).Color = CV(k1, :);
% end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, compose('C_%d', 1:size(c1, 2)), 'Location', 'North')
end
toc
end
The fitted values are reasonably good, and the code did not take more than a few seconds (actually about 3.514 seconds) to run.
.
Kishen Mahadevan
Kishen Mahadevan on 21 Nov 2024
I would like to add that in situations when the data is uniformly sampled, you can also use grey box estimation techniques from System Identification Toolbox to estimate these parameters.
Here is an example that highlights the steps of estimating coefficients of ODEs that describe the model dynamics to fit a given response trajectory.
Estimate Coefficients of ODEs to Fit Given Solution

Sign in to comment.

More Answers (0)

Categories

Find more on Grey-Box Model Estimation in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!