NDSolve behaves worse than ParametricNDSolve in FindFit

Question

I am trying to fit some model parameters to some data.

Ts = 0.001; duty = 0.2;
Ubat = 14; Ud = 0.6; tEnd = 0.005; Td = 8 10^-7;
pars = {R1 -> 11, L1 -> 0.01};

This is the RL cuircuit model driven by a PWM voltage.

u[t_, d_, T_, Ub_, Ud_] := If[Mod[t, T] < d T, Ub, -Ud]
eqn[R1_, L1_] = {i1'[t] == (u[t] - R1 i1[t])/L1};
ic[i01_] = {i1[0] == i01};

Make some noisy data.

times = N[Range[0, tEnd, Td 20]];
noisyData = 
Transpose[{times,sol[times] + RandomVariate[NormalDistribution[0, 0.01], Length[times]]}];
Show[Plot[sol[t], {t, 0, tEnd}], lp = ListPlot[noisyData, PlotRange -> All,  PlotStyle -> Red, PlotLegends -> {"data"}]]

Than I compare the two different approches:

With normal NDSolve and "memoization":

modelRLNDS[R1_?NumberQ,L1_?NumberQ,i01_?NumberQ]:=
modelRLNDS[R1,L1,i01]=(i1/.NDSolve[Flatten[{eqn[R1,L1]/.u[t]->u[t,duty,Ts,Ubat,Ud],ic[i01]}],i1,{t,tEnd},AccuracyGoal >10][[1]])
AbsoluteTiming[sol1=modelRLNDS[R1,L1,0.14]/.pars]

with ParametricNDSolve:

modelRLParaNDS=ParametricNDSolveValue[Flatten[{eqn[R1,L1]/.u[t]->u[t,duty,Ts,Ubat,Ud],ic[i01]}],i1,{t,tEnd},{R1,L1,i01},AccuracyGoal->10];
AbsoluteTiming[sol2=modelRLParaNDS[R1,L1,0.14]/.pars]

Both yield the same result in the same amount of time. But the models behave different in FindFit:

AbsoluteTiming[
 Block[{count = 0}, {fit = 
    FindFit[noisyData, 
     modelRLNDS[R1, L1, i01][x], {{R1, 1}, {L1, 0.01}, {i01, 0}}, x, 
     EvaluationMonitor :> count++, AccuracyGoal -> 5], count}]]

AbsoluteTiming[
 Block[{count = 0}, {fit2 = 
    FindFit[noisyData, 
     modelRLParaNDS[R1, L1, i01][x], {{R1, 1}, {L1, 0.01}, {i01, 0}}, 
     x, EvaluationMonitor :> count++, AccuracyGoal -> 5], count}]]

The solution with the parametric model is not only correct but it is also faster (only 6 function evaluations)!

Show[Plot[{modelRLNDS[R1, L1, i01][x] /. fit, 
   modelRLParaNDS[R1, L1, i01][x] /. fit2, 
   modelRLParaNDS[1, 0.1, 0][x]}, {x, 0, tEnd}, 
  PlotRange -> {{0, tEnd}, {0, 0.4}}, 
  PlotLegends -> {"fit1", "fit2", "starting values"}], lp]

Did I do something wrong with the normal NDSolve?

Answer 1

OK, the Update of the question turns out to be interesting.

Did I do something wrong with the normal NDSolve?

No, you've done nothing wrong, but you should notice the "normal NDSolve" approach and the ParametricNDSolve approach is different. modelRLNDS is completely a black box, while modelRLParaNDS isn't. For example, modelRLParaNDS can be used to "symbolicly" calculate the derivative of R1, L1, i01 without difficulty:

D[modelRLParaNDS[R1, L1, i01][x], R1]
% /. {R1 -> 1, L1 -> 0.01, i01 -> 0}
Plot[%, {x, 0, tEnd}]

while D[modelRLNDS[R1, L1, i01][x], R1] cannot:

D[modelRLNDS[R1, L1, i01][x], R1]
% /. {R1 -> 1, L1 -> 0.01, i01 -> 0}
Plot[%, {x, 0, tEnd}]

It's not clearly mentioned, but the document of Gradient suggests that, when calculating gradient, FindFit performs symbolic differentiation if the model can be symbolically differentiated, and use finite difference if the model is a blackbox, so it's not surprising that the "normal NDSolve" model and the ParametricNDSolve model gives different fitting result. Still, fit can be improved, by setting "DifferenceOrder" -> 2:

AbsoluteTiming[
 Block[{count = 0}, {fit = 
    FindFit[noisyData, modelRLNDS[R1, L1, i01][x], {{R1, 1}, {L1, 0.01}, {i01, 0}}, x, 
     EvaluationMonitor :> count++, 
     Gradient -> {"FiniteDifference", "DifferenceOrder" -> 2}], count}]]
(* {0.222677, {{R1 -> 11.0405, L1 -> 0.0100116, i01 -> 0.141926}, 99}} *)