Global data fitting for a chemical kinetics problem?

Questions
, , ,
1

Dear all
The problem of global data analysis is a common problem among experimenters, especially in chemical kinetics.
I have different kinetics (variation of absorbance as a function of time) at different wavelengths. I would like to adjust a set of 10 kinetics. These kinetics are governed by the same kinetic law, the same rate constants (here k3,k4,k5 ) but the absorption coefficients (ea, eh and ef) depend on the wavelength. I usually do these fits with Igor (Wavemetrics) but without any analytical solution I turned to Mathematica.
I tried the MultiNonlinearModelFit function

Clear[k3, k4, k5, ea1, eh1, ef1, ea2, eh2, ef2]
sol = ParametricNDSolve[{Derivative[1][ag][t] == -k3*ag[t], 
    Derivative[1][gox][t] == k3*ag[t], 
    Derivative[1][h][t] == 
     k3*ag[t] - 2*k4*h[t]*h[t] + 2*k5*ff[t]*ff[t], 
    Derivative[1][ff][t] == 2*k4*h[t]*h[t] - 2*k5*ff[t]*ff[t], 
    ag[0] == 0.000006, h[0] == 0, gox[0] == 0, ff[0] == 0}, {ag, gox, 
    h, ff}, {t, 0, 2000}, {k3, k4, k5}];
{agFunc, goxFunc, hFunc, fFunc} = {ag, gox, h, f} /. sol;
model1[k3_, k4_, k5_, ea1_, eh1_, ef1_][
   t_] := (((agFunc[k3, k4, k5][t])*ea1) + (eh1*
      hFunc[k3, k4, k5][t]) + (ef1*fFunc[k3, k4, k5][t]));
model2[k3_, k4_, k5_, ea2_, eh2_, ef2_][
   t_] := (((agFunc[k3, k4, k5][t])*ea2) + (eh2*
      hFunc[k3, k4, k5][t]) + (ef2*fFunc[k3, k4, k5][t]));
fit = MultiNonlinearModelFit[{dat1, 
    dat2}, {model1[k3, k4, k5, ea1, eh1, ef1][t], 
    model2[k3, k4, k5, ea2, eh2, ef2][t]}, {k3, k4, k5, ea1, eh1, ef1, 
    ea2, eh2, ef2}, {t}];
fit["BestFitParameters"]

I have a concern with the output I don’t understand.

ParametricNDSolve::ndsz: At t == 3.127085338387159`, step size is effectively zero; singularity or stiff system suspected.

InterpolatingFunction::dmval: Input value {3.2} lies outside the range of data in the interpolating function. Extrapolation will be used.

…

NonlinearModelFit::nrlnum: The function value {-0.039527+1. f[1.,1.,1.][0.2],-0.0396433+1. f[1.,1.,1.][0.4],-0.0396255+1. f[1.,1.,1.][0.6],-0.0394311+1. f[1.,1.,1.][0.8],-0.0393989+1. f[1.,1.,1.][1.],-0.0394361+1. f[1.,1.,1.][1.2],-0.0394022+1. f[1.,1.,1.][1.4],-0.039186+1. f[1.,1.,1.][1.6],-0.0392537+1. f[1.,1.,1.][1.8],-0.0392642+1. f[1.,1.,1.][2.],<<32>>,-0.0384155+1. f[1.,1.,1.][8.6],-0.0383621+1. f[1.,1.,1.][8.8],-0.0385297+1. f[1.,1.,1.][9.],-0.0384413+1. f[1.,1.,1.][9.2],-0.0384188+1. f[1.,1.,1.][9.4],-0.0382994+1. f[1.,1.,1.][9.6],-0.0383953+1. f[1.,1.,1.][9.8],-0.0381658+1. f[1.,1.,1.][10.],<<850>>} is not a list of real numbers with dimensions {900} at {k3,k4,k5,ea1,eh1,ef1,ea2,eh2,ef2} = {1.,1.,1.,1.,1.,1.,1.,1.,1.}.

Yet it provides me with parameters which for some are far from the expected values:
{k3 → 76620.6, k4 → -23293.5, k5 → 10237.6, ea1 → 169.232,
eh1 → 81.862, ef1 → 14461.1, ea2 → 157.347, eh2 → 80.0374,
ef2 → 14761.5}

Parts of the data are

dat2 = {{0.2, 0.0366929}, {0.4, 0.0368321}, {0.6, 0.0368786}, {0.8, 
   0.0367551}, {1, 0.036677}, {1.2, 0.036687}, {1.4, 0.0365921}, {1.6,
    0.0365438}, {1.8, 0.0364348}, {2, 0.0364835}, {2.2, 
   0.036489}, {2.4, 0.0363618}, {2.6, 0.0361545}, {2.8, 
   0.0362492}, {3, 0.0362564}, {3.2, 0.0362539}, {3.4, 
   0.0361269}, {3.6, 0.0362652}, {3.8, 0.0361286}, {4, 
   0.0361022}, {4.2, 0.0360444}, {4.4, 0.0360916}, {4.6, 
   0.0360261}, {4.8, 0.0360464}, {5, 0.0359304}, {5.2, 
   0.0358969}, {5.4, 0.0359568}, {5.6, 0.036002}, {5.8, 
   0.0358563}, {6, 0.0360431}, {6.2, 0.0358403}, {6.4, 
   0.0358596}, {6.6, 0.0358963}, {6.8, 0.0359946}, {7, 
   0.0358273}, {7.2, 0.0359167}, {7.4, 0.0357477}, {7.6, 
   0.0357798}, {7.8, 0.0358253}, {8, 0.0358452}, {8.2, 
   0.0358645}, {8.4, 0.035836}, {8.6, 0.0358447}, {8.8, 
   0.0357238}, {9, 0.0358621}, {9.2, 0.0358738}, {9.4, 
   0.0357445}, {9.6, 0.0356904}, {9.8, 0.0357929}, {10, 0.0355983}}

and

dat1 = {{0.2, 0.039533}, {0.4, 0.0396493}, {0.6, 0.0396315}, {0.8, 
   0.0394371}, {1, 0.0394049}, {1.2, 0.0394421}, {1.4, 
   0.0394082}, {1.6, 0.039192}, {1.8, 0.0392597}, {2, 
   0.0392702}, {2.2, 0.039207}, {2.4, 0.0391561}, {2.6, 
   0.0390186}, {2.8, 0.038951}, {3, 0.0390279}, {3.2, 
   0.0389623}, {3.4, 0.0389681}, {3.6, 0.0389333}, {3.8, 
   0.0387812}, {4, 0.0388705}, {4.2, 0.0388193}, {4.4, 
   0.0387575}, {4.6, 0.0386211}, {4.8, 0.0388169}, {5, 
   0.0387259}, {5.2, 0.038528}, {5.4, 0.0387643}, {5.6, 
   0.03884}, {5.8, 0.0385277}, {6, 0.0385109}, {6.2, 0.0386881}, {6.4,
    0.0385117}, {6.6, 0.038618}, {6.8, 0.0386693}, {7, 
   0.0385695}, {7.2, 0.0385362}, {7.4, 0.0384269}, {7.6, 
   0.038459}, {7.8, 0.038489}, {8, 0.0384817}, {8.2, 0.0384481}, {8.4,
    0.0384566}, {8.6, 0.0384215}, {8.8, 0.0383681}, {9, 
   0.0385357}, {9.2, 0.0384473}, {9.4, 0.0384248}, {9.6, 
   0.0383054}, {9.8, 0.0384013}, {10, 0.0381718}}

other question, but perhaps not far from the problem: is it possible to reduce the search area for the right parameters by specifying an existence interval for each parameter. E.g. 0.001<k3<0.01, 1<k4<400, 500<ea, etc?

Thank you in advance for your advice and help.

Hélène

2

What version do you have that has a function called MultiNonlinearModelFit?

Found it: https://resources.wolframcloud.com/FunctionRepository/resources/MultiNonlinearModelFit.

And you do have a typo:

{agFunc, goxFunc, hFunc, fFunc} = {ag, gox, h, f} /. sol;

should be

{agFunc, goxFunc, hFunc, fFunc} = {ag, gox, h, ff} /. sol;
3

Dear Jim
The version I used is the one I found at the link you sent me, that is to say

Options[MultiNonlinearModelFit] = {AccuracyGoal -> Automatic, 
   ConfidenceLevel -> 19/20, EvaluationMonitor -> None, 
   Gradient -> Automatic, MaxIterations -> Automatic,
   	Method -> Automatic, PrecisionGoal -> Automatic, 
   StepMonitor -> None, Tolerance -> Automatic, 
   VarianceEstimatorFunction -> Automatic,
   	Weights -> Automatic, WorkingPrecision -> Automatic, 
   "DatasetIndexSymbol" -> \[FormalN]
   }; (* \[Equal] Join[Options[NonlinearModelFit], \
{"DatasetIndexSymbol" -> \[FormalN]}] *)

MultiNonlinearModelFit[datasets_, form_, fitParams_, 
   independents : Except[_List], opts : OptionsPattern[]] := 
      MultiNonlinearModelFit[datasets, form, 
   fitParams, {independents}, opts];
 
MultiNonlinearModelFit[datasets_, form : Except[{__Rule}, _List], 
   fitParams_, independents_, opts : OptionsPattern[]] := 
      MultiNonlinearModelFit[
   datasets, <|"Expressions" -> form, "Constraints" -> True|>, 
   fitParams, independents, opts];
 
MultiNonlinearModelFit[
       datasets : {__?(MatrixQ[#1, NumericQ] &)}, 
       KeyValuePattern[{"Expressions" -> expressions_List, 
     "Constraints" -> constraints_}], 
       fitParams_List, 
       independents_List,
       opts : OptionsPattern[]
   ] := Module[{
        fitfun, weights,
        numSets = Length[expressions], 
        augmentedData = Join @@ MapIndexed[
               Join[ConstantArray[N[#2], Length[#1]], #1, 2] &,
               datasets
           ], 
        indexSymbol = OptionValue["DatasetIndexSymbol"]
    },
       Condition[
            fitfun = With[{
                   
       conditions = 
        Join @@ ({#1, expressions[[#1]]} & ) /@ Range[numSets]
               }, 
                  Switch @@ Prepend[conditions, Round[indexSymbol]]
              ]; 
            weights = Replace[
                  OptionValue[Weights],
                  {
                       (list_List)?(VectorQ[#1, NumericQ] & ) /; 
         Length[list] === numSets :> 
                            
        Join @@ MapThread[ConstantArray, {list, Length /@ datasets}], 
                       
       list : {__?(VectorQ[#1, NumericQ] & )} /; 
         Length /@ list === Length /@ datasets :>
                            Join @@ list, 
                       
       "InverseLengthWeights" :> 
        Join @@ (ConstantArray[1./#1, #1] & ) /@ Length /@ datasets
                   }
              ]; 
            NonlinearModelFit[
                 augmentedData,
                 
     If[TrueQ[constraints], fitfun, {fitfun, constraints}], 
                 fitParams,
                 Flatten[{indexSymbol, independents}],
                 Weights -> weights, 
                 
     Sequence @@ FilterRules[{opts}, Options[NonlinearModelFit]]
             ]
            ,
            numSets === Length[datasets]
        ]
   ];

Thanks for the typo, I had forgotten to correct it in my message.
Hélène

4

I think the main issue is that both models are overparameterized for the available data (and secondarily better starting values are needed). A plot of the data shows just slight curvature over values of t for both datasets suggesting that maybe just 3 parameters for each curve might be warranted. Or at least that having 6 parameters per curve is too much for the available data.

I want to emphasize that I’m not saying the models are wrong. It’s just that the available data (meaning values of t from 0 to 10) doesn’t allow for the estimation of all of the parameters. Maybe if there are more than two datasets, then more parameters might be able to be estimated.

Here is what works:

fit = MultiNonlinearModelFit[{dat1, dat2}, {model1[k3, k4, k5, ea1, eh1, ef1][t], 
    model2[k3, k4, k5, ea2, eh2, ef2][t]},
   {{k3, 0.22}, {k4, 0}, {k5, 0}, {ea1, 6600}, {eh1, 6400}, {ef1, 0}, {ea2, 6200}, {eh2, 6000}, {ef2, 0}}, {t}];
mle = fit["BestFitParameters"]
(* {k3 -> 0.2188, k4 -> -5.24086, k5 -> 0., ea1 -> 6623.65, 
 eh1 -> 6364.94, ef1 -> 0.000336612, ea2 -> 6152.56, eh2 -> 5927.64, ef2 -> 0.000419775} *)
fit["CorrelationMatrix"] // MatrixForm

correlationmatrix.png
correlationmatrix.png721×155 5.46 KB

fit["ParameterTable"]

parametertable2.png
parametertable2.png354×191 4.6 KB

Note that the P-values suggest that k4, k5, ef1, and ef2 are not statistically significantly different from zero. However, the huge standard errors suggest that those parameters are not estimated very well. In short, there’s not much one can say about those 4 parameters.

Show[ListPlot[{dat1, dat2}],
 Plot[{model1[k3, 0, 0, ea1, eh1, 0][t] /. mle,
   model2[k3, 0, 0, ea2, eh2, 0][t] /. mle}, {t, 0, 10}]]

twoplots.png
twoplots.png367×214 3.57 KB

While MultiNonlinearModelFit is very handy for multiple datasets with common parameters, it does assume that the variability about the individual curves is the same (although with the Weight option one could allow for knowing that the variability is, say, twice as much in one dataset as another but I suspect that is a rare occurrence.) I think an explicit statement in the documentation about the assumption of the equality of the variances is warranted.

5

Indeed, you’re right, I only sent you the beginning of the data : the data goes up to 850 and not 10.

I managed to do again what you did. Thank you so much for that breakthrough.I guess {k3, 0.22} means that the k3 values start from 0.22. Does it ?

Hélène

6

Hi Jim,
first of all thank you for your kind attention. Now I succeed. I’ve simplified my model a bit. K5 was not very useful ! I still have to implement this with my 10 kinetics…

Once again, many thanks!
Hélène