2
$\begingroup$

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

enter image description here

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]

enter image description here

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]

enter image description here

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

enter image description here

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

enter image description here

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]

enter image description here

Did I do something wrong with the normal NDSolve?

$\endgroup$
  • 2
    $\begingroup$ The problem with the fit is, that your model is evaluated again and again for every point [x]. To speed up the fit you could use "memoization": model[R1_?NumericQ, L1_?NumericQ] := model[R1, L1] = i1 /. NDSolve[ Flatten[{eqn[R1, L1] /. u[t] -> u[t, duty, Ts, Ubat, Ud], ic}], i1, {t, tEnd}][[1]]. But your fit depends strongly on the starting values ... which would be another question $\endgroup$ – grbl Feb 7 '17 at 9:31
  • 1
    $\begingroup$ Guys, please don't vote to close this question, the update of the question is no longer simple, see my answer for more details. BTW @ohmsweetohm I think you can clean up your question a bit. $\endgroup$ – xzczd Feb 8 '17 at 3:45
3
$\begingroup$

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

Mathematica graphics

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

Mathematica graphics

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}} *)
$\endgroup$
  • $\begingroup$ Thank you for your answer! Would "memoization" of the ParmetricNDSolve model help to improve the performance? $\endgroup$ – OhmSweetOhm Feb 8 '17 at 18:17
  • $\begingroup$ @OhmSweetOhm I don't think so, because as far as I can tell ParametricFunction already owns memorization ability. (Just try e.g. modelRLParaNDS[1, 0.02, 0]//AbsoluteTiming twice and check the timing. ) $\endgroup$ – xzczd Feb 9 '17 at 2:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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