Parameter Estimation for a System of Differential Equations
Show older comments
The following set of ODE's have been previously used in a research paper to produce best fit curves for concentration vs time data for 3 reactants [3], [4] and [5]:
I am trying to reproduce these best fit curves. Usually I would try to solve these ODE's and then take the sum of the square error's of each with the actual y data and then minimze the output in order to solve for k1, k'-1 and k3. However, solving these ODE's has proven to be extremely chaotic (as demonstrated previously by Walter Roberson) and I'm not sure if my whole approach is best given the complexities. Is there another way to go about doing this?
This is the concentration and time data:
%Consider [3] = a, [4] = b and [5] = c
%x values for a & b
t = [0 2.52 4.42 11.99 22.08 32.18 53.00 73.82 99.05 124.92 155.21 191.17 227.13 278.23 380.44];
T = t';
%y values for a
a_ydata = [0.26 0.27 0.27 0.25 0.21 0.17 0.16 0.13 0.10 0.08 0.05 0.04 0.03 0.02 0.01];
A_Ydata = a_ydata';
%y values for b
b_ydata = [0 0 0.01 0.03 0.04 0.05 0.08 0.11 0.14 0.16 0.15 0.17 0.17 0.18 0.18];
B_Ydata = b_ydata';
%x and y values for c
x_c = [0 99.05 124.92 155.21 191.17 227.13 278.23 380.44];
T_C = x_c';
y_c = [0 0.01 0.015 0.018 0.023 0.024 0.028 0.029];
C_Ydata = y_c';
If anyone is able to help me, I would greatly appreciate it!
Accepted Answer
Star Strider
on 9 Nov 2020
0 votes
See the identically-titled Parameter Estimation for a System of Differential Equations for a working solution. You will need to adapt it to your data.
7 Comments
BobbyJoe
on 9 Nov 2020
Star Strider
on 9 Nov 2020
You may not be doing anything wrong. Your system of differential equations is likely encountering a singularity (-Inf or +Inf) at that point. It will likely be necessary to check the system of differential equations to be certain they are correctly written. It also could be as simple as choosing different initial parameter estimates for ‘theta0’.
I would choose:
theta0 = rand(3,1);
instead, to see if that works better, since I suspect the actual kinetic constants are likely less than 1, if my experience is relevant.
If you are still having problems after further experimentation, I will use the ga (genetic algorithm) function with your system to see if it can estimate appropriate parameters for it.
BobbyJoe
on 9 Nov 2020
Star Strider
on 9 Nov 2020
Edited: Star Strider
on 10 Nov 2020
I just now had time to use ga on your system.
I got these parameter estimates:
Rate Constants:
Theta(1) = 0.00922689
Theta(2) = 0.00038813
Theta(3) = 0.01192939
and this plot:
That appears to be a reasonably decent fit! Use these parameter estimates as ‘theta0’ in your code with lsqcurvefit possibly to get even more accurate and precise estimates.
My ga code for your problem:
function parameterestimation_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(theta,t)
c0=[0.26;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
dcdt=zeros(3,1);
dcdt(1)=-theta(1).*c(1)+theta(2).*c(2)-theta(3).*c(1).*c(2);
dcdt(2)=theta(1).*c(1)-theta(2).*c(2)-theta(3).*c(1).*c(2);
dcdt(3)=theta(3).*c(1).*c(2);
dC=dcdt;
end
C=Cv;
end
T = [0 2.52 4.42 11.99 22.08 32.18 53.00 73.82 99.05 124.92 155.21 191.17 227.13 278.23 380.44];
t = T';
%y values for a
a_ydata = [0.26 0.27 0.27 0.25 0.21 0.17 0.16 0.13 0.10 0.08 0.05 0.04 0.03 0.02 0.01];
A_Ydata = a_ydata';
%y values for b
b_ydata = [0 0 0.01 0.03 0.04 0.05 0.08 0.11 0.14 0.16 0.15 0.17 0.17 0.18 0.18];
B_Ydata = b_ydata';
c_ydata = [0 0 0.001 0.001 0.0015 0.002 0.0034 0.0069 0.0089 0.011 0.0168 0.0188 0.0232 0.0242 0.029];
C_Ydata = c_ydata';
c = [A_Ydata B_Ydata C_Ydata];
ftns = @(theta) norm(c-kinetics(theta,t));
PopSz = 500;
Parms = 3;
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([0 0 0 0 0 GA_Time], '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) = %11.8f\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')
lgdc = compose('C_%d(t)', 1:size(c,2));
legend(hlp, lgdc, 'Location','N')
format short eng
end
.
BobbyJoe
on 9 Nov 2020
Star Strider
on 9 Nov 2020
My pleasure!
Unfortunately, no.
All the vectors must have the same number of elements in order to form the matrix.
The correct way to do that is to interpolate it:
c_ydata = interp1(x_c, y_c, t)
This has the disadvantage of ‘creating’ data where none previously existed, however I doubt if that would have any significant effect on the estimated parameters.
I commented-out the original ‘c_ydata’ you provided and added these lines to the function I posted:
x_c = [0 99.05 124.92 155.21 191.17 227.13 278.23 380.44];
y_c = [0 0.01 0.015 0.018 0.023 0.024 0.028 0.029];
c_ydata = interp1(x_c, y_c, t);
C_Ydata = c_ydata(:);
and got these parameters:
Rate Constants:
Theta(1) = 0.00912635
Theta(2) = 0.00086229
Theta(3) = 0.01418786
with a slightly better fit. Use any interpolation method you want, I used the default 'linear' method here. The 'pchip' method may be better. I will let you explore that, since it did not appear to make much difference in the tests I ran with it.
If my Answer helped you solve your problem, please Accept it!
.
Star Strider
on 10 Nov 2020
As always, my pleasure!
To change it in your computer, just copy it from my previous Comment and then paste it manually to a new tab in your MATLAB Editor. You can then name the function anything you want, and it will automatically be saved with the function name you provide.
If you want me to change it in the Comment I posted, I can easily do that. I changed it to: parameterestimation_GA. You should see that change when you next acess this thread. I will change it in my system as well, in the event that I want to post it as an example in the future.
More Answers (0)
Categories
Find more on Nonlinear Grey-Box Models in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
