How to use FindFit to do kinetic simulation?

Question

I'd like to simulate the chemical kinetic model by fitting the experimental data. Here are my codes for simulation.

data = {{0.`, 0.003207613963640383`}, {0.256348`, 
0.751307591324145`}, {0.5004879999999998`, 
0.946609520707414`}, {0.7568359999999998`, 
0.9997534786329835`}, {0.9887700000000001`, 
0.9788588889284163`}, {2.0019530000000003`, 
0.7107228417712506`}, {3.00293`, 0.4859859783085131`}, {4.003906`,
 0.329989302539268`}, {5.004883`, 
0.2236874613858439`}, {6.005858999999999`, 
0.15155994391149338`}, {7.006836`, 
0.10294001208237556`}, {8.007815`, 
0.07060326258953398`}, {9.008785`, 
0.048587435830738386`}, {10.009765`, 0.033384617610655296`}};

Clear[knn, kf, knf, a];
ttime = 10.253905;
knf = 0.4;
a = 2.36;
model = ParametricNDSolveValue[{
S1'[t] == (-kf*S1[t] - knn*S1[t]),
NN'[t] == (knn*S1[t] - knf*NN[t]),
S[t] == a*NN[t],
S1[0] == 1, NN[0] == 0}, S, {t, 0, ttime}, {kf, knn}]

fit = FindFit[data, model[kf, knn][t], {{kf, 1.5}, {knn, 2}}, t]
Plot[model[kf, knn][t] /. fit, {t, 0, ttime},Epilog -> {Red, Point[data]}]

I have tried the condition that kf = 1.5, knn = 2, a = 2.36 can roughly fit the data by using manually change the parameter. However, I get an error like this.

ParametricNDSolveValue:: Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions

I don't know how to fix this problem. I hope if I take {kf, knn, a} these 3 parameters to freely fitting with suitable initial value, I can also correctly fit the data and the rate constants are physically reasonable.(0 < kf, knn < 10, a has no constrained )

I would appreciate any help in rebuilding the code.

Many Thanks!

Answer 1

A quick "fix" would be to remove S from the ParametricNDSolveValue command:

Clear[knn, kf, knf, a];
ttime = 10.253905;
knf = 0.4;
a = 2.36;
model = ParametricNDSolveValue[{S1'[t] == (-kf*S1[t] - knn*S1[t]), 
  NN'[t] == (knn*S1[t] - knf*NN[t]), S1[0] == 1, NN[0] == 0}, 
  NN, {t, 0, ttime}, {kf, knn}];

fit = FindFit[data, model[kf, knn][t], {{kf, 1.5}, {knn, 2}}, t]
(* {kf -> -1.17998, knn -> 4.35971} *)

Plot[model[kf, knn][t] /. fit, {t, 0, ttime}, Epilog -> {Red, Point[data]}]

A better fix would be to directly use S instead of NN, incorporate the estimation of a, and use NonlinearModelFit:

a =.;
model2 = ParametricNDSolveValue[{S1'[t] == (-kf*S1[t] - knn*S1[t]), 
    S'[t] == a (knn*S1[t] - knf S[t]/a), S1[0] == 1, S[0] == 0}, 
   S, {t, 0, ttime}, {a, kf, knn}];

fit2 = NonlinearModelFit[data, model2[a, kf, knn][t], {{a, 1}, {kf, -1.5}, {knn, 4.4}}, t];
fit2["BestFitParameters"]
(* {a -> 1.0003133640741626, kf -> -1.1786150346863764, knn -> 4.3583459462046354} *)

Plot[model2[a, kf, knn][t] /. fit2["BestFitParameters"], {t, 0, 
  ttime}, Epilog -> {Red, Point[data]}] 

With NonlinearModelFit you can get measure of precision for the estimated parameters (standard errors, confidence limits, correlations among parameters, etc.) and an assortment of model checking features. One essential feature is a look at the residuals (as you are assuming a particular error distribution):

ListPlot[Transpose[{fit2["PredictedResponse"], fit2["FitResiduals"]}],
  Frame -> True, FrameLabel -> {"Predicted", "Residual"}]

One can see a definite pattern remaining suggesting that more is going on that what the current model can describe. Also there appears to be a single data point that needs checking out.