estimation of parameters using lsqnonlin

hello all, i am in the process of estimating parameters for my model

i wrote the following sample codes..pls tell me whether they are right. i ran the script i got the result as

Initial point is a local minimum.

Optimization completed because the size of the gradient at the initial point

is less than 1e-4 times the default value of the function tolerance.

Pls let me know whether my approach is right or any changes to be made....pls help me i am new to matlab...

1. My model :

function dc = kinetics(p,t,c)

k1=p(1);

k2=p(2);

k3=p(3);

k4=p(4);

k5=p(5);

k6=p(6);

dc =zeros(5,1);

dcdt(1) = (k1+k2+k3)*c(1);

dcdt(2) = (k1*c(1))-(k4*c(2));

dcdt(3) = (k2*c(1))-(k5*c(3));

dcdt(4) = (k3*c(1))-(k6*c(4));

dcdt(5)= k6*c(4);

end

2. my objective function to estimate parameters :

function obj= interest (p,t,c1meas,c2meas,c3meas,c4meas,c5meas)

global c1meas c2meas c3meas c4meas c5meas

cdata =[c1meas ,c2meas,c3meas,c4meas,c5meas]

%simulate model

c=simulate(p);

obj = c-cdata;

end

3. Simulate function containing ode45 to solve differential equations

function c = simulate(p)

global t c1meas c2meas c3meas c4meas c5meas

c=zeros(length(t),5);

c0=[1,0,0,0,0];

for i=1:length(t)-1

ts=[t(i),t(i+1)];

[tsol,csol]= ode45(@(t,c)kinetics2(p,t,c),ts,c0);

c(i+1,1)=c0(1);

c(i+1,2)=c0(2);

c(i+1,3)=c0(3);

c(i+1,4)=c0(4);

c(i+1,5)=c0(5);

end

end

4. my final script for calling all the above functions and using lsqnonlin

global t c1meas c2meas c3meas c4meas c5meas

t= [0.88 ; 0.96; 0.98;1.04 ; 1.05];

c1meas=[0.211 ;0.066 ;0.17 ; 0.455; 0.088];

c2meas =[0.666 ;0.165 ;0.083 ;0.047 ;0.009];

c3meas=[0.302 ;0.093; 0.155; 0.341 ;0.094];

c4meas=[0.237;0.084; 0.177; 0.404 ;0.082];

c5meas=[0.686 ;0.16 ;0.072; 0.041; 0.008];

%parameters initial guess

k1=0;

k2=0;

k3=1;

k4=1;

k5=1;

k6=1;

p0=[k1,k2,k3,k4,k5,k6];

p= lsqnonlin(@(p) interest(p),p0)

% show final objective

disp(['Final SSE Objective: ' num2str(interest(p))]);

% optimized parameter values

k1=p(1);

k2=p(2);

k3=p(3);

k4=p(4);

k5=p(5);

k6=p(6);

disp(['k1: ' num2str(k1)])

disp(['k2: ' num2str(k2)])

disp(['k3: ' num2str(k3)])

disp(['k4: ' num2str(k4)])

disp(['k5: ' num2str(k5)])

disp(['k6: ' num2str(k6)])

% calculate model with updated parameters

ci = simulate(p0);

cp = simulate(p);

Answer 1

[Torsten]

1.

p= lsqnonlin(@(p) interest(p),p0)

Use

p=lsqnonlin(@(p)interest(p,t,c1meas,c2meas,c3meas,c4meas,c5meas),p0)

instead.

2.

In "simulate", you have to assign csol to c, not c0.

[vinutha paravastu]

Thank you Torsten

1. i will modify isononlin function based on your suggestion.

2. in simulate function as i am using ode45 i have to use intital condition co ??

as ode45 is in the for loop and i wrote a code such that co value changes

[Torsten]

function csol = simulate(p)
global t 
c0 = [1,0,0,0,0];
[~,csol] = ode45(@(tf,c)kinetics2(p,tf,c),t,c0);
end

[vinutha paravastu]

thanku...

i will try...