Parameter Estimation for a System of Differential Equations

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!

Answer 1

[Star Strider]

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.

[BobbyJoe]

I tried to do that with the following:

function Igor_Moura

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];
theta0=[1;1;1];

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);

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(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)', 'Location','N')

end

But then I get this error:

>> Igor_Moura
Warning: Failure at t=3.932757e+00.  Unable to meet integration tolerances without reducing the step size below the
smallest value allowed (7.105427e-15) at time t. 
> In ode45 (line 360)
  In Igor_Moura/kinetics (line 12)
  In lsqcurvefit/objective (line 278)
  In snls (line 336)
  In lsqncommon (line 164)
  In lsqcurvefit (line 271)
  In Igor_Moura (line 41) 
Matrix dimensions must agree.

Error in lsqcurvefit/objective (line 279)
         F = F - YDATA;

Error in snls (line 336)
            newfvec = feval(funfcn{3},xcurr,varargin{:});

Error in lsqncommon (line 164)
            snls(funfcn,xC,lb,ub,flags.verbosity,options,defaultopt,initVals.F,initVals.J,caller, ...

Error in lsqcurvefit (line 271)
    lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,allDefaultOpts,caller,...

Error in Igor_Moura (line 41)
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);

I'm somewhat new to MATLAB so I'm not sure what it is I'm doing wrong, any help would be greatly appreciated!

[Star Strider]

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]

I tried using

theta0 = rand(3,1);

But I still get the following errors:

>> Igor_Moura
Warning: Failure at t=4.881305e+00.  Unable to meet integration tolerances without reducing the step size below the
smallest value allowed (1.421085e-14) at time t. 
> In ode45 (line 360)
  In Igor_Moura/kinetics (line 12)
  In lsqcurvefit/objective (line 278)
  In snls (line 336)
  In lsqncommon (line 164)
  In lsqcurvefit (line 271)
  In Igor_Moura (line 41) 
Matrix dimensions must agree.

Error in lsqcurvefit/objective (line 279)
         F = F - YDATA;

Error in snls (line 336)
            newfvec = feval(funfcn{3},xcurr,varargin{:});

Error in lsqncommon (line 164)
            snls(funfcn,xC,lb,ub,flags.verbosity,options,defaultopt,initVals.F,initVals.J,caller, ...

Error in lsqcurvefit (line 271)
    lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,allDefaultOpts,caller,...

Error in Igor_Moura (line 41)
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);

Any idea what to do?

[Star Strider]

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]

WOW!!! Thank you SO much Star Strider!!! You truly deserve MVP recognition!!! I seriously appreicate all your time and effort in helping me figure this out!

I have one question to ask: In the original post when I first posted my x & y values, there were less values for 'c' compared to 'a' & 'b'. Is it possible to calculate the 'k' values even when this is the case? If so, may I ask how? Really appreciate it!

[Star Strider]

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]

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.