# Estimating parameters on system of differential equations

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: (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. 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!

• 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. Commented May 8, 2019 at 0:36
• Also, the activation energy for Cc'[t] equation is "g". It probably should be "e3". Commented May 8, 2019 at 0:49

## 1 Answer

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"}]
`

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

• Tim, thank you very much for taking the time to fix my mistakes! Seems to be working fine now. Commented May 8, 2019 at 5:59
• @user3741793 You are welcome. Good luck with your fits! Commented May 8, 2019 at 11:41