lsqnonlin curve fitting to the output of ODE45

[James]

Dear Matlab Community,

I am trying to fit the last vector of a coupled ODE (5 ode) to a set of simulated data by optimizing the set of parameters in the coupled ODE. I could be getting lsqnonlin syntax wrongly, so please advise.

This is the file (runode1.m) that runs the simulated data to be fitted:

function runode1
tspan = 0:0.1:100;
init_cond1 = [.1];

[t, x] = ode45(@ode_fun1, tspan, init_cond1);

save('ode1data');

figure; plot(t, x(:, 1),'b');;

end

function dxdt = ode_fun1(t, x)
dxdt = zeros(1, 1);
n=2;
k=.18; %
beta=0.36;

dx1dt = beta * (x(1)^n / (k^n + x(1)^n))- 0.18*x(1);
dxdt(1) = dx1dt;
end


The data to be fitted can be loaded from ode1data as (t,x).
Next this is the coupled ODE system (fit_simp.m). I would like to fit the last vector, in this case sol(:,5) to the simulated data above (t,x) optimizing the param.


function diff = fit_simp(param)

init_cond2 = [2 2 0.0015 0.0005 0.1];
tspan = 0:0.1:100;

global kf1 kr2 kfd1 krd1 kf2

kf1=param(1);
kr2=param(2);
kfd1=param(3);
krd1=param(4);
kf2=param(5);

sol=ode45(@myeqn, tspan, init_cond2);
diff = sol(:,5);

function dydt=myeqn(y)
dydt(1) = -kf1*y(1)*y(2) + kr1*y(3);
dydt(2) = -kf1*y(1)*y(2) + kr1*y(3);
dydt(3) = kf1*y(1)*y(2) - kr1*y(3) - 2*kfd1*y(4) + 2*krd1*y(4);
dydt(4) = 2*kfd1*y(3)^2 - krd1*y(4) - kf2*y(4);
dydt(5) = kf2*y(4);
end
end


To call lsgnonlin i use the following command:
t=tspan;
param0=[1 1 1 1 1]';
param=lsqnonlin(@fit_simp,param0,[],[],[],t,x);

The error seems to point that my number of input arguments is not just sol(:,5). I would appreciate if someone could point out my mistakes?
Thanks in advance.

Yours sincerely,
James

[Torsten]

"James" wrote in message <kndt51$rin$1...@newscl01ah.mathworks.com>...

If t and x are vectors of the same length, call lsqnonlin as

param=lsqnonlin(@(param)fit_simp(param,t,x),param0);

and organize fit_simp as

function diff = fit_simp(param,t,x)

global kf1 kr2 kfd1 krd1 kf2

kf1=param(1);
kr2=param(2);
kfd1=param(3);
krd1=param(4);
kf2=param(5);

init_cond2 = [2 2 0.0015 0.0005 0.1];

sol=ode45(@myeqn, t, init_cond2);
diff = sol(:,5)-x;

function dydt=myeqn(y)
dydt = zeros(5,1);

dydt(1) = -kf1*y(1)*y(2) + kr1*y(3);
dydt(2) = -kf1*y(1)*y(2) + kr1*y(3);
dydt(3) = kf1*y(1)*y(2) - kr1*y(3) - 2*kfd1*y(4) + 2*krd1*y(4);
dydt(4) = 2*kfd1*y(3)^2 - krd1*y(4) - kf2*y(4);
dydt(5) = kf2*y(4);
end
end

Best wishes
Torsten.

[James]

"Torsten" wrote in message <knf5kl$iia$1...@newscl01ah.mathworks.com>...

Thanks Torsten for the reply. Vectors t and x are of the same length. I tried your recommendation, with a little tweak, it almost work! However an error event took place due to tolerance issue with the integration process due to tolerance.
I could not understand why the coupled ODE integrates smoothly when ran by itself but faces tolerance issue when ran with lsqnonlin.
Appreciate your inputs on this.

Warning: RelTol has been increased to 2.22045e-014.
> In funfun\private\odearguments at 128
In ode45 at 114
In fit_simp at 15
In @(param)fit_simp(param,t,x)
In lsqnonlin at 199
.
.
.

Here's the fit_simp.m that I have tweaked from your recommendation:

function diff = fit_simp(param,t,x)

global kf1 kr1 kfd1 krd1 kf2

kf1=param(1);
kr1=param(2);
kfd1=param(3);
krd1=param(4);
kf2=param(5);

tspan=0:0.1:100;

init = [2 2 0.0015 0.0005 0.1];
options = odeset('RelTol',1e-200,'AbsTol',[1e-200 1e-200 1e-200 1e-200 1e-200]);
[t,y]=ode45(@myeqn, tspan, init,options);
save('test');
diff = y(:,5);

function dydt=myeqn(t,y)
dydt = zeros(5, 1);

[Steven_Lord]

"James " <hsuk...@gmail.com> wrote in message
news:kng7ic$cm$1...@newscl01ah.mathworks.com...

> "Torsten" wrote in message <knf5kl$iia$1...@newscl01ah.mathworks.com>...
>> "James" wrote in message <kndt51$rin$1...@newscl01ah.mathworks.com>...

*snip*

> Thanks Torsten for the reply. Vectors t and x are of the same length. I
> tried your recommendation, with a little tweak, it almost work! However
> an error event took place due to tolerance issue with the integration
> process due to tolerance. I could not understand why the coupled ODE
> integrates smoothly when ran by itself but faces tolerance issue when ran
> with lsqnonlin.
> Appreciate your inputs on this.
>
> Warning: RelTol has been increased to 2.22045e-014.
>> In funfun\private\odearguments at 128
> In ode45 at 114
> In fit_simp at 15
> In @(param)fit_simp(param,t,x)
> In lsqnonlin at 199 .
>

> Here's the fit_simp.m that I have tweaked from your recommendation:
>
> function diff = fit_simp(param,t,x)
>
> global kf1 kr1 kfd1 krd1 kf2
>
> kf1=param(1);
> kr1=param(2);
> kfd1=param(3);
> krd1=param(4);
> kf2=param(5);
>
> tspan=0:0.1:100;
>
> init = [2 2 0.0015 0.0005 0.1];
> options = odeset('RelTol',1e-200,'AbsTol',[1e-200 1e-200 1e-200 1e-200
> 1e-200]);

Don't use tolerances this small. By the definitions in the documentation:

http://www.mathworks.com/help/matlab/ref/odeset.html

your RelTol value is asking for 1e-198% accuracy, or

0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001%
accuracy.

Leaving RelTol at its default value of 1e-3, or at most reducing it down to
1e-6 or so, should be sufficient for most problems. The warning you received
was ODE45 telling you it was ignoring your request for it to use a RelTol of
1e-200 and using 100*eps instead. Similarly I would not reduce AbsTol below
1e-10, 1e-12 without a strong reason and a solution of very small magnitude
(and even then, I'd probably rescale the problem instead, like using 100
meters instead of 0.1 kilometers as the quantities in my problem.)

*snip rest of code*

--
Steve Lord
sl...@mathworks.com
To contact Technical Support use the Contact Us link on
http://www.mathworks.com

[James]

"Steven_Lord" <sl...@mathworks.com> wrote in message <kngbv2$et8$1...@newscl01ah.mathworks.com>...

Thanks Steve, I have tried from 1e-3 till 1e-20, the latter I have suspected of "overkill" tolerance set. Having said that, thank for the advice. I will stick to a reasonable tolerance. I think there is some underlying issue or mistake that I have made when I call ode45 from lsqnonlin. I m still trying to fix the problem as we speak. Thank you all for reading and posting useful hint. I m much closer to solving this problem!

James Ooi

[Bjorn Gustavsson]

"James" wrote in message <kngk1s$d41$1...@newscl01ah.mathworks.com>...

Might one problem be (become) that lsqnonlin runs into some corner of parameter space where your ODE becomes stiff (or worse)? If so you might have to constrain the search-space - or use one of the other ode-integration functions.

[James]

"Bjorn Gustavsson" <bj...@irf.se> wrote in message <kngshj$8ta$1...@newscl01ah.mathworks.com>...

I did try ode23s thinking along the same line as you but it did really help. What actually works (partially) is that I made "param" a global variable. Then lsqnonlin was able to call fit_simp.m.
However lsqnonlin did not do what it is intend to, ie. to optimize/ or fit to the data t,x supplied while varying the param. I see 6 iterations that lsqnonlin call fit_simp.m and every iteration did nothing but just print the supplied parameters. At the 6th iteration lsqnonlin basically give up, I am thinking of changing the lsqnonlin stopping criteria. Any suggestions?

"Initial point is a local minimum.

Optimization completed because the size of the gradient at the initial point
is less than the default value of the function tolerance."


here's my latest fit_simp.m
function diff = fit_simp(param,t,x)

global kf1 kr1 kfd1 krd1 kf2 param

kf1=param(1);
kr1=param(2);
kfd1=param(3);
krd1=param(4);
kf2=param(5);

% param(1)=kf1;
% param(2)=kr1;
% param(3)=kfd1;
% param(4)=krd1;
% param(5)=kf2;

tspan=0:0.1:100;

init = [2 2 0.0015 0.0005 0.1];

% options = odeset('RelTol',1e-3,'AbsTol',[1e-6 1e-6 1e-6 1e-6 1e-6]);
[t,y]=ode45(@myeqn, tspan, init);

save('test');
diff = y(:,5);

figure;plot(t,y(:,5),'g',t,x,'r');

[Marc]

"James" wrote in message <kngvom$h6p$1...@newscl01ah.mathworks.com>...

I have not had good results with nlinfit, lsqnonlin when fitting systems of odes. I would set this up like an optimization problem, finding the least squares yourself. This way you can also play with the weights. Then try fminsearch. It's slower but the nelder mead simplex does not require taking differentials, so I have found it to converge with much better results.

Also, I would try ode15s. If I have some time tomorrow i will give your system a shot.

I am not sure why you have the params as global. Since you are nesting the functions, this should not be needed.

[Torsten]

"James" wrote in message <kng7ic$cm$1...@newscl01ah.mathworks.com>...

Didn't you notice in my partial code that you have to return
diff=y(:,5)-x
to lsqnonlin ?

Best wishes
Torsten.

[James]

"Torsten" wrote in message <knhpdu$na5$1...@newscl01ah.mathworks.com>...

Hi Torsten,
Yes I did, and i did try it. I know it makes sense as lsqnonlin is minimizing the difference between the ouput of ODE45 and X. I still get 6 itereation without making the difference any smaller. Sorry about my previous fit_simp.m being incomplete.

Here's what I have tried and still work in progress:

%fit_simp.m
function diff = fit_simp(param,t,x)

global kf1 kr1 kfd1 krd1 kf2 param

kf1=param(1);
kr1=param(2);
kfd1=param(3);
krd1=param(4);
kf2=param(5);

tspan=0:0.1:100;

init = [2 2 0.0015 0.0005 0.1];

[t,y]=ode45(@myeqn, tspan, init);
save('test');
diff = y(:,5)-x;

figure;plot(t,y(:,5),'g',t,x,'r');

[James]

"James" wrote in message <knit3q$a1d$1...@newscl01ah.mathworks.com>...

Dear all,

Thank you all for reading and providing me suggestions/corrections.
I would like to report that I have successfully utilize the lsqnonlin for my problem.
Here I want to recap so that the solution may be used for others in the future.

First you need your data to be fitted to. [t,x] where t and x are 1001x1 matrix respectively. Then you need to set your initial parameter guesstimate =>

param0 = [0.150000000000000
0.0100000000000000
1
0.0100000000000000
2]

Then call the lsqnonlin using:
newparam=lsqnonlin(@(param)fit_simp(param,t,x),param0,[0 0 0 0 0],[]);

Here is the code that works:

%fit_simp.m
function diff = fit_simp(param,t,x)

kf1=param(1);
kr1=param(2);
kfd1=param(3);
krd1=param(4);
kf2=param(5);

tspan=0:0.1:100;

init = [2 2 0.0015 0.0005 .1];
[t,y]=ode23s(@myeqn, tspan, init);

save('test');
diff = y(:,5)-x;

function dydt=myeqn(~,y)

dydt = zeros(5, 1);
dydt(1) = -kf1*y(1)*y(2) + kr1*y(3);
dydt(2) = -kf1*y(1)*y(2) + kr1*y(3);
dydt(3) = kf1*y(1)*y(2) - kr1*y(3) - 2*kfd1*y(4) + 2*krd1*y(4);
dydt(4) = 2*kfd1*y(3)^2 - krd1*y(4) - kf2*y(4);
dydt(5) = kf2*y(4);
end
end

I hope this solution can help others. Thank you!