estimation of parameters using lsqnonlin
2 ビュー (過去 30 日間)
古いコメントを表示
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...
- 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 件のコメント
Alan Weiss
2019 年 8 月 26 日
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.
Alan Weiss
MATLAB mathematical toolbox documentation
回答 (0 件)
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!