MATLAB: How to properly model a kinetic reaction system from experimental data using lsqcurvefit

kineticslsqcurvefitMATLABparameter estimation

Hello,

I am currently trying to model chemical kinetics using MATLAB.

Currently, I am having a code that works but gives me a set on non-sense results. So there must be a mistake in the method, or in my way of thinking.

I will introduce the problem, the available data and my MATLAB integration. I hope you can help me find the mistakes and obtain a fit for my experimental data.

Thanks in advance.

Problem introduction:

This system consists of four reactions. Compound c1 reacts to compound c2. Compound c2 can react further to c3, although I am currently neglecting the reaction of c2 to c3. The challenge is that both c1 and c2 can undergo polymerization reactions to p1 and p2. This of course has an influence on their concentrations and therefore on kinetics. Fitting just reaction 1 (c1 to c2) and neglecting reactions 3 and 4 leads to extreme reaction orders. Also, during experiments the mass balance can never close properly if these reactions are to be neglected.

I have shared my code at the bottom of this question so you can see how I implemented the three reactions into MATLAB.

Available data:

Over a certain time period we have the concentration profiles of compounds c1 and c2. The amount of c3 is negligable, p1 and p2 can not be analyzed.

These profiles look like this:

This data is figurated however close to the experimental data. It should be well sufficient for this modeling challenge.

MATLAB Integration:

In MATLAB, I have used the lsqcurvefit method as described here: https://nl.mathworks.com/matlabcentral/answers/43439-monod-kinetics-and-curve-fitting. Thanks Star Strider for this starting point. I can now create a code that runs. However, the results from lsqcurvefit do not match the experimental data at all. Moreover, the results are different for every run.

My goal is to plot in figure 2 the experimental data and the model fit. As you will probably see, it is a mess although I am not sure why.

Your time and effort is very much apprieciated.

My code:

clear, clc
%Reaction system
%c1 --> c2 --> c3 (c2 --> c3 for now neglected)
%c1 --> p1
%c2 --> p2
%Reaction rates are described as r = k*c^n
%As c2 --> c3 is found negligable, data is available only for c1 and c2
%The last two reaction equations describe polymerization reactions, there
%is no data available for both p1 and p2.
%Data
time = [10 20 40 60 180];
c1 = [0.508 0.442 0.385 0.354 0.258];
c2 = [0.246 0.301 0.359 0.398 0.489];
c1_0 = 0.711;
c2_0 = 0;
%Plot

time_plot = [0,time];
c1_plot = [c1_0,c1];
c2_plot = [c2_0,c2];
figure(1)
plot(time_plot,c1_plot,'bo')
grid on
hold on
plot(time_plot,c2_plot,'ro')
hold off
xlabel('Time [min]')
ylabel('Concentration [mol/L]')
legend('c1','c2')
%Data processing
x_in = time';
y_in = [c1',c2'];
%Starting point for lsqcurvefit
p0 = [0.1,1,0.1,1,0.1,1];
lb=zeros(6,1);
ub=4*ones(6,1);
%Lsqcurvefit
[p,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@diff1,p0,x_in,y_in,lb,ub);
FitData = diff1(p,time,1);
c1_out = FitData(:,1);
c2_out = FitData(:,2);
%Plot
c1_out_plot = [c1_0,c1_out'];
c2_out_plot = [c2_0,c2_out'];
figure(2)
plot(time_plot,c1_plot,'bo')
grid on
hold on
plot(time_plot,c2_plot,'ro')
plot(time_plot,c1_out_plot,'b')
plot(time_plot,c2_out_plot,'r')
hold off
xlabel('Time [min]')
ylabel('Concentration [mol/L]')
legend('c1','c2')
function C = diff1(p,time,~)
%dc1/dt = -(k1*c1^n1+k3*c1^k3)
%dc2/dt = k1*c1^n1-(k2*c2^k2+k4*c2^n4)
%dc3/dt = k2*c2^n2 = 0
%variables: c1 = x(1), c2 = x(2)
%parameters: k1 = p(1), n1 = p(2), k3 = p(3), 
%n3 = p(4), k4 = p(5), n4 = p(6)
%for now we neglect reaction r2 of c2 to c3.
x0 = rand(2,1);
[t,Cout] = ode45(@DifEq, time, x0);
      function dC = DifEq(t, x)
          xdot = zeros(2,1);
          xdot(1) = -(p(1).*x(1).^p(2)+p(3).*x(1)^p(4));
          xdot(2) = p(1).*x(1).^p(2)-p(5).*x(2).^p(6);        
          dC = xdot;
      end
C = Cout;
end

Best Answer

Thank you for referring to my ’Monod Kinetics’ Answer! I’m glad you found it useful!

There are two problems with your code.

First, ‘xdot(1)’ should be:

xdot(1) = -(p(1).*x(1).^p(2)-p(3).*x(1)^p(4));
                            ↑

with the sign of the second term being negative (corrected here).

Second, it is a good idea to always include the initial values as parameters, as a rule, the last parameters. So here:

p0 = [0.1,1,0.1,1,0.1,1,c1(1),c2(1)];

and then incorporating those changes in ‘diff1’:

function C = diff1(p,time,~)
%dc1/dt = -(k1*c1^n1+k3*c1^k3)
%dc2/dt = k1*c1^n1-(k2*c2^k2+k4*c2^n4)
%dc3/dt = k2*c2^n2 = 0
%variables: c1 = x(1), c2 = x(2)
%parameters: k1 = p(1), n1 = p(2), k3 = p(3), 
%n3 = p(4), k4 = p(5), n4 = p(6)
%for now we neglect reaction r2 of c2 to c3.
x0 = p(7:8);
[t,Cout] = ode45(@DifEq, time, x0);
      function dC = DifEq(t, x)
          xdot = zeros(2,1);
          xdot(1) = -(p(1).*x(1).^p(2)-p(3).*x(1)^p(4));
          xdot(2) = p(1).*x(1).^p(2)-p(5).*x(2).^p(6);        
          dC = xdot;
      end
C = Cout;
end

the estimated parameters are:

p =
   0.105408507696630
   3.999984376449200
   0.000078501021266
   0.007282181170239
   0.002202264694465
   3.902023967475535
   0.502316713470111
   0.244896706878439

and your model fits your data almost exactly!

Related Solutions

MATLAB: Fitting data with multiple inputs, ODE equation, and lsqnonlin

The fundamental problem is that lsqnonlin ~= lsqcurvefit! They are definitely not interchangable, and have different argument lists and calling conventions, the ones used here being for lsqcurvefit, so that change was straightforward (after checking to be sure everything matched correctly).

Beyond that, there is no ‘x(3)’ since there is no third differential equation. Since ‘x(3)’ should be temperature (that I here assume is a funciton of time), I created:

Temp = interp1(x_in(:,1), x_in(:,2),t);                                               % <— ADDED & Changed ‘x(3)’ To ‘Temp’

to interpolate temperature with time, and then made the substitutions in the differential equations. If temperature is not a function of time, it will be necessary to create a function that expresses temperature in a way that is meaningful for the differential equations.

I also set the lower bounds on the parameters to zero, since concentrations and kinetic constants, at least in this system, appear to always be greater than zero.

This version runs without errors:

%Reaction system
%c1 --> c2 
%c1 --> p1
%c2 --> p2
%Reaction rates are described as r = k*c^n
%As c2 --> c3 is found negligable, data is available only for c1 and c2
%The data for p1 and p2 is unknown.
%Data
time = [10 20 40 60 100 140 180];
c1 = [0.508 0.442 0.385 0.354 0.299 0.267 0.258];
c2 = [0.246 0.301 0.359 0.398 0.422 0.467 0.489];
T_o = [400 423 433 456 467 468 463];
c1_0 = 0.711;
c2_0 = 0;
T_o_0 = 400;
R = 8.314;
%Data processing
x_in = [time(:),T_o(:)];
y_in = [c1(:),c2(:)];
%Lsqnonlin
%[p,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqnonlin(@diff1,p0,x_in,y_in,lb,ub);
%FitData = diff1(p,time,1);
%c1_out = FitData(:,1);
%c2_out = FitData(:,2);
% Options for fitting
options = optimoptions('lsqcurvefit', ...
    'Display', 'iter', ...
    'Algorithm', 'levenberg-marquardt', ...
    'UseParallel', false, ...
    'StepTolerance', 1e-6);
%initial value for iteration
B0 = [0.1,1,0.1,1,0.1,1,0.1,1,0.1,1,0.1,1,c1(1),c2(1),T_o(1)];
% lsqnonlin
optimizedParams = lsqcurvefit(@diff1, B0, x_in, y_in, zeros(size(B0)), [], options)
function C = diff1(B,x_in)
%dc1/dt = -(k1*c1^n1+k3*c1^k3)
%dc2/dt = k1*c1^n1-(k2*c2^k2+k4*c2^n4)
%where
%k1 = A1*exp(-E1/RT_o)
%k2 = A2*exp(-E2/RT_o)
%k3 = A3*exp(-E3/RT_o)
%k4 = A4*exp(-E4/RT_o)
%
%variables: c1 = x(1), c2 = x(2), T_o = x(3)
%parameters: 
%A1 = B(1), E1 = B(2), n1 = B(3)
%A2 = B(4), E2 = B(5), n2 = B(6)
%A3 = B(7), E3 = B(8), n3 = B(9)
%A4 = B(10), E4 = B(11), n4 = B(12)
time = x_in(:,1);
x0 = B(11:12);
[t,Cout] = ode45(@DifEq, time, x0);
      function dC = DifEq(t, x)
          xdot = zeros(2,1);
          Temp = interp1(x_in(:,1), x_in(:,2),t);                                               % <— ADDED & Changed 'x(3)' To 'Temp'
          xdot(1) = -(B(1)*exp(-B(2)/R/Temp).*x(1).^B(3)-B(7)*exp(-B(8)/R/Temp).*x(1)^B(9));
          xdot(2) = B(4)*exp(-B(5)/R/Temp).*x(1).^B(6)-B(10)*exp(-B(11)/R/Temp).*x(2).^B(12);
          dC = xdot;
      end
C = Cout;
end
%Plot
time_plot = [time];
c1_plot = [c1];
c2_plot = [c2];
T_o_plot = [T_o];
Cest = diff1(optimizedParams,x_in)
%figure 1, C1 and C2 vs time
figure(1)
plot(time_plot,c1_plot,'bo')
grid on
hold on
plot(time_plot,c2_plot,'ro')
plot(x_in(:,1), Cest, '--')
hold off
xlabel('Time [min]')
ylabel('Concentration [mol/L]')
legend('c1','c2')
%figure 2, C1 and C2 vs T
figure(2)
plot(T_o_plot,c1_plot,'bo')
grid on
hold on
plot(T_o_plot,c2_plot,'ro')
plot(x_in(:,2), Cest, '--')
hold off
xlabel('Temperature [K]')
ylabel('Concentration [mol/L]')
legend('c1','c2')

I shifted the plots to the end in order to incorporate the fitted equations, and added those.

It will be necessary to experiment with different initial parameter estimates to get a decent fit, since the current estimates do not create that. If you have the Golbal Optimization Toolbox, use the ga (genetic algorithm) function for this. It will generally find a reasonable fit, although the 'InitialPopulationMatrix' may need to be designed specifically for these parameters. I can hel with that if necessary.

MATLAB: Solve symbol four equations for c1 c2 c3 c4

Thanks, I solve it,

clc
clear
syms a b x y c1 c2 c3 c4 T1 T2 T3 T4  
T=c1+c2*x+c3*y+c4*x*y;
%%
r1=c1+c2*a+c3*b+c4*a*b-T1;
r2=c1-c2*a+c3*b-c4*a*b-T2;
r3=c1-c2*a-c3*b+c4*a*b-T3;
r4=c1+c2*a-c3*b-c4*a*b-T4;
%%Find C1 and C4
p13=r1+r3;
p24=r2+r4;
n1=solve(p13==0,c1);
n2=solve(p24==0,c1);
C4=solve(n1==n2,c4)
w1=solve(p13==0,c4);
w2=solve(p24==0,c4);
C1=solve(w1==w2,c1)
r1=subs(r1,c1,C1);
r1=subs(r1,c4,C4);
r2=subs(r2,c1,C1);
r2=subs(r2,c4,C4);
r3=subs(r3,c1,C1);
r3=subs(r3,c4,C4);
r4=subs(r4,c1,C1);
r4=subs(r4,c4,C4);
%%Find C3
q12=r1+r2;
q34=r3+r4;
C3=solve(q12==0,c3)
r1=subs(r1,c3,C3);
r2=subs(r2,c3,C3);
r3=subs(r3,c3,C3);
r4=subs(r4,c3,C3);
%% find C2
C2=solve(r1==0,c2)
r1=subs(r1,c2,C2);
r2=subs(r2,c2,C2);
r3=subs(r3,c2,C2);
r4=subs(r4,c2,C2);
C1
C2
C3
C4

Best Answer

Related Solutions

MATLAB: Fitting data with multiple inputs, ODE equation, and lsqnonlin

MATLAB: Solve symbol four equations for c1 c2 c3 c4

Related Question