1
$\begingroup$

I'm trying to minimizing the sum of squared errors on a 2-equation system of ODEs and empirical data. I keep getting an error and I think the root of the problem is my invocation of ReplaceAll.

Here is my code:

myode1[kp1_, kp2_, θ_] := 
   ParametricNDSolveValue[{
     p1'[t] == 1/HRT (0.2 p1in - p1[t]) - kp1*Exp[θ (temp - 20)] p1[t],
     p2'[t] == 1/HRT (0.8 p1in - p2[t]) - kp2*Exp[θ (temp - 20)]  p2[t] + 
       kp1*Exp[θ (temp - 20)] p1[t],
     p1[0] == 0.2 p1in, p2[0] == 0.8 p1in}, 
     {p1, p2}, {t, 0, 350}, {kp1, kp2, θ}];

sumsquare1[kp1_?NumericQ, kp2_?NumericQ, θ_?NumericQ] := 
  Sqrt[Sum[(effpartdata[[i]] - p2[t] /. myode1[kp1, kp2, θ] /. 
    {t -> rtnTime[[i]], temp -> tempCdata[[i]], HRT -> HRTdata[[i]], 
     p1in -> totalCODin[[i]]})^2, {i, 1, 25}]] // Quiet

fitter = NMinimize[sumsquare1[kp1, kp2, θ], {kp1, kp2, θ}]

The arrays: HRTdata, totalCODin, tempCdata, effpartdata and rtnTime are my data arrays. Added them below. I've seen a couple other data fitting questions (e.g. here & here), but mine differs in that I have multiple data trails and as of right now I've been unable to tweak the other codes to accomodate that.

Any ideas?

EDIT: Tried using ParametricNDSolveValue and calling Quiet on my sumsquare1 function to no avail. thanks @Guess who it is

partDATA1 = {{1.`, 1.`, 8.88889`, 2652.`, 156.`, 261.`, 59.`, 2391.`, 97.`, 0.959`}, {1.`, 1.`, 8.889`, 2652.`, 170.`, 261.`, 51.`, 2391.`, 119.`, 0.950`}, {0.5`, 0.166`, 12.88`, 2544.`, 663.`, 493.`, 317.`, 2051.`, 346.`, 0.831`}, {1.`, 1.`, 7.22`, 2014.`, 100.`, 134.`, 32.`, 1880.`, 68.`, 0.963`}, {1.`, 1.`, 13.166`, 1902.`, 66.`, 116.`, 8.`, 1786.`, 58.`, 0.967`}, {1.`, 1.`, 13.166`, 1902.`, 48.`, 159.`, 10.`, 1743.`, 38.`, .978`}, {0.5`, 0.166`, 11.66`, 1964.`, 349.`, 224.`, 134.`, 1740.`, 215.`, 0.8764367816091954`}, {0.5`, 0.166`, 14.722`, 1568.`, 406.`, 203.`, 147.`, 1365.`, 259.`, 0.81`}, {0.5`, 0.166`, 14.722`, 1568.`, 406.`, 203.`, 143.`, 1365.`, 263.`, 0.8073260073260073`}, {1.`, 1.`, 6.833333333333332`, 1434.`, 92.`, 163.`, 62.`, 1271.`, 30.`, 0.97639653815893`}, {1.`, 1.`, 6.833333333333332`, 1434.`, 110.`, 163.`, 58.`, 1271.`, 52.`, 0.959087332808812`}, {1.`, 1.`, 7.000000000000001`, 1073.`, 150.`, 150.`, 70.`, 923.`, 80.`, 0.9133261105092091`}, {1.`, 1.`, 7.000000000000001`, 1073.`, 116.`, 150.`, 74.`, 923.`, 42.`, 0.9544962080173348`}, {1.`, 0.166`, 15.55`, 932.`, 361.`, 202.`, 138.`, 730.`, 223.`, 0.694`}, {1.`, 1.`, 7722`, 937.`, 78.`, 273.`, 71.`, 664.`, 7.`, 0.9892`}, {1.`, 1.`, 7.22`, 937.`, 97.`, 273.`, 65.`, 664.`, 32.`, 0.951`}, {1.`, 1.`, 13.333`, 479.`, 132.`, 116.`, 48.`, 363.`, 84.`, 0.768`}, {1.`, 0.66`, 20.557`, 670.`, 150.`, 310.`, 74.`, 360.`, 76.`, 0.88`}, {0.5`, 0.166`, 12.8`, 663.`, 306.`, 317.`, 233.`, 346.`, 73.`, 0.789`}, {1.`, 0.33`, 19.33`, 500.`, 200.`, 234.`, 117.`, 266.`, 83.`,0.687`}, {0.5`, 0.66`, 14.22`, 406.`, 196.`, 143.`, 92.`, 263.`, 104.`,0.604`}, {0.5`, 0.166`, 14.22`, 406.`, 271.`, 147.`, 102.`, 259.`, 169.`,0.3474`}, {1.`, 0.66`, 20.55`, 328.`, 148.`, 124.`, 64.`, 204.`, 84.`,0.58`}, {1.`, 1.`, 13.05`, 391.`, 75.`,230.`, 105.`, 161.`, 10.`, 0.937`}, {1.`,0.33`, 19.33`, 244.`, 122.`, 234.`,120.`, 10.`, 2.`, 0.8`}};
HRTdata = partDATA1[[All, 1]];
totalCODin = partDATA1[[All, 4]];
tempCdata = partDATA1[[All, 3]];
effpartdata = partDATA1[[All, 9]];
rtnTime = partDATA1[[All, 2]];
Improve this question
$\endgroup$
2
  • $\begingroup$ The new function ParametricNDSolve[] seems to be useful for this situation. $\endgroup$
    – J. M.'s persistent exhaustion ♦
    May 4, 2015 at 19:29
  • $\begingroup$ @Guesswhoitis. I gave ParametricNDSolve[] a go with no success. the error is coming up when I use NMinimize or FindMinimum. I suspect it is because of All the parameters I am trying to /. ? $\endgroup$
    – E3labs
    May 4, 2015 at 21:57

2 Answers 2

Reset to default
4
$\begingroup$

[Edit notice: It turns out the system can be integrated exactly, which leads to faster performance.]

Using DSolve

I replaced the approximate coefficients 0.2 and 0.8 by the exact numbers 2/10 and 8/10; otherwise, things that should cancel out do not and lead to false singularities.

myode1 = Function[{kp1, kp2, θ, temp, HRT, p1in}, 
   Evaluate@
    First@DSolve[{p1'[t] == 
        1/HRT (2/10 p1in - p1[t]) - kp1*Exp[θ (temp - 20)] p1[t], 
       p2'[t] == 
        1/HRT (8/10 p1in - p2[t]) - kp2*Exp[θ (temp - 20)] p2[t] + 
         kp1*Exp[θ (temp - 20)] p1[t],
       p1[0] == 2/10 p1in, p2[0] == 8/10 p1in}, {p1, p2}, t]];

Clear[sumsquare1C];
sumsquare1[kp1_, kp2_, θ_] =
  Sqrt[Sum[(effpartdata[[i]] - p2[t] /. 
        myode1[kp1, kp2, θ, tempCdata[[i]], HRTdata[[i]], 
         totalCODin[[i]]] /. t -> rtnTime[[i]])^2, {i, 1, 
     Length[rtnTime]}]];

fitter = NMinimize[{sumsquare1C[kp1, kp2, θ],
    {0 < kp1 < 25, 1 < kp2 < 10, 0.1 < θ < 4}}, {kp1, kp2, θ}] // 
  AbsoluteTiming
(*
  {3.53719, {1289.68, {kp1 -> 0., kp2 -> 10., θ -> 0.1}}}
*)

Original answer

Use ParametricNDSolve instead of ParametricNDSolveValue. It returns replacement rules for replacing p1 and p2 by their corresponding ParametricFunction. Note the arguments to a ParametricFunction are the parameters and you do not need them to be arguments to myode1.

There are parameters HRT, temp, and p1in that are not declared as parameters for ParametricNDSolve, but they should be.

myode1 = ParametricNDSolve[
   {p1'[t] == 1/HRT (0.2 p1in - p1[t]) - kp1*Exp[θ (temp - 20)] p1[t], 
    p2'[t] == 1/HRT (0.8 p1in - p2[t]) - kp2*Exp[θ (temp - 20)] p2[t] + 
      kp1*Exp[θ (temp - 20)] p1[t],
    p1[0] == 0.2 p1in, p2[0] == 0.8 p1in},
   {p1, p2}, {t, 0, 350}, {kp1, kp2, θ, temp, HRT, p1in}];

sumsquare1[kp1_?NumericQ, kp2_?NumericQ, θ_?NumericQ] := 
 Sqrt[Sum[
   (effpartdata[[i]] -
       p2[kp1, kp2, θ, tempCdata[[i]], HRTdata[[i]], totalCODin[[i]]][t] /.
     myode1 /. t -> rtnTime[[i]])^2,
   {i, 1, 9}]]
fitter = NMinimize[
  sumsquare1[kp1, kp2, θ], {kp1, kp2, θ}, 
  Method -> {"NelderMead", 
    "InitialPoints" -> {{0, 1, -0.1}, {1, 0, -0.1}, {0.1, 1, 0}, {1, 1, -1}}}]
(*
  {72.1155, {kp1 -> -0.2426, kp2 -> 9.28932, \[Theta] -> -0.0641051}}
*)

I get some ParametricNDSolve::mxst errors that may be of concern. Perhaps the initial points are not well chosen. I'll leave that to the OP.

Improve this answer
$\endgroup$
4
  • $\begingroup$ As always-Thank you so much! Your guys' help goes farther than you know. If you don't mind, can you elaborate on your use of ParametricNDSolve instead of ParametricNDSolveValue? And is that the reason why @DanielLichtblau 's response runs So Much faster than the above? Or is it because of the algorithmic differences in the function sumsquare1? $\endgroup$
    – E3labs
    May 5, 2015 at 18:20
  • $\begingroup$ @e3labs Which one is faster? I found my sumsquare1 about 4 times faster. What was your input exactly? $\endgroup$
    – Michael E2
    May 5, 2015 at 18:56
  • $\begingroup$ Daniel's is faster. I used exactly what you had there. Except my dataset has a length of 25 (updated in OP). And I added a constraint on NMinimize: 10<kp1<25, 1<kp2<10, 0.1<theta<4. BTW I tried using {i,1,1} instead of 25 (or 9) and it took about 10 seconds to run. Could it be because it cycles through all of the data points? $\endgroup$
    – E3labs
    May 6, 2015 at 12:57
  • $\begingroup$ @e3labs Table[sumsquare1[kp1, kp2, theta], {kp1, 0, 25}, {kp2, 1, 10}, {theta, 0.1, 4}]; // AbsoluteTiming` takes 4.5 ± 0.2 sec with my definition and 23-24 sec with Daniel's -- I have no idea why, though. The length Length[rtnTime] is 9, so I can't test 25. Comparing values, there are some slight differences (mainly less than 10^-12 but some as high as 6*10^-10. It could be than minimizing a numerical approximation is sensitive to such slight oscillations. $\endgroup$
    – Michael E2
    May 6, 2015 at 14:12
2
$\begingroup$

Here is one variant that seems to evaluate to a number.

myode1[p1in_?NumericQ, HRT_?NumericQ, temp_?NumericQ] := 
  ParametricNDSolveValue[{p1'[t] == 
     1/HRT (0.2 p1in - p1[t]) - kp1*Exp[θ (temp - 20)] p1[t], 
    p2'[t] == 
     1/HRT (0.8 p1in - p2[t]) - kp2*Exp[θ (temp - 20)] p2[t] + 
      kp1*Exp[θ (temp - 20)] p1[t], p1[0] == 0.2 p1in, 
    p2[0] == 0.8 p1in}, {p1[t], p2[t]}, {t, 0, 350}, {kp1, 
    kp2, θ}];

sumsquare1[kp1_?NumericQ, kp2_?NumericQ, θ_?NumericQ] := 
 Sqrt[Total @
   Table[(effpartdata[[i]] - (With[{temp = tempCdata[[i]], 
           HRT = HRTdata[[i]], p1in = totalCODin[[i]]}, 
          myode1[p1in, HRT, temp][kp1, kp2, θ][[2]]] /. {t -> 
           rtnTime[[i]]}))^2, {i, 1, Length[rtnTime]}]]
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.