Fitting parameters in ODE for a kinetic reaction

Hello. Ok, so I'm new to matlab and I've a question regarding parameter estimation for a kinetic model, for a plug flow reactor with 0.7m.

I have 2 different reactants and their concentrations are c1 (Oxygen) and c2(Hydrogen) with the porpuse of forming c3(Water). The initial concentration of c1 and c2 are known, being c2 constant along the reactor. The only data provided is the initial and final concentration of oxygen for the different experiments.

In this study 5 experiments were perfomed at different tempartures were (475;521;562;579;593 kelvin), and the only data available is the final concentration of oxygen in the system in function of the temperatures, that the experiment was conducted.

This means that , the initial concentration of c1 is 3971

and at the end of the experimente 1 (475k), the concentration of c1 is 3605.8

at the end of the experiment 2 (521k), the concetration of c1 is 2304.9 and then on respectively.

With the respective mass balance being determined in the following way:

    dcdx= -((S * rho_cat)/Q)* theta(1) * exp((-theta(2))/(R*Te))*c(1)*c(3) ;

we want to estimate the different "theta" .

We tried to use the answer provided by Star Stride in Monod kinetics and curve fitting, however due to our data being only the final concentration of oxygen in function of the temperature and not the lenght of the reactor, we couldnt solve it.

Hope someone can help us, i leave the code attached the question

Answer 1

[Star Strider]

My parameter estimation code will only work for dynamic integrations with respect to time.

If you want to estimate parameters based on the final time, integrate the differential equations and then use one of the parameter estimation functions (probably lsqcurvefit) to estimate the parameters.

The integration:

syms c1(t) c2(t) S rho_cat Q a1 a2 a3 a4 c10 c20 R T t 
Eq1 = diff(c1) == -((S * rho_cat) / Q) * a1* exp(-a2 / (R * T)) * c1 * 77;
Eq2 = diff(c2) == -((S * rho_cat) / Q) * a3 * exp(-a4 / (R * T)) * c2 * (77.^(1.5)) ;
[C1,C2] = dsolve(Eq1,Eq2, c1(0) == c10, c2(0) == c20)
C1C2 = matlabFunction([C1;C2], 'Vars',{[a1,a2,a3,a4],t,c10,c20,S,rho_cat,Q, R, T})

producing:

C1 =
c10*exp(-(77*S*a1*rho_cat*t*exp(-a2/(R*T)))/Q)

C2 =
c20*exp(-(5943276032390893*S*a3*rho_cat*t*exp(-a4/(R*T)))/(8796093022208*Q))

C1C2 = @(in1,t,c10,c20,S,rho_cat,Q,R,T)[c10.*exp((S.*in1(:,1).*rho_cat.*t.*exp(-in1(:,2)./(R.*T)).*-7.7e+1)./Q);c20.*exp((S.*in1(:,3).*rho_cat.*t.*exp(-in1(:,4)./(R.*T)).*(-6.756722578291934e+2))./Q)]

Note that ‘in1’ is the parameter vector of the ‘a’ values (in order). You need to create another anonymous function to put ‘C1C2’ into a form that lsqcurvefit can use, likely something like:

@(in1,T) C1C2(in1,t,c10,c20,S,rho_cat,Q,R,T)

Remember to include the initial conditions ‘c10’ and ‘c20’ and ‘t’ (probalbly the end time), as well as the rest of the parameters.

I leave the rest to you.

[Hetune Manmohandas]

Hey Star Strider

Thanks a lot for your answer !

I just didnt understand the anonymous function. The we creat a first anonymous function (C1C2) and compact the variables (a1:a4) into "in1" .

Then we create another anonymous function, to call in the lsqcurvefit

fitt = @(in1,T) C1C2(in1,t,c10,c20,S,rho_cat,Q,R,T)

After that we add our constants: t, c10, c20, S, rho_cat, Q, R

And the x-data and y-data for the curve fitting and the respective command for the lsqcurvefit, can it be due the way we provide the C1 and T? We tried to change it but nothing worked so far

C1=[3605.8 5318.6 
    2304.9 5318 
    0 4514.5 
    0 3936.2 
    0 2973 ]; 

T = [475
    521 
    562
    579
    593];

%guess
in[1;1;1;1]


[in1,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(fitt,in,T,C1);

[Star Strider]

As always, my pleasure!

There were originally two columns of ‘c’ data, your original code provided for both of them, and ‘C1C2’ regressed both of them.

If you are only fitting ‘C1’, then you can only estimate ‘a1’ and ‘a2’. The integration is now:

syms c1(t) c2(t) S rho_cat Q a1 a2 a3 a4 c10 c20 R T t 
Eq1 = diff(c1) == -((S * rho_cat) / Q) * a1* exp(-a2 / (R * T)) * c1 * 77;
C1 = dsolve(Eq1, c1(0) == c10)
C1fcn = matlabFunction(C1, 'Vars',{[a1,a2,a3,a4],t,c10,S,rho_cat,Q, R, T})

The integrated expression for ‘C1’ is:

C1 =
c10*exp(-(77*S*a1*rho_cat*t*exp(-a2/(R*T)))/Q)

and the anonymous function to estimate the parameters of ‘C1’ is:

C1fcn =  @(in1,t,c10,S,rho_cat,Q,R,T) c10.*exp((S.*in1(:,1).*rho_cat.*t.*exp(-in1(:,2)./(R.*T)).*-7.7e+1)./Q);

In the lsqcurvefit call, this anonymous function appears as:

@(in1,T) (in1,t,c10,S,rho_cat,Q,R,T)

since the objective is to estimate the parameters as a function of ‘C1’ and ‘T’. That the parameter vector is the first argument and the independent variable is the second argument is all that lsqcurvefit needs to know. The ‘C1fcn’ function will gather the other variables from its workspace. This has to do with what lsqcurvefit expects the function to accept as the arguments it provides in order to do thed optimisation, and their order in the argument list.

[Hetune Manmohandas]

If in the lsqcurvefit call, the anonymous function (C1fcn) appear as:

@(in1,T) (in1,t,c10,S,rho_cat,Q,R,T)

Then we should be able to use lsqcurvefit with only this command: (because the command already read it as inf in function of T)

C1=[3605.8  
    2304.9  
    0 
    0  
    0]; 

T = [475
    521 
    562
    579
    593];
%guess
in=[1,1]

[in1,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(C1fcn,in,T,C1);

However I obtain this:

Not enough input arguments.

Error in symengine>@(in1,t,c10,S,rho_cat,Q,R,T)c10.*exp((S.*in1(:,1).*rho_cat.*t.*exp(-in1(:,2)./(R.*T)).*-7.7e1)./Q)

Error in lsqcurvefit (line 213)
            initVals.F = feval(funfcn_x_xdata{3},xCurrent,XDATA,varargin{:});

Error in starsolo (line 39)
[in1,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(C1fcn,in,T,C1);

Caused by:
    Failure in initial objective function evaluation. LSQCURVEFIT cannot continue.

Sorry im kind of a newby in matlab, but from what i understand from this command is :

    x = lsqcurvefit(fun,x0,xdata,ydata)´
    
    %fun = C1fcn
    %x0 = in (random sample for the parameteres to be estimated)
    %xdata = temperature ( independent varible)
    %ydata = C1 ( final concentration)

sorry, i been stuck in this deadpoint allday

[Star Strider]

There was a typographical error in my previous Comment that I did not catch before I posted it.

The anonymous function to use with lsqcurvefit should be:

@(in1,T) C1fcn(in1,t,c10,S,rho_cat,Q,R,T)

.

[Hetune Manmohandas]

Wow Star cant believe i missed that!!

The code is running perfectly! Thanks a lot for your help, now i will try to do how you mention first, having both Equation in one code and calcultating all the a(x) !!

Thanks a lot! you really make understanding and use matlab easier for us (newbies)

Have a great day, thanks a lot <3

(i will still bother you if it doent work, if its okey :P)

[Star Strider]

‘Wow Star cant believe i missed that!!’

That makes two of us!

‘The code is running perfectly! Thanks a lot for your help, now i will try to do how you mention first, having both Equation in one code and calcultating all the a(x) !!’

The ‘C1C2’ function should work as written. I did not have all the necessary parameters, so I could not test it. If you have problems with it, please provide all of the parameters and the data, including the initial conditions and the value for ‘t’.

‘Thanks a lot! you really make understanding and use matlab easier for us (newbies)

Have a great day, thanks a lot <3’

Thank you! As always, my pleasure!

‘(i will still bother you if it doent work, if its okey :P)’

Of course it is!

.

[Hetune Manmohandas]

Hey Star Strider, sorry to bother you but if you are available can you have a look in the code that i attached. (it has all the variable, concetration and t)

Everything runs except the lsqcurvefit

Im having douth in how to present the respective values for the c1 and c2. since the curve fit only accepts one variable for x-data and y-data. I decided to put the values of c1 and c2into one matrix, where the first column is c1 data and second column is the c2 data

c=[3605.8            5318.6 
    2304.9           5318 
    0                4514.5 
    0                3936.2 
    0                2973 ]; 

I dont think the lsqcurvefit is reading the data in this way

[in1,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(myfunct,in,T,c,[],[],options);

Error using lsqcurvefit (line 258)
LSQCURVEFIT requires all values returned by functions to be of
data type double.

[Star Strider]

The ‘C1C2’ function is set up to produce row vectors, not column vectors, so it is necessary to transpose ‘c’ and ‘T’. "(This is easier than editing ‘C1C2’.)

This ran without error and produced what appear to be good results when I ran it (R2020a). Note the order of the statements has changed, so that ‘myfunct’ appears first after all the constants have been assigned:

C1C2 = @(in1,t,c10,c20,S,rho_cat,Q,R,T)[c10.*exp((S.*in1(:,1).*rho_cat.*t.*exp(-in1(:,2)./(R.*T)).*-7.7e+1)./Q);c20.*exp((S.*in1(:,3).*rho_cat.*t.*exp(-in1(:,4)./(R.*T)).*(-6.756722578291934e+2))./Q)]

c10=3971.5;
c20=5318.8;
t=0.7;
d_reactor = 0.148; % m (algae campaign data)
L_reactor = 0.7; % m (algae campaign data)
m_cat = 0.0486; % kg (algae campaign data)
R = 8.314 ; % J/mol*K
Q = 6.9*10^(-9); % m^3/s    based on TEPE2 report (103.68 mL / h)
S = (d_reactor^2 / 4)*pi;
rho_cat = m_cat / (S * L_reactor);

c=[3605.8 5318.6 
    2304.9 5318 
    0 4514.5 
    0 3936.2 
    0 2973 ]; 


T = [475 
    521 
    562 
    579 
    593 ];

T = T.';
c = c.';
%guess
in=[0.00001,100000,0.0001,100000]

myfunct = @(in1,T) C1C2(in1,t,c10,c20,S,rho_cat,Q,R,T);

% Q1 = myfunct(in,T);

options = optimset('Display','iter','TolX', 1e-10, 'TolFun', 1e-10, 'MaxFunEvals', 4000, 'MaxIter', 4000);
[in1,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(myfunct,in,T,c,[],[],options);

figure
hp = plot(T,c,'ko', T,myfunct(in1,T),'b-');
legend([hp(1) hp(3)], 'real', 'fitting')
fprintf(1,'\tRate Constants:\n')

for k1 = 1:length(in1)
    fprintf(1, '\t\ta(%d) = %14.5f\n', k1, in1(k1))
end

producing:

	Rate Constants:
		a(1) =       94.80747
		a(2) =   109265.44454
		a(3) =        7.84732
		a(4) =   122732.35669

and a good curve fit. I changed the format in the loop slightly, to make it more readable, and updated the legend arguments so they will show what you want them to show.

This produces a slightly smoother curve:

Tv = linspace(min(T), max(T), 50);

figure
hp = plot(T,c,'ko', Tv,myfunct(in1,Tv),'b-');
legend([hp(1) hp(3)], 'real', 'fitting')
fprintf(1,'\tRate Constants:\n')

.

[Shailesh Pathak]

Dear Strider

I am getting imaginary number while solving the ode for the estimation of parameter,

Can you please help me? I shall be very thankful to you.

k(1)=1.00955;

k(2)=0.6488;

k(3)=1.88562;

t=[0;1.8;2.8;3.2;3.8;4.4]; % x axis values

cfit3 = kinetics(k,t)

function C=kinetics(k,t)

c0=[16.8;150];

[~,Cv]=ode113(@DifEq,t,c0);

function dC=DifEq(t,c)

rate=zeros(2,1);

rate(1,1)=-7.8*c(1)*(c(2))/((1+(k(2)*c(1))+(k(3)*((c(2))^0.5))));

rate(2,1)=-12.8*c(1)*(c(2))/((1+(k(2)*c(1))+(k(3)*((c(2))^0.5))));

dC=rate;

end

C=Cv;

end

[Torsten]

The kinetic parameters seem to be such that c(2) becomes negative, and the root of a negative number gives a complex number.

[Shailesh Pathak]

Thanks Torsten. Is there any way to avoid these imaginary numbers?

[Torsten]

Try ode15s instead of ode113 and set RelTol and AbsTol to small numbers, e.g 1e-8.

[Star Strider]

@Torsten — Thank you!

[Shailesh Pathak]

Thanks Torsten. It works for me.

Answer 2

[Hetune Manmohandas]

Thanks a lot Star Strider! you were a great help, eternally gratefull

Best Regards

H

[Star Strider]

As always, my pleasure!

This was an interesting problem!