1
$\begingroup$

I am looking for some help estimating the parameters of a series of differential equations to assist fitting the experimental data I have collected.

The Problem:
The equations are based off of the work of fellow researchers: awesome equations (PREDICTIONS OF CONCENTRATION IN THE PYROLYSIS OF BIOMASS MATERIALS--I V. K. SRIVASTAVA and R. K. JALAN, 1994)

My fitting to my experimental data been pretty bad thus far using MATLAB. bad bad fit :( If you want my MATLAB code as well I would be happy to share. Reading around online leads me to believe this is mainly due to improper scaling of data in addition to bad initial guesses. Therefore, this is my attempt to better approximate the initial guesses for the parameters using Mathematicas manipulate function.

My approach:
I first started out defining constants and defining equation (6) from the research paper above. Rg is ideal gas constant ($kJ*mol^{-1}*K^{-1}$), HR is heating rate ($K*min^{-1}$), n are the reaction rates (n1=0,n2=n3=1.5 from literature), T is temperature ($K$).

Rg = 8.314*10^-3;
HR = 5; 
n1 = 0;
n2 = 1.5;
n3 = 1.5;
T = HR*t + T0;

I then defined the system of differential equations in the research paper above (1),(2),(3),(4) and the initial conditions (5) as so:

system = {
   Cb'[t] == -a1*exp (-e1/(Rg*T))*(Cb[t]^n1), 
   Cbp'[t] == 
    a1*exp (-e1/(Rg*T))*Cb[t]^n1 - a2*exp (-e2/(Rg*T))*Cbp[t]^n2 - 
     a3*exp (-e3/(Rg*T))*Cbp[t]^n3, 
   Cg'[t]  == a2*exp (-e2/(Rg*T))*Cbp[t]^n2,
   Cc'[t]  == a3*exp (-g/(Rg*T))*Cbp[t]^n3,
   Cb[0] == 1, Cbp[0] == 0, Cg[0] == 0, Cc[0] == 0
   };

Finally, I attempt to find the numerical solution to the set of ODE's by using Mathematicas ParametricNDSolve. I set the time interval to be from 0 --> 120 minutes, define my parameters and try to manipulate the equation parameters to find something close to resembling my data. I left the initial temperature (T0) as a parameter, though the reaction from experimental data began at around 325K. At a constant heating rate of 5K/min, the time at that point would be ~65 min into heating. 

parfun = ParametricNDSolveValue[
   system, {Cb, Cbp, Cg, Cc}, {t, 0, 120}, {a1, e1, a2, e2, a3, e3,
     T0}];

Manipulate[
 Plot[
  Evaluate@Through[parfun[a1, e1, a2, e2, a3, e3, T0][t]], {t, 0, 
   120}, PlotLegends -> {"Cb[t]", "Cb+[t]", "Cg[t]", "Cc[t]"}
  ],
 {{a1, 10}, 0.000001, 1000, Appearance -> "Labeled"}, {{e1, 200}, 
  0.000001, 5000, Appearance -> "Labeled"}, {{a2, 10}, 0.000001, 1000,
   Appearance -> "Labeled"}, {{e2, 1}, 0.000001, 5000, 
  Appearance -> "Labeled"}, {{a3, 10}, 0.000001, 1000, 
  Appearance -> "Labeled"}, {{e3, 1}, 0.000001, 5000, 
  Appearance -> "Labeled"}, {{T0, 325}, 270, 1000, 
  Appearance -> "Labeled"}
 ]

My data is represented as Temperature in Kelvin in the first column and Equation (1) in the research paper for the second column. 

data = {{328.284, 0.964936}, {333.399, 0.958594}, {338.515, 
  0.951508}, {343.63, 0.945054}, {348.745, 0.93767}, {353.86, 
  0.930621}, {358.975, 0.924241}, {364.091, 0.918532}, {369.206, 
  0.913828}, {374.321, 0.909571}, {379.436, 0.906699}, {384.552, 
  0.902843}, {389.667, 0.899256}, {394.356, 0.895557}, {399.897, 
  0.893198}, {405.013, 0.890616}, {410.128, 0.888927}, {415.243, 
  0.886568}, {420.358, 0.883986}, {425.473, 0.881068}, {430.589, 
  0.879044}, {435.704, 0.875904}, {440.819, 0.873992}, {445.934, 
  0.872079}, {451.05, 0.868269}, {456.165, 0.864346}, {461.28, 
  0.859754}, {466.395, 0.854045}, {471.511, 0.847107}, {476.626, 
  0.839836}, {481.741, 0.83232}, {486.856, 0.823817}, {491.545, 
  0.816521}, {495.808, 0.808806}, {499.644, 0.800254}, {503.054, 
  0.793412}, {506.465, 0.784241}, {509.875, 0.776555}, {512.859, 
  0.768897}, {515.842, 0.760681}, {518.826, 0.752716}, {521.81, 
  0.744918}, {524.794, 0.736366}, {527.778, 0.727981}, {531.188, 
  0.719429}, {534.598, 0.711211}, {538.008, 0.702269}, {541.419, 
  0.694721}, {544.829, 0.687876}, {548.665, 0.680476}, {551.649, 
  0.672473}, {555.205, 0.664087}, {557.069, 0.655255}, {561.027, 
  0.648548}, {564.863, 0.639606}, {568.274, 0.630329}, {570.831, 
  0.622772}, {573.815, 0.614992}, {576.799, 0.607214}, {579.868, 
  0.598194}, {582.426, 0.58951}, {584.642, 0.580256}, {587.456, 
  0.570534}, {589.928, 0.561144}, {592.656, 0.55343}, {595.129, 
  0.545915}, {598.539, 0.538331}, {601.949, 0.532287}, {605.785, 
  0.526689}, {610.474, 0.52098}, {615.589, 0.515829}, {620.705, 
  0.510343}, {625.82, 0.504746}, {630.935, 0.500154}, {636.05, 
  0.496231}, {641.166, 0.492197}, {646.281, 0.48794}, {651.396, 
  0.484241}, {656.511, 0.480095}, {661.627, 0.476285}, {666.742, 
  0.473032}, {671.857, 0.469446}, {676.972, 0.463736}, {682.087, 
  0.458808}, {686.947, 0.454328}, {692.318, 0.449735}, {697.433, 
  0.445366}, {702.548, 0.440997}, {707.664, 0.43607}, {712.779, 
  0.431478}, {717.894, 0.426662}, {723.009, 0.421958}, {728.124, 
  0.4177}, {733.24, 0.414785}, {738.355, 0.411864}, {743.47, 
  0.409951}, {747.875, 0.407258}, {754.553, 0.404675}, {759.668, 
  0.40254}, {764.784, 0.400404}, {769.899, 0.398715}, {775.014, 
  0.396133}, {780.129, 0.394891}, {785.245, 0.392644}, {790.36, 
  0.390173}, {795.475, 0.388038}, {800.59, 0.386461}, {805.705, 
  0.384437}, {810.821, 0.382748}, {815.936, 0.379831}, {821.051, 
  0.378923}, {826.166, 0.376788}, {831.282, 0.374541}, {836.397, 
  0.372293}, {841.512, 0.371051}, {846.627, 0.369251}, {851.743, 
  0.36812}, {856.858, 0.367102}, {861.973, 0.365189}, {867.088, 
  0.363054}, {872.203, 0.361253}, {877.319, 0.359229}, {882.434, 
  0.357094}, {887.549, 0.35574}, {892.664, 0.353046}, {897.78, 
  0.351246}, {902.895, 0.34911}, {908.01, 0.347399}, {913.125, 
  0.343313}, {918.24, 0.340358}, {923.441, 0.334761}, {928.471, 
  0.330393}, {933.586, 0.325219}};
pdata = ListPlot[data, PlotStyle -> {PointSize[Medium], Red}];
Show[pdata, 
 AxesLabel -> {HoldForm[Temperature K], HoldForm[Concentration]}, 
 PlotLabel -> None, LabelStyle -> {GrayLevel[0]}]

I keep getting issues trying to estimate the parameters. In particular: 

Encountered non-numerical value for a derivative at t == 0.`

I thought originally it was because the timestep was starting at zero in the numerical approximation, so I changed the times to start at 0.000001, as starting at t=0 would would lead to division by zero if the T0 was zero along with t, however this did not fix my issues.

I think my code requires a relatively simple fix, if someone who uses Mathematica frequently could take a look I don't think it would take too long to fix. Thank you!

Improve this question
$\endgroup$
2
  • 1
    $\begingroup$ The Mathematica function for exponential is Exp[] (case sensitive) or simply E^(). You will want to replace "exp" appropriately. I still had issues after the replacement. $\endgroup$
    – Tim Laska
    May 8, 2019 at 0:36
  • $\begingroup$ Also, the activation energy for Cc'[t] equation is "g". It probably should be "e3". $\endgroup$
    – Tim Laska
    May 8, 2019 at 0:49

1 Answer 1

Reset to default
2
$\begingroup$

I corrected the typos leading to following Mathematica code:

Rg = 8.314*10^-3;
HR = 5;
n1 = 0;
n2 = 1.5;
n3 = 1.5;
T = HR*t + T0;
system = {Cb'[t] == -a1*E^(-e1/(Rg*T))*(Cb[t]^n1), 
   Cbp'[t] == 
    a1*E^ (-e1/(Rg*T))*Cb[t]^n1 - a2*E^ (-e2/(Rg*T))*Cbp[t]^n2 - 
     a3*E^ (-e3/(Rg*T))*Cbp[t]^n3, 
   Cg'[t] == a2*E^ (-e2/(Rg*T))*Cbp[t]^n2, 
   Cc'[t] == a3*E^ (-e3/(Rg*T))*Cbp[t]^n3, Cb[0] == 1, Cbp[0] == 0, 
   Cg[0] == 0, Cc[0] == 0};
parfun = ParametricNDSolveValue[
   system, {Cb, Cbp, Cg, Cc}, {t, 0, 120}, {a1, e1, a2, e2, a3, e3, 
    T0}];
Manipulate[
 Plot[Evaluate@Through[parfun[a1, e1, a2, e2, a3, e3, T0][t]], {t, 0, 
   120}, PlotLegends -> {"Cb[t]", "Cb+[t]", "Cg[t]", "Cc[t]"}], {{a1, 
   8}, 0.000001, 1000, Appearance -> "Labeled"}, {{e1, 40}, 0.000001, 
  100, Appearance -> "Labeled"}, {{a2, 294}, 0.000001, 500, 
  Appearance -> "Labeled"}, {{e2, 17}, 0.000001, 100, 
  Appearance -> "Labeled"}, {{a3, 486}, 0.000001, 1000, 
  Appearance -> "Labeled"}, {{e3, 32}, 0.000001, 100, 
  Appearance -> "Labeled"}, {{T0, 325}, 270, 1000, 
  Appearance -> "Labeled"}]

TGA Plot

There are problems with your reaction mechanism, because the zero order kinetics for Cb cause it to go negative.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Tim, thank you very much for taking the time to fix my mistakes! Seems to be working fine now. $\endgroup$
    – user3741793
    May 8, 2019 at 5:59
  • $\begingroup$ @user3741793 You are welcome. Good luck with your fits! $\endgroup$
    – Tim Laska
    May 8, 2019 at 11:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.