Parameter estimation using fminsearch and ode45

Dear all, I have been trying to estimate 3 parameters that exist in one ordinary differential equation by fitting it to experimental data. I used fminsearch and sum of least squares for minimisation and ode45 to solve my ode. The curve fitting is more or less ok though it takes long time but the parameters that I found at the end (a,b,c) are far from the expected. Could you please debugg my code? Here is my code with its three parts. Help is much appriciated. Thank you

%function ode
    function dy = nmodel(t,y,a,b,c)
    T=423+((10/60).*(t-a));
    K=b.*exp(-c/(0.0083144621.*T));
    dy=K.*(1-y);

%least_square  
function val=nlstsq(k)
global a;
global b;
global c;
a=k(1); b=k(2);c=k(3);
mytime=[0 1  45  90  135  180  225  270  315  360  405  450  495  540  585  630  675  720  765  810  855  900  945  990  1035  1080  1125  1170  1215  1260  1305  1350  1395  1440  1485  1530  1575  1620  1665  1710  1755  1800  1845  1890  1935  1980  2025  2070  2115  2160  2205  2250  2295  2340  2385  2430  2475  2520  2565  2610  2655  2700  2745  2790  2835  2880  2925  2970  3015  3060  3105  3150  3195  3240  3285  3330  3375  3420  3465  3510  3555  3600  3645  3690  3735  3780  3825  3870  3915  3960  4005  4050  4095  4140  4185  4230  4275  4320  4365  4410  4455  4477]';

 mydata(:,1)=[0 -0.000002  -0.000052  -0.00001  -0.00001  0.000033  0.000004  0.000153  0.000287  0.00035  0.000465  0.000722  0.001095  0.001648  0.002188  0.003105  0.004489  0.006287  0.008965  0.013339  0.019986  0.033034  0.05822  0.11007  0.209152  0.374219  0.593259  0.797734  0.907849  0.93628  0.943445  0.948166  0.951732  0.954762  0.957273  0.959441  0.961251  0.963055  0.964656  0.966058  0.96702  0.968447  0.969607  0.970655  0.971693  0.973015  0.973763  0.97464  0.975294  0.976192  0.976843  0.977653  0.978123  0.978825  0.979442  0.980173  0.980687  0.981281  0.981765  0.982131  0.98271  0.983251  0.983876  0.98423  0.984645  0.985196  0.985716  0.986188  0.986665  0.986943  0.98748  0.988026  0.988459  0.98909  0.989243  0.989685  0.990413  0.990708  0.990814  0.991289  0.991675  0.992083  0.992535  0.992876  0.993231  0.99361  0.99402  0.994146  0.994903  0.99506  0.995384  0.995717  0.996489  0.996583  0.99704  0.997304  0.997796  0.998238  0.998767  0.999368  0.999836  1];

y0(1) = 0; 
modelfun=@(t,x)nmodel(t,x,a,b,c);
[t ycalc]=ode45(modelfun,mytime,y0);
resid = (ycalc-mydata).*(ycalc-mydata);

plot(mytime,mydata,'k');
hold on
plot(t,ycalc,'ro');
hold off
val = sum(sum(resid));

%main
clc
 clear all
 global a;
 global b;
 global c;
a = 1; b = 1; c = 5;
nlstsq([a b c])

k=[a b c];
format long e
options = optimset('TolFun',1e-2,'TolX',1e-2,'MaxFunEvals',500,'MaxIter',500,'Display','Iter');
%options=optimset('MaxFunEvals',1000,'MaxIter',1000,'Display','Iter');
[k nlstsq]=fminsearch('nlstsq',k,options )

[Matt J]

You need to format your code (use the "{} Code" toolbar icon) to make it readable.

In the meantime, if the curve fit is okay but the parameters are different than you expect, it may mean that there are multiple solutions, (or approximately so).

[Abex]

Dear Matt, Many thanks. How can I access all possible solutions? In addition why is it so highly dependet on the initial values that I put?

[Ryan G]

I understand you may be limited in tools, but this is the problem simulink design optimization was created for. Essentially you could put your differential equation in Simulink and use SDO to solve for the 3 parameters. It's a bit easier to do visually than programtically in my opinion.

[Abex]

Can any one suggest me how to pass additional arguments over fminsearch to obtain optimum parameters which can minimise the objective function globally?

[Matt J]

FMINSEARCH is not a Global Optimization Toolbox function. Do you have that toolbox? If you use a local solver like fminsearch, your ability to minimize globally will inevitably depend on the quality of your single initial guess.

[Abex]

Dear Matt, Thank you. I have no Global Optimization Toolbox function and that is why I am wondering around. Can you suggest me the smart way of selecting initial guess? Is that good idea to take initial values close to the real value of the parameters?

Answer 1

[Matt J]

You should tune your TolFun and TolX parameter (1e-2 looks pretty big) by running some tests with simulated data where you know the true a,b,c parameters. Start with a moderately close initial guess [a0,b0,c0] to the known solution and make TolFun,TolX smaller and smaller until you get a good approximation of the solution.

Answer 2

[Matt J]

Can you suggest me the smart way of selecting initial guess? Is that good idea to take initial values close to the real value of the parameters?

Yes, choosing an initial guess close to the real values is generally thought to increase chances of avoiding local minima, as well as to speed up the optimization. My only idea for coming up with an initial guess is to ignore the errors in mydata and approximate log(K) as

 log(K) = log(diff(mydata)./diff(mytime)./(1-mydata))

You might perhaps want to filter/smooth mydata to get rid of noise content. Then try to get an approximate algebraic solution to the equation

 log(K)=log(b)-c/(0.0083144621.*(423+((10/60).*(t-a))));

for the unknowns a,c, and log(b). Note that log(K) will have a peak/valley at t=a, depending on whether c is positive or negative.. You could therefore look for peaks in log(K) to get an initial guess of "a". Once you have a value for "a", the equations reduce to a linear equation system in the remaining unknowns c and log(b), and you can use a linear solver to estimate them.

Another possibility would be to use

http://www.mathworks.com/matlabcentral/fileexchange/10093-fminspleas

which is able to process linear parameters like c and log(b) separately from nonlinear parameters like "a".