Change domain of kinetic functions in order to fit with ParametricNDSolve

Question

I'm trying to solve kinetic differential equation depended on time with three parameters: preexponential factor a, reaction order n and temperature T. My experimental data depends on temperature, not time, so my approach is: solve ParametricNDSolve with few temperature values, then collect values for a certain time t for all teperatures and then fit result to experiment. When I assume that preexp factor and reaction order as numbers, it works well:

temps = Table[T, {T, 300, 600, 25}]; (*temperature range*)
Ea = 90000;  (*activation energy in J/mol*)
R = 8.31;  (*gas constant*)

sol = ParametricNDSolve[{c'[t] == -a*Exp[-Ea/(R*#)]*c[t]^n, 
     c[0] == 1}, c, {t, 0, 10}, {a, n}] & /@ temps;

c1 = c[10^9, 1.5] /. sol; (*assuming a=10^9 and n=1.5*)
c2 = c1[[#]][10] /. sol[[#]] & /@ Range[Length[temps]]; (*collect values for time=10 at all temperatures*)

conv = (c2[[1]] - c2[[#]])/c2[[1]] & /@Range[Length[c2]]; (*construct conversion values at all temperatures*)
data = {temps, conv} // Transpose;
noisedata = {temps,data[[#, 2]] + RandomReal[{-0.02, 0.02}] & /@Range[Length[conv]]} // Transpose;

My problem is: how to construct the code to fit parameters a and n to noisedata ? I've tried the same approach:

c3 = c[a, n] /. sol;
c4 = c3[[#]][10] /. sol[[#]] & /@ Range[Length[temps]];
conv2 = (c4[[1]] - c4[[#]])/c4[[1]] & /@ Range[Length[c4]];
data2 = {temps, conv2} // Transpose;
f = Interpolation[data2] (* interpolate to obrain continous function *)

Interpolation with parameters as before plot the correct answer:

Show[Plot[f[T] /. {a -> 10^9, n -> 1.5}, {T, 300, 600}], 
 ListPlot[data]]

But it fails when it comes to fit:

fit = NonlinearModelFit[noisedata, f[T], {{a, 10^9}, {n, 1.4}}, T]

(* NonlinearModelFit::fmgz: Encountered a gradient that is effectively zero. The result returned
 may not be a minimum; it may be a maximum or a saddle point.

or

NonlinearModelFit::nrlnum: The function value {0.0123183 +0. I,(...)} is not a list of real numbers with
 dimensions {13} at {a,n} = {1.*10^13,1.4}. *)

Which is not true, because applying these number to Interpolate gives real results. Aby other possibilities to solve that puzzle?

Answer 1

I think the main problem here is that you are trying to fit to an interpolating function, which no longer explicitly depends on your parameters, instead of to a model function containing those parameters.

Below I've refactored your code rather significantly because I think you were making your life harder than it needs to be :-) The main change is the use of ParametricNDSolveValue to directly calculate the conversion rate at $t=10$, rather than doing that by hand afterwards. This also provides a parametric function object that can be used as a model in your fitting.

Let's start by generating some exact solutions, as you did:

Ea = 90000; R = 8.31;
c10 = 
  ParametricNDSolveValue[
    {c'[t] == -a Exp[-Ea/(R T)]c[t]^n, c[0] == 1},
    (c[0] - c[10])/c[0], {t, 0, 15}, {a, n, T}
  ]

Plot[c10[10^9, 1.5, T], {T, 300, 600}]

Let's fabricate some noisy data now:

noisy=
  Table[
    {T, c10[10^9, 1.5, T] + RandomReal[{-0.03, 0.03}]},
    {T, 300, 600, 25}
  ];

ListPlot[noisy, Frame -> True]

Now let's fit the noisy data to the model provided by c10[a, n, T], where $a$ and $E_a$ are parameters to be fit, and $T$ is the independent variable:

fit = NonlinearModelFit[noisy, c10[a, n, T], {{a, 5 10^7}, {n, 4}}, T];
fit["BestFitParameters"]

(* Out: {a -> 1.02589*^9, n -> 1.53829} *)

Above I've used relatively "bad" starting values to show that the fit is somewhat resilient as well, at least with this relatively low amount of noise. Finally let's compare the fitted model with the noisy data it was fitted to, to show that it is a good fit:

Plot[fit[T], {T, 300, 600}, Epilog -> Point[noisy], Frame -> True]