MATLAB Answers

can not find parameter estimation

2 views (last 30 days)
Mouna YAHYA
Mouna YAHYA on 1 Sep 2020
Hi
I am working on the following code
function Mouna
function dydt=kinetics(p,t,y)
y0=[26 16.4 14.6 0.069 0.84 1.1 0.05 12.1 12.1 0.48 0.2 520 1 0.21 0.1];
[t,y]=ode45(@kinetics,t,y0);
%
dydt=zeros(15,1); %initialize
%Hydrogen production
%growth coefficient
u1=p1*y(2)/(p2+y(2));
%Differential equations
dydt(1)= -p(3)*y(1)-p(4)*y(3)*y(1)+p(3)*p(5)*p(6); %S0 g/L
dydt(2)= -p(3)*y(2)+p(4)*y(3)*y(1)-u1*y(3)/p(7); %S1 g/L
dydt(3)=u1*y(3)-p(3)*y(3); %X1 g/L
dydt(4)=u1*y(3)/p(8)-p(3)*y(4); %Pr1 g/L
dydt(5)=u1*y(3)/p(9)-p(3)*y(5); %But1 g/L
dydt(6)=u1*y(3)/p(10)-p(3)*y(6); %Ac1 g/L
%Hydrogen rate
y(7)=p(11)*u1*y(3); %H2 g/L
%**********************************************************************************************************************
%Methane production
%kinetic equations
ufBut2=p(12)*y(11)/(p(13)+y(11));
ufPr2=p(14)*y(10)/(p(15)+y(10));
qfBut2=p(16)*y(11)/(p(13)+y(11));
qfPr2=p(17)*y(10)/(p(15)+y(10));
um=p(18)*y(14)/(p(19)+y(14));
qm=p(20)*y(14)/(p(21)+y(14));
ue=p(22)*y(14)/(p(23)+y(14));
qe=p(24)*y(14)/(p(23)+y(14));
%VFA differential equations
dydt(8)=ufPr2*y(8)-p(25)*y(8)-p(26)*y(8); %XfPr2
dydt(9)=ufBut2*y(9)-p(25)*y(9)-p(26)*y(9); %XfBut2
dydt(9)= -qfPr2*y(8)+p(27)*(y(5)-y(10)); %Pr2
dydt(9)= -qfBut2*y(9)+p(27)*(y(6)-y(11)); %But2
dydt(12)=um*y(12)-p(28)*y(12)-p(29)*y(12); %Xm
dydt(13)=ue*y(13)-p(30)*y(13)-p(31)*y(13); %Xe
dydt(14)= -qe*y(14)-qm*y(12)+p(27)*(y(7)-y(14)); %Ac2
%Methane rate
y(15)=p(32)*qm*(y(13)+y(14))*p(33);
end
t=[26 16.4 14.6 0.069 0.84 1.1 0.05 12.1 12.1 0.4 0.2 520 1 0.21 0.1];
yexp=[26 16.4 12.4 0.118 1.888 2.64 1.64 12.1 12.1 0.55 0.3 12.1 12.1 0.2 2.58
2402 16.02 11.87 0.145 1.678 2.53 1.6 11.56 11.56 0.24 0.4 11.56 11.56 0.195 2.58
20.8 14.89 11.02 0.145 1.972 2.75 1.63 10.72 10.72 0.34 0.15 10.72 10.72 0.211 2.46
18 13.34 10.59 0.152 1.982 2.725 1.68 10.01 10.01 0.22 0.12 10.01 10.01 0.19 2.52
17.65 11.79 9.63 0.138 1.89 2.73 1.7 9.32 9.32 0.28 0.11 9.32 9.32 0.216 2.56
16.98 9.8 8.88 0.145 1.806 2.65 1.7 8.67 8.67 0.24 0.14 8.67 8.67 0.22 2.54
15.4 7.8 8 0.143 1.982 2.77 1.69 7.85 7.85 0.2 0.13 7.85 7.85 0.22 2.65
14.2 6.2 7.21 0.16 1.962 2.79 1.65 6.95 6.95 0.2 0.15 6.95 6.95 0.221 2.62
13 4.9 6.53 0.163 2.184 2.86 1.8 6.02 6.02 0.2 0.27 6.02 6.02 0.288 2.61
10.4 3.2 5.03 0.173 2.037 284 1.81 5.34 5.34 0.34 0.24 5.34 5.34 0.293 2.6];
p0=[0.586;3.914;3.914;1;1;16.4;0.08;4.2;2.1;1.2;0.05;0.05;0.22;0.05;0.22;12;14;0.3;80;0.05;20;0.05;20;0.04;0.1;4.12;0.01;0.2;0.01;0.2;122;2.2;2.2];
p=lsqcurvefit(@kinetics,p0,t,yexp);
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(p)
fprintf(1, '\t\tp(%d) = %8.5f\n', k1, p(k1))
end
tv = linspace(min(t), max(t));
yfit = kinetics(p, tv);
end
I got
Caused by:
Failure in initial user-supplied objective function evaluation.
LSQCURVEFIT cannot continue.

  0 Comments

Sign in to comment.

Answers (0)

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!