MATLAB Answers

0

estimation of parameters using lsqnonlin

Asked by vinutha paravastu on 26 Aug 2019
Latest activity Commented on by vinutha paravastu on 26 Aug 2019
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);

  0 Comments

Sign in to comment.

Tags

1 Answer

Answer by Torsten
on 26 Aug 2019
Edited by Torsten
on 26 Aug 2019
 Accepted Answer

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.

  3 Comments

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
on 26 Aug 2019
function csol = simulate(p)
global t
c0 = [1,0,0,0,0];
[~,csol] = ode45(@(tf,c)kinetics2(p,tf,c),t,c0);
end
thanku...
i will try...

Sign in to comment.