Message Boards Message Boards

0
|
7714 Views
|
5 Replies
|
1 Total Likes
View groups...
Share
Share this post:

Fitting model with experimental data from a kinetic (batch culture)

Posted 8 years ago

Hi every body, I am trying to fit a model (batch culture) with experimental data (biomass ,subtrate uptake and ethanol concentration) but I can´t get it. Actually the model is a nonlinear ordinary diffrential equation system with 4 equations. I guess that the problem is the way that I am importing the data. I have doubts about the terms: a)coordinates, (when the answers is : number of coordinates is not equal to number of variables) b)list of real numbers with dimensions {8} at {Ke,YxsG,ms,YpsG,[Alpha]a,[Beta]b},

So, I would like to someone can help me. If is possible, of course.

Best Regards. Ricardo Sánchez

Attachments:
POSTED BY: Ricardo Sánchez
5 Replies
Posted 8 years ago

(1) There really is no reason to attach a notebook. The code is short enough that it can go directly into a post.

(2) There is no evidence of any debugging effort whatever. The NDSolve result, for example, will not evaluate as utilized inside FindFit.

(3) There is only one independent variable, and it is t. As there are three simultaneous fits needed, we emulate that with a discrete variable n that takes on values 1,2,3. For this to work we need to rearrange the date list a bit.

Cinetica2 = {{0, 0.04`, 16.04`, 0.3`}, {2, 0.19`, 13.52`, 0.51`}, {4, 0.39`, 10.99`, 1.22`}, {6, 1.11`, 8.23`, 2.96`}, 
  {8, 1.96`, 0.37`,  6.08`}, {10, 2.09`, 0.31`, 5.82`}, {12, 2.41`, 0.28`,  5.69`}, {14, 2.41`, 0.26`, 5.62`}};

Rearrange the data:

clists = Flatten[Map[Thread[{First[#], Range[3], Rest[#]}] &, Cinetica2], 1]

(* Out[165]= {{0, 1, 0.04}, {0, 2, 16.04}, {0, 3, 0.3}, {2, 1, 0.19}, {2,
   2, 13.52}, {2, 3, 0.51}, {4, 1, 0.39}, {4, 2, 10.99}, {4, 3, 
  1.22}, {6, 1, 1.11}, {6, 2, 8.23}, {6, 3, 2.96}, {8, 1, 1.96}, {8, 
  2, 0.37}, {8, 3, 6.08}, {10, 1, 2.09}, {10, 2, 0.31}, {10, 3, 
  5.82}, {12, 1, 2.41}, {12, 2, 0.28}, {12, 3, 5.69}, {14, 1, 
  2.41}, {14, 2, 0.26}, {14, 3, 5.62}} *)

Redo the function in such a way that it will actually evaluate. One can and should test that it evaluates to a number, when fed appropriate input.

Clear[LoteGModel];
LoteGModel[Ke_?NumberQ, YxsG_?NumberQ, ms_?NumberQ, 
  YpsG_?NumberQ, \[Alpha]a_?NumberQ, \[Beta]b_?NumberQ] := 
 LoteGModel[Ke, YxsG, ms, YpsG, \[Alpha]a, \[Beta]b] =
  NDSolveValue[{{Y'[t] == (\[Mu]ma S[t]/(Ke + S[t]))*Log[Kx/Y[t]]*
       Y[t], Y[0] == 0.04},
    {S'[t] == -(((\[Mu]ma S[t]/(Ke + S[t]))*Log[Kx/Y[t]]*Y[t])/
        YxsG) - ms Y[
         t] - (\[Alpha]a*((\[Mu]ma S[t]/(Ke + S[t]))*Log[Kx/Y[t]]*
           Y[t]) + \[Beta]b*Y[t])/YpsG, 
     S[0] == 15 }, {Et'[
       t] == \[Alpha]a*( (\[Mu]ma S[t]/(Ke + S[t]))*Log[Kx/Y[t]]*
          Y[t]) + \[Beta]b*Y[t], Et[0] == 0.3}},
   {Y, S, Et}, {t, 15}]

Clear[func]
func[Ke_?NumberQ, YxsG_?NumberQ, ms_?NumberQ, 
  YpsG_?NumberQ, \[Alpha]a_?NumberQ, \[Beta]b_?NumberQ, t_?NumberQ, 
  n_?NumberQ] := 
 LoteGModel[Ke, YxsG, ms, YpsG, \[Alpha]a, \[Beta]b][[n]][t]

Now run the fit. Warnings indicate there may be trouble.

fitGomp = 
 FindFit[clists, 
  func[Ke, YxsG, ms, YpsG, \[Alpha]a, \[Beta]b, t, 
   n], {{Ke, 0.025}, {YxsG, 0.44}, {ms, 0.036}, {YpsG, 
    0.48}, {\[Alpha]a, 0.9}, {\[Beta]b, 0.08}}, {t, n}]

During evaluation of In[175]:= NDSolveValue::ndsz: At t == 7.1147359866114295`, step size is effectively zero; singularity or stiff system suspected. >>

During evaluation of In[175]:= InterpolatingFunction::dmval: Input value {8.} lies outside the range of data in the interpolating function. Extrapolation will be used. >>

During evaluation of In[175]:= InterpolatingFunction::dmval: Input value {8.} lies outside the range of data in the interpolating function. Extrapolation will be used. >>

During evaluation of In[175]:= InterpolatingFunction::dmval: Input value {8.} lies outside the range of data in the interpolating function. Extrapolation will be used. >>

During evaluation of In[175]:= General::stop: Further output of InterpolatingFunction::dmval will be suppressed during this calculation. >>

During evaluation of In[175]:= NDSolveValue::ndsz: At t == 9.45455859284449`, step size is effectively zero; singularity or stiff system suspected. >>

During evaluation of In[175]:= NDSolveValue::ndsz: At t == 9.558957684515919`, step size is effectively zero; singularity or stiff system suspected. >>

During evaluation of In[175]:= General::stop: Further output of NDSolveValue::ndsz will be suppressed during this calculation. >>

During evaluation of In[175]:= FindFit::sszero: The step size in the search has become less than the tolerance prescribed by the PrecisionGoal option, but the gradient is larger than the tolerance specified by the AccuracyGoal option. There is a possibility that the method has stalled at a point that is not a local minimum. >>

(* Out[175]= {Ke -> 8.046436990185919*10^-7, YxsG -> 0.4891746182309673, 
 ms -> 0.1809314600576526, 
 YpsG -> 0.5974714908245637, \[Alpha]a -> 
  2.426821441437095, \[Beta]b -> 0.08677711620145376} *)

One can assess that this is at least a plausible result by comparing explicit evaluations at the fitted parameter values with clists.

func2 = func[Ke, YxsG, ms, YpsG, \[Alpha]a, \[Beta]b, t, n] /. fitGomp;

Flatten[Table[{t, n, func2}, {t, 0, 14, 2}, {n, 1, 3}], 1]

(* Out[186]= {{0, 1, 0.04}, {0, 2, 15.}, {0, 3, 0.3}, {2, 1, 
  0.2770802783180547}, {2, 2, 13.46442074844808}, {2, 3, 
  0.898749873186528}, {4, 1, 0.768809615766145}, {4, 2, 
  10.13049265182656}, {4, 3, 2.180254427838512}, {6, 1, 
  1.316805932379483}, {6, 2, 6.101006017102625}, {6, 3, 
  3.691954237495139}, {8, 1, 1.748846456442421}, {8, 2, 
  2.454657792902989}, {8, 3, 5.008687899807296}, {10, 1, 
  1.948546275680985}, {10, 2, -3.578836840383028*10^-7}, {10, 3, 
  5.82196503324831}, {12, 1, 1.751107445453094}, {12, 
  2, -2.785847756390413*10^-7}, {12, 3, 5.663556487512911}, {14, 1, 
  1.573674295875138}, {14, 2, -2.280475064379094*10^-7}, {14, 3, 
  5.52119884657305}} *)

Here is a measure of relative errors, that is, discrepancies between data and fitted model.

Map[{Most[#[[1]]], 
   Abs[Last[#[[1]]] - Last[#[[2]]]]/
    Sqrt[Last[#[[1]]]^2 + Last[#[[2]]]^2]} &, 
 Transpose[{clists, evals}]]

(* Out[189]= {{{0, 1}, 0.}, {{0, 2}, 0.0473568702539026}, {{0, 3}, 
  0.}, {{2, 1}, 0.2591933469342784}, {{2, 2}, 
  0.002912820154303103}, {{2, 3}, 0.3761965972571338}, {{4, 1}, 
  0.439417580307384}, {{4, 2}, 0.05750441950496325}, {{4, 3}, 
  0.3843506227871235}, {{6, 1}, 0.120080287441452}, {{6, 2}, 
  0.2078129017888882}, {{6, 3}, 0.1546807481166434}, {{8, 1}, 
  0.0803844078497986}, {{8, 2}, 0.8397795549170736}, {{8, 3}, 
  0.1359982682960941}, {{10, 1}, 0.04950374839530081}, {{10, 2}, 
  1.00000115446283}, {{10, 3}, 0.0002387033995335336}, {{12, 1}, 
  0.2211786530788786}, {{12, 2}, 1.000000994945133}, {{12, 3}, 
  0.003293828907684501}, {{14, 1}, 0.2905634359494454}, {{14, 2}, 
  1.000000877105409}, {{14, 3}, 0.01254087995090153}} *)
POSTED BY: Updating Name

Hi, I really apreciate your answer and your help. i find out that it works!

But, let me to give you a better reply for this. Tomorrow, I will send you another message.

Best Regards. Ricardo Sánchez

POSTED BY: Ricardo Sánchez

Could you please let us know real name of the author behind Updating Name - so we can fix it.

POSTED BY: Moderation Team

Sorry, I do not have idea who is.

POSTED BY: Ricardo Sánchez

Hi updating name,

Sorry for delaying.

Here you can see that actually the data are a (x,y) pair. In my case the data are : (t, X) for the first differential equation, (t, S) for the second one, and (t,Eta) for the third one. So I would like to try with the answer that you give before, but changing the data with (x,y) pair. Like you said before, of course is a " three simultaneous fit" . Maybe is goind to be helpful for both, and thankyou again. Here: the wolfram demosntration.

In[23]:= Cinetica2 = ({ {0, 0.04, 16.04, 0.30}, {2, 0.19, 13.52, 0.51}, {4, 0.39, 10.99, 1.22}, {6, 1.11, 8.23, 2.96}, {8, 1.96, 0.37, 6.08}, {10, 2.09, 0.31, 5.82}, {12, 2.41, 0.28, 5.69}, {14, 2.41, 0.26, 5.62} });

In[24]:= DataX = {{0, 0.04}, {2, 0.19}, {4, 0.39}, {6, 1.11}, {8, 1.96}, {10, 2.09}, {12, 2.41}, {14, 2.41}}; DataS = {{0, 16.04}, {2, 13.52}, {4, 10.99}, {6, 8.23}, {8, 0.37}, {10, 0.31}, {12, 0.28}, {14, 0.26}}; DataEta = {{0, 0.30}, {2, 0.51}, {4, 1.22}, {6, 2.96}, {8, 6.08}, {10, 5.82}, {12, 5.69}, {14, 5.62}};

In[27]:= Manipulate[ Module[{Gompertz, Ks, Yxsg, Ypsg, Yps, Ms, Yo2, Coe}, Ks = 0.025; Yxsg = 0.44; Ypsg = 0.40; Yps = 0.4; Ms = 0.036; Yo2 = 1.14; Coe = 0.007; [Mu]v = [Mu]max Su[t]/(Ks + Su[t])*(Cox[t]/(ko + Cox[t])); Rb1 = ([Mu]vLog[Kx/Yb[t]]Yb[t]); REt1 = [Alpha]Rb1 + [Beta]Yb[t]; Gompertz = NDSolve[{D[Yb[t], {t, 1}] == ([Mu]vLog[Kx/Yb[t]]Yb[t]), D[Su[t], {t, 1}] == -(Rb1/Yxsg) - Ms Yb[t] - REt1/Ypsg, D[Eta[t], {t, 1}] == REt1, D[Cox[t], {t, 1}] == KlA(Coe - Cox[t]) - ([Mu]v/Yo2)Yb[t], Yb[0] == 0.1, Su[0] == 15, Eta[0] == 0.36, Cox[0] == 0.007}, {Yb[t], Su[t], Eta[t], Cox[t]}, {t, 0, 15}]; Multicolumn[{ Show[Plot[Evaluate[Yb[t] /. Gompertz], {t, 0, 15}, PlotRange -> All, Frame -> True, PlotStyle -> {Blue, Thick}, FrameLabel -> {"time (h)", "Concentration(g/L)"}, PlotLabel -> "Biomass Production", PlotLegends -> {"[Biomass]"}, PlotTheme -> "Scientific", GridLines -> Automatic, ImageSize -> Small], ListPlot[DataX, Filling -> Axis, PlotMarkers -> Style["[FilledCircle]", 12, Blue], PlotStyle -> Red, PlotLegends -> {"Experimental data"}]], Show[Plot[Evaluate[Su[t] /. Gompertz], {t, 0, 15}, PlotRange -> All, Frame -> True, PlotStyle -> {Orange, Thick}, FrameLabel -> {"time (h)", "Concentration(g/L)"}, PlotLabel -> "Substrate Uptake", PlotLegends -> {"[Ethanol]"}, PlotTheme -> "Scientific", GridLines -> Automatic, ImageSize -> Small], ListPlot[DataS, Filling -> Axis, PlotMarkers -> Style["[FilledCircle]", 12, Orange], PlotStyle -> Red, PlotLegends -> {"Experimental data"}]], Show[Plot[Evaluate[Eta[t] /. Gompertz], {t, 0, 15}, PlotRange -> All, Frame -> True, PlotStyle -> {Gray, Thick}, FrameLabel -> {"time (h)", "Concentration(g/L)"}, PlotLabel -> "Ethanol Formation", PlotLegends -> {"[Ethanol]"}, PlotTheme -> "Scientific", GridLines -> Automatic, ImageSize -> Small], ListPlot[DataEta, Filling -> Axis, PlotMarkers -> Style["[FilledCircle]", 12, Black], PlotStyle -> Black, PlotLegends -> {"Experimental data"}]]}, 1]], Item[Style["Metaboites Formation", Bold, Blue, Italic, 12], Alignment -> Center], Item[Style["Ethanol", Bold], Alignment -> Center], {{[Alpha], 0.25, "alfa"}, 0, 1, 0.01, Appearance -> "Open"}, {{[Beta], 0.01, "beta"}, 0, 1, 0.01, Appearance -> "Open"}, Delimiter, Item[Style["oxygen transfer coefficient", Bold], Alignment -> Center], {{KlA, 450, "Kla"}, 0, 500, 10, Appearance -> "Labeled"}, Delimiter, Item[Style["Growth rate", Bold], Alignment -> Center], {{[Mu]max, 0.9, "[Mu]max"}, 0.05, 1, 0.05, Appearance -> "Labeled"}, Delimiter, Item[Style["Constant oxygen affinity", Bold], Alignment -> Center], {{ko, 0.021, "kc"}, 0, 0.03, 0.001, Appearance -> "Labeled"}, Delimiter, Item[Style["Gompertz Model", Bold, Blue, Italic, 12], Alignment -> Left], Item[Style["Maximum biomass concentration", Bold], Alignment -> Center], {{Kx, 3, "K"}, 1, 3.5, 0.5, Appearance -> "Labeled"}, TrackedSymbols :> {[Alpha], [Beta], [Mu]max, ko, Kx, KlA}, SynchronousUpdating -> False, ControlPlacement -> Left]

Out[27]= Manipulate[Module[{Gompertz$, Ks$, Yxsg$, Ypsg$, Yps$, Ms$, \ Yo2$, Coe$}, Ks$ = 0.025; Yxsg$ = 0.44; Ypsg$ = 0.4; Yps$ = 0.4; Ms$ \ = 0.036; Yo2$ = 1.14; Coe$ = 0.007; [Mu]v = [Mu]max(Su[t]/(Ks$ + Su[t]))(Cox[t]/(ko \ + Cox[t])); Rb1 = [Mu]vLog[Kx/Yb[t]]Yb[t]; REt1 = [Alpha]*Rb1 + \ [Beta]*Yb[t]; Gompertz$ = NDSolve[{D[Yb[t], {t, 1}] == \ [Mu]vLog[Kx/Yb[t]]Yb[t], D[Su[t], {t, 1}] == -(Rb1/Yxsg$) - \ Ms$Yb[t] - REt1/Ypsg$, D[Eta[t], {t, 1}] == REt1, D[Cox[t], {t, 1}] == KlA(Coe$ - \ Cox[t]) - ([Mu]v/Yo2$)*Yb[t], Yb[0] == 0.1, Su[0] == 15, Eta[0] == \ 0.36, Cox[0] == 0.007}, {Yb[t], Su[t], Eta[t], Cox[t]}, {t, 0, 15}]; Multicolumn[{Show[Plot[Evaluate[Yb[t] /. Gompertz$], {t, 0, 15}, \ PlotRange -> All, Frame -> True, PlotStyle -> {Blue, Thick}, FrameLabel -> {"time (h)", "Concentration(g/L)"}, PlotLabel \ -> "Biomass Production", PlotLegends -> {"[Biomass]"}, PlotTheme -> \ "Scientific", GridLines -> Automatic, ImageSize -> Small], ListPlot[DataX, \ Filling -> Axis, PlotMarkers -> Style["[FilledCircle]", 12, Blue], \ PlotStyle -> Red, PlotLegends -> {"Experimental data"}]], \ Show[Plot[Evaluate[Su[t] /. Gompertz$], {t, 0, 15}, PlotRange -> All, \ Frame -> True, PlotStyle -> {Orange, Thick}, FrameLabel -> {"time (h)", \ "Concentration(g/L)"}, PlotLabel -> "Substrate Uptake", PlotLegends \ -> {"[Ethanol]"}, PlotTheme -> "Scientific", GridLines -> Automatic, ImageSize \ -> Small], ListPlot[DataS, Filling -> Axis, PlotMarkers -> Style["\ [FilledCircle]", 12, Orange], PlotStyle -> Red, PlotLegends -> {"Experimental data"}]], \ Show[Plot[Evaluate[Eta[t] /. Gompertz$], {t, 0, 15}, PlotRange -> \ All, Frame -> True, PlotStyle -> {Gray, Thick}, FrameLabel -> {"time (h)", \ "Concentration(g/L)"}, PlotLabel -> "Ethanol Formation", PlotLegends \ -> {"[Ethanol]"}, PlotTheme -> "Scientific", GridLines -> Automatic, ImageSize \ -> Small], ListPlot[DataEta, Filling -> Axis, PlotMarkers -> Style["[FilledCircle]", 12, Black], PlotStyle \ -> Black, PlotLegends -> {"Experimental data"}]]}, 1]], Item[Style["Metaboites Formation", Bold, RGBColor[0, 0, 1], Italic, \ 12], Alignment -> Center], Item[Style["Ethanol", Bold], Alignment -> \ Center], {{[Alpha], 0.25, "alfa"}, 0, 1, 0.01, Appearance -> "Open"}, {{\ [Beta], 0.01, "beta"}, 0, 1, 0.01, Appearance -> "Open"}, Delimiter, Item[Style["oxygen transfer coefficient", Bold], Alignment -> \ Center], {{KlA, 450, "Kla"}, 0, 500, 10, Appearance -> "Labeled"}, \ Delimiter, Item[Style["Growth rate", Bold], Alignment -> Center], {{[Mu]max, \ 0.9, "[Mu]max"}, 0.05, 1, 0.05, Appearance -> "Labeled"}, Delimiter, Item[Style["Constant oxygen affinity", Bold], Alignment -> Center], \ {{ko, 0.021, "kc"}, 0, 0.03, 0.001, Appearance -> "Labeled"}, \ Delimiter, Item[Style["Gompertz Model", Bold, RGBColor[0, 0, 1], Italic, 12], \ Alignment -> Left], Item[Style["Maximum biomass concentration", Bold], Alignment -> Center], {{Kx, 3, "K"}, 1, 3.5, 0.5, Appearance -> \ "Labeled"}, TrackedSymbols :> {[Alpha], [Beta], [Mu]max, ko, Kx, \ KlA}, SynchronousUpdating -> False, ControlPlacement -> Left]

POSTED BY: Ricardo Sánchez
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract