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