Trouble fitting ode45

Have anyone find that fitting using optimization toolbox (lsqcurvefit, curvefit, etc) doesn't fit well with ode45 or the ode solvers in general? I have been playing with this for months, and I have came to a conclusion that the fitting from optimization toolbox cannot fit with ode solutions from ode library that matlab provided. I found that (maybe), the fitting can't fit is due to the variable step size (h) in the ode solver.

[Star Strider]

To correct for the variable step size, use the independent variable of your data as the ‘tspan’ argument to the ODE solver you are using. Also, experiment with different solvers. For example with a ‘stiff’ problem, use ode15s.

[Arbol]

I correct with tspan for the right increasing time interval, this still doesn't help. I also tried other ODE solvers too, probably all of them, and no luck :(. The variable time step is nice though, but the fitting library couldn't fit nicely.

Answer 1

[John D'Errico]

Fitting the solution of an ODE solver can be difficult, since the ODE solver is an adaptive method, with a tolerance on the result. That means there will be tiny variations in the results, to within that tolerance.

But the fitting tool must now differentiate the result. It does so by making tiny changes, and looking at the result. After all, a derivative is simply a limit of deltay/deltax, for small deltax. There is a problem in this. Remember, we said that we expect tiny, uncontrolled variations in y, within the convergence tolerance. But divide them by a tiny value of deltax, and they may now be highly significant. In fact, those tiny variations in y may actually dominate the Jacobian estimate.

This can be an issue with optimizations on top of other methods, such as numerical integration, rootfinders, etc. Essentially, anytime you layer different adaptive methods, I would want to take great care this is not an issue.

There may be other reasons for a fitting problem. Do we know that you have posed a reasonable model for this fit?

One big reason for curve fitting problems is often that the user has posed a model that has no chance of fitting the data they give to the fitting tool. A second reason for failure is a terribly poor choice of starting values for the parameters.

But we are given no clue as to what the problem is, only that you have had problems in fitting.

[Arbol]

Hi John, If you can read my comment that I left for Walter that would be great! I have also tried multisearch to look into different points, but this didn't do much either. I'm aware that I have to have good model that can fit my data, this model, I have gone over multiple times this last two months and I'm pretty sure that it is the correct model to fit with.

Here is my fit if you want to take a look at, it is a two compartmental model:

tu is 225x3; tdatas=tu(:,1),tissdata=tu(:,3);

param0=[0.0001 0.0001  0.000011    0.00001];
fun = @(p,tdatas) objfunction(p,tdatas,tu);
z=lsqcurvefit(fun,param0,tdatas(:,1),...
    tissdata,[0 0 0 0],[1 1 1 1],options);

%%Functions
%%%Model
function Tissue= objfunction(p,texp,tu)
init = [0 0];
[~,y]=ode15s(@(t,x) myeqn(t,x,p,tu), texp, init);
P=y(:,1).*p(2);
I=y(:,2).*p(3);
Tissue = I+P;
end

function dx=myeqn(t,x,p,tu)
F=p(1);fp=p(2);fis=p(3);PS=p(4);
tdata=tu(:,1); adata=tu(:,2);
output=interpn(tdata,adata,t); 
dx=zeros(2,1);
dx(1,1)= (F./fp).*(output./(1-0.45)-x(1))- (PS./fp).*(x(1)-x(2));
dx(2,1)=  (PS./fis).*(x(1)-x(2))

end

Futhermore, I have some results from my fit. Here: This does look well fitted but the result is not consistent.

Now, let's compare to my RK4 (fixed step) ODE solver (it takes a way way longer to fit) and it can give consistent result no matter what my starting point is:

If you want to try, I also attached my data in in here. Where tu=load('test.txt').

[Arbol]

Well, I take back it can't fit the ode solver, it can, but the fitting is not consistent and not as good as the fixed-step time.

Answer 2

[Walter Roberson]

My experience with your system from your earlier question https://www.mathworks.com/matlabcentral/answers/344290-is-this-the-correct-optimization-tool-lsqcurvefit-multistart was that you were asking for high precision output for a calculation that tends to get stuck in local minima with some notable differences between the "basins of attraction". So you spend a lot of time figuring out precisely the best coefficient for an area where the function might be twice as large as the best value. That is a waste of time. You would be better off surveying multiple areas with lower precision on the calculations, looking for the basins of attraction, and then only narrowing in on the best few.

This has nothing to do with the variable step-size on the solver: your goal is global optimization but you are spinning around a lot in local optimization.

If you had reason to expect that a small change in parameters could lead to a big change in function value (which is the case for some problems), then you probably should not be using interp1 inside your ode, as that interp1 amounts to a guess about how the function behaves over that interval: an expectation of a large change in resulting function value is implicitly an admission that the data being interpolated over is unlikely to be well represented by a series of straight line segments. If you had a model of the underlying function of the data you are interpolating, you would probably be better off doing a series of taylor expansions to at least a couple of degrees and using those in a piecewise fashion. Perhaps you could even justify piecewise cubic spline (which would save some trouble as that kind of interpolation is built in.)

[Arbol]

Hey Walter, thank you for your follow up. I have tried global search with multistart and patternsearch, nonlinear, or whatever matlab has to offer, the fitting is not quite good. I even tried to weight the data. You do have a good point about my interp function, this is not what I see when I try to fit with a fixed step size that I made. Furthermore, I want to change my interp function to use in my RK4 fixed step.

Whenever I try to fit the ode solver, it would not converge. However, the fixed model that I create would converge to a local minimum. But the fixed model (RK4) that I made takes a lot more time than the ode solvers that matlab provide. I'm slowly getting there.

I have some figures of the ODE45 vs the fixed step that I fitted, look in the above comment. In John's section.

Answer 3

[Arbol]

Hey everyone, I have figured out how to fit my data. You guys are right about the model. I was able to fit one data set. But other data set I couldn't fit because of different units that they are in. I found that changing the odeset also can influence in the fitting after reading the Optimization Toolbox Manual! Thank you for all of your follow up!