MATLAB Answers

PARAMETER estimation of kinetic and adsorbtion constant of langmuir hinshelwood haugen watson model

10 ビュー (過去 30 日間)
GOVIND BHANDARI
GOVIND BHANDARI 2021 年 9 月 28 日
回答済み: Alex Sha 2021 年 10 月 19 日
please help me calculated parameter for different model atttached in literature. i tried for model E but could not get my parameter values close to anywhere mentioned in literature. Ch2 is consant = 2.47e-3
star strider please help me

回答 (2 件)

GOVIND BHANDARI
GOVIND BHANDARI 2021 年 10 月 4 日
function model4
% 2016 12 03
% NOTES:
%
% 1. The ‘theta’ (parameter) argument has to be first in your
% ‘kinetics’ funciton,
% 2. You need to return ALL the values from ‘DifEq’ since you are fitting
% all the values
function C=kinetics(theta,t)
c0=[7.00E-02;0.00E+00;0.00E+00;0.00E+00];
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
dcdt=zeros(4,1);
dcdt(1)=(-theta(1)*c(1)*0.00273-theta(2)*c(1)*0.00273)/(1+(theta(6)*c(1)+theta(7)*c(2)+theta(8)*c(3)+theta(9)*c(4)+theta(10)*0.00273))^3;
dcdt(2)= (theta(1)*c(1)*0.00273-theta(3).*c(2)*0.00273-theta(5)*c(2))/(1+(theta(6)*c(1)+theta(7)*c(2)+theta(8)*c(3)+theta(9)*c(4)+theta(10)*0.00273))^3;
dcdt(3)= (theta(2)*c(1)*0.00273-theta(4)*c(2)*0.00273+theta(5)*c(2))/(1+(theta(6)*c(1)+theta(7)*c(2)+theta(8)*c(3)+theta(9)*c(4)+theta(10)*0.00273))^3;
dcdt(4)= (theta(3)*c(2)*0.00273+theta(4)*c(3)*0.00273)/(1+(theta(6)*c(1)+theta(7)*c(2)+theta(8)*c(3)+theta(9)*c(4)+theta(10)*0.00273))^3;
dC=dcdt;
end
C=Cv;
end
t=[20
40
60
90
120
180
240];
c=[5.69E-02 1.36E-03 8.25E-03 3.48E-03
4.28E-02 9.79E-03 9.50E-03 7.94E-03
2.93E-02 2.85E-03 1.53E-02 2.26E-02
1.46E-02 4.69E-03 2.00E-02 3.28E-02
1.19E-03 3.52E-03 1.10E-02 5.43E-02
1.75E-03 0.00E+00 0.00E+00 6.10E-02
1.00E-09 1.00E-09 1.00E-09 7.00E-02];
theta0=[1.86e-1,1.70e-1,1.28,1.51e-1,1.0e-1,15.90,2.60e-1,1.87e-2,3.25,2.33e-2];
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);
figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit);
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'Location','N')
end
THE PROBLEM IS ORANGE STAR AND RED LINE ARE NOT FITTING THE TOP SHOULD BE AROUND 0.02

Alex Sha
Alex Sha 2021 年 10 月 19 日
if want all parameters to be positive:
Root of Mean Square Error (RMSE): 0.0035717712315911
Sum of Squared Residual: 0.00035721139246301
Correlation Coef. (R): 0.892442040492228
R-Square: 0.796452795637932
Parameter Best Estimate
-------------------- -------------
theta1 17.172493498211
theta2 21.5160808499238
theta3 38.9519219861006
theta4 51.8990082731777
theta5 5.17756884421286E-16
theta6 17.526757560822
theta7 6.33562009831634E-21
theta8 1.27030055679559E-14
theta9 9.92523484202955
theta10 11.9701814364323
while, if all parameters are free ranges, the solutions will be not unique, one of them like below:
Root of Mean Square Error (RMSE): 0.00271818383364327
Sum of Squared Residual: 0.000206878653897429
Correlation Coef. (R): 0.894998138423552
R-Square: 0.801021667781624
Parameter Best Estimate
-------------------- -------------
theta1 0.306719359089766
theta2 0.808060248987747
theta3 4.89340760623797
theta4 -0.0316673300900705
theta5 -0.00965914856097772
theta6 -79.2660953582632
theta7 -52.8331941956074
theta8 -94.6702336579354
theta9 -83.3482760040795
theta10 1914.31778883528

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by