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

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. May 4, 2015 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 /. ? May 4, 2015 at 21:57

[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}}}
*)


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, θ},
"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? May 5, 2015 at 18:20
• @e3labs Which one is faster? I found my sumsquare1 about 4 times faster. What was your input exactly? May 5, 2015 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? May 6, 2015 at 12:57
• @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. May 6, 2015 at 14:12

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