MATLAB Answers

estimation of parameters using lsqnonlin

4 views (last 30 days)
vinutha paravastu
vinutha paravastu on 25 Aug 2019
Commented: Alan Weiss 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);

  1 Comment

Alan Weiss
Alan Weiss on 26 Aug 2019
I suspect that you are not passing t properly; it is declared as global in some places, but you do not declare it global in the kinetics function or interest function.
Generally, you should avoid using global variables. Learn to parametrize functions.
Alan Weiss
MATLAB mathematical toolbox documentation

Sign in to comment.

Answers (0)

Tags


Translated by