Minimize the sum-of-squares errors for a system of ODEs

Question

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]];

The new function ParametricNDSolve[] seems to be useful for this situation. — J. M.'s missing motivation, May 04 '15 at 19:29
@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 /. ? — E3labs, May 04 '15 at 21:57

Michael E2 · Accepted Answer · 2015-05-06T15:26:06.483

[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.

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? — E3labs, May 05 '15 at 18:20
@e3labs Which one is faster? I found my sumsquare1 about 4 times faster. What was your input exactly? — Michael E2, May 05 '15 at 18:56
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? — E3labs, May 06 '15 at 12:57
@e3labs Table[sumsquare1[kp1, kp2, theta], {kp1, 0, 25}, {kp2, 1, 10}, {theta, 0.1, 4}]; // AbsoluteTimingtakes 4.5 ± 0.2 sec with my definition and 23-24 sec with Daniel's -- I have no idea why, though. The lengthLength[rtnTime]is 9, so I can't test 25. Comparing values, there are some slight differences (mainly less than10^-12but some as high as6*10^-10`. It could be than minimizing a numerical approximation is sensitive to such slight oscillations. — Michael E2, May 06 '15 at 14:12

score 2 · Answer 2 · edited May 05 '15 at 02:36

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]}]]

Minimize the sum-of-squares errors for a system of ODEs

2 Answers2

Using DSolve

Original answer

Using `DSolve`