I need to optimize kinetic parameter in following equations:
dca/dw=rA
dcb/dw=rB
dcc/dw=rC
dcd/dw=rD
dce/dw=rE
dcg/dw=rG
dP/dw=(-alfa / 2) * (var(8) / T0) * (var(7) / (var(7) / P0)) * (Ft / Ft0);
dT/dW=(((U*a*(Ta-var(8))/roB*10)+((-dHr1)*(-rA1)+(-dHr2)*(-rA2)+(-dHr3)*(-rC3)+(-dHr4)*(-rA4))))/(var(1)*cpA+var(2)*cpB+var(3)*cpC+var(4)*cpD+var(5)*cpE+var(6)*cpG);
where:
rA = -(rA1 + rA2 + rA4);
rB1 = 3.5 * rA1;
rB2 = 6.5 * rA2;
rB3 = 3 * rC3;
rB4 = 15 / 6 * rA4;
rB = -(rB1 + rB2 + rB3 + rB4);
rC1 =rA1;
rC = rC1+rC3;
rD1 = 4 * rA1;
rD2 = 5 * rA2;
rD3 = rC3;
rD4 = 14 / 6 * rA4;
rD = rD1 + rD2 + rD3 + rD4;
rE2 =4 * rA2;
rE3 =4 * rC3;
rE = rE2+rE3;
rG =8 / 6 * rA4;
I used this:
and this:
to make function:
function C=kinetics(par,W)
var0=[0.001808;0.02258704;0;0.00057911;0.00003447;0;134;431.15];
[W,Cv]=ode45(@DifEq,W,var0);
function dC=DifEq(W,var)
dcdt=zeros(8,1);
T0=431.15;
a=190.48;
Ta=673.15;
du=0.021;
U=1.22*((var(8)-Ta)/du)^(1/4);
P0=134;
pi = 3.14;
dv = 25 / 1000;
delta = 2 / 1000;
Hkat = 3.25;
roB = 0.965517241/ ((pi / 4) * ((dv - 2 * delta) ^ 2 * Hkat));
Ft=var(1)+var(2)+var(3)+var(4)+var(5)+var(6);
Ct0=P0/(8.314*T0);
Ca= Ct0 * (var(1) / Ft) * (T0 / var(8)) * (var(7) / P0);
Cb= Ct0 * (var(2) / Ft) * (T0 / var(8)) * (var(7) / P0);
Cc= Ct0 * (var(3) / Ft) * (T0 / var(8)) * (var(7) / P0);
Cd= Ct0 * (var(4) / Ft) * (T0 / var(8)) * (var(7) / P0);
Ce= Ct0 * (var(5) / Ft) * (T0 / var(8)) * (var(7) / P0);
Cg= Ct0 * (var(6) / Ft) * (T0 / var(8)) * (var(7) / P0);
FA0= 0.001808;
FB0= 0.02258704;
FC0= 0;
FD0= 0.00057911;
FE0= 0.00003447;
FG0 = 0;
Ft0=FA0+FB0+FC0+FD0+FE0+FG0;
G = 0.0007189;
epsilon = 0.45;
roO = 1.160147;
Dp = 0.003461538;
Ac = 0.00035;
viskA = ((261.5278 + 12.55191 * (var(7) / 100) - 0.55997 * (var(8) - 273.15)) / 1000000) * 3600;
yA = 0.016405;
Ma = 58;
viskB = ((17.87683461 + 4.238289936 * (var(7) / 100) + 0.016904644 * (var(8) - 273.15)) / 1000000) * 3600;
yB = 0.773086;
Mb = 28;
viskC = ((28.30049 + 5.541702 * (var(7) / 100) - 0.00358 * (var(8) - 273.15)) / 1000000) * 3600;
yC = 0.000313;
Mc = 44;
viskE = ((730.3278 + 11.11223 * (var(7) / 100) - 1.11653 * (var(8) - 273.15)) / 1000000) * 3600;
yE = 0.005255;
Me = 18;
viskD = ((-5.38734 + 2.683369 * (var(7) / 100) + 0.877204 * (var(8) - 273.15)) / 1000000) * 3600;
yD = 0.204942;
Md = 32;
viskosmj = (((viskA * yA * (Ma ^ 0.5)) + (viskB * yB * (Mb ^ 0.5)) + (viskC * yC * (Mc ^ 0.5)) + (viskD * yD * (Md ^ 0.5)) + (viskE * yE * (Me ^ 0.5))) / ((yA * (Ma ^ 0.5)) + (yB * (Mb ^ 0.5)) + (yC * (Mc ^ 0.5)) + (yD * (Md ^ 0.5)) + (yE * (Me ^ 0.5)))) / 1000;
alfa = ((2 * G * (1 - epsilon)) / (roO * Dp * epsilon ^ 3 * roB * Ac * P0)) * ((150 * (1 - epsilon) * viskosmj) / Dp + 1.75 * G) * (760 / (60 * 10));
cpA1 = 9.487;
cpA2 = 3.313e-1;
cpA3 = -0.0001108;
cpA4 = -0.000000002821;
cpB1 = 28.106;
cpB2 = -0.00000368;
cpB3 = 1.475e-5;
cpB4 = -0.00000001065;
cpC1 = -13.075;
cpC2 = 0.3484;
cpC3 = -0.0002184;
cpC4 = 0.00000004839;
cpD1 = 32.243;
cpD2 = 0.001923;
cpD3 = 0.00001055;
cpD4 = -0.000000003596;
cpE1 = 19.975;
cpE2 = 0.07343;
cpE3 = -0.00005601;
cpE4 = 0.00000001715;
cpG1 = 1.742;
cpG2 = 3.19e-1;
cpG3 = -0.0002352;
cpG4 = 6.975e-8;
cpA = cpA1 + cpA2 * var(8) + cpA3 * var(8) ^ 2 + cpA4 * (var(8) ^ 3);
cpB = cpB1 + cpB2 * var(8) + cpB3 * var(8) ^ 2 + cpB4 * (var(8) ^ 3);
cpC = cpC1 + cpC2 * var(8) + cpC3 * var(8) ^ 2 + cpC4 * (var(8) ^ 3);
cpD = cpD1 + cpD2 * var(8) + cpD3 * var(8) ^ 2 + cpD4 * (var(8) ^ 3);
cpE = cpE1 + cpE2 * var(8) + cpE3 * var(8) ^ 2 + cpE4 * (var(8) ^ 3);
cpG = cpG1 + cpG2 * var(8) + cpG3 * var(8) ^ 2 + cpG4 * (var(8) ^ 3);
rA1=((par(1) * par(2) * Ca * (Cb^par(3)))) / (1 + par(2) * Ca);
rA2=(par(4) * (Cb ^ par(3)));
rC3=(par(5) * Cc * ((Cb ^ par(6)))/ (Ca ^ par(7)));
rA4 =((par(8) * par(9)*Ca* (Cb ^ par(10)))) / (1 + par(9)*Ca);
dHr1 = -1240070 + (8.039 * (var(8) - 298.15) + (0.024805 / 2) * (var(8) ^ 2 - 298.15 ^ 2) - (0.00012 / 3) * (var(8) ^ 3 - 298.15 ^ 3) + (7.4102e-8 / 4) * (var(8) ^ 4 - 298.15 ^ 4));
dHr2 = -2646628.37 + (48.939 * (var(8) - 298.15) - (0.0279411 / 2) * (var(8) ^ 2 - 298.15 ^ 2) - (0.0001551 / 3) * (var(8) ^ 3 - 298.15 ^ 3) + (1.22801e-7 / 4) * (var(8) ^ 4 - 298.15 ^ 4));
dHr3 = -1428408.55 + (40.9 * (var(8) - 298.15) - (0.026385 / 2) * (var(8) ^ 2 - 298.15 ^ 2) - (1.19333 / 3) * (var(8) ^ 3 - 298.15 ^ 3) + (1.2141E-8 / 4) * (var(8) ^ 4 - 298.15 ^ 4));
dHr4 = -5322.22+(1005.905*(var(8)-298.15)+4.595567/2*(var(8)^2-298.15^2)-0.00224/3*(var(8)^3-298.15^3)+4.37E-7/4*(var(8)^4-298.15^4));
rA = -(rA1 + rA2 + rA4);
rB1 = 3.5 * rA1;
rB2 = 6.5 * rA2;
rB3 = 3 * rC3;
rB4 = 15 / 6 * rA4;
rB = -(rB1 + rB2 + rB3 + rB4);
rC1 =rA1;
rC = rC1-rC3;
rD1 = 4 * rA1;
rD2 = 5 * rA2;
rD3 = rC3;
rD4 = 14 / 6 * rA4;
rD = rD1 + rD2 + rD3 + rD4;
rE2 =4 * rA2;
rE3 =4 * rC3;
rE = rE2+rE3;
rG =8 / 6 * rA4;
dcdt(1)=rA;
dcdt(2)=rB;
dcdt(3)=rC;
dcdt(4)=rD;
dcdt(5)=rE;
dcdt(6)=rG;
dcdt(7)=(-alfa / 2) * (var(8) / T0) * (var(7) / (var(7) / P0)) * (Ft / Ft0);
dcdt(8)=(((U*a*(Ta-var(8))/roB*10)+((-dHr1)*(-rA1)+(-dHr2)*(-rA2)+(-dHr3)*(-rC3)+(-dHr4)*(-rA4))))/(var(1)*cpA+var(2)*cpB+var(3)*cpC+var(4)*cpD+var(5)*cpE+var(6)*cpG);
dC=dcdt;
end
C=Cv;
end
and i made script:
W=[0
0.118832891
0.148541114
0.178249337
0.20795756
0.237665782
0.267374005
0.297082228
0.326790451
0.356498674
0.356498675
0.386206896
0.386206897
0.415915119
0.415915120
0.445623342
0.475331565
0.505039788
0.53474801
0.564456233
0.594164456
0.623872679
0.653580902
0.683289124
0.712997347
0.74270557
0.772413793
0.831830238
0.89124668
0.95066313
0.965517241];
var=[0.001808 0.022587 0 0.000579 0.00003447 0 134 431.15
0.0017568 0.0225492 0.000048 0.00162 0.0003141 0.000979 131.75 685.15
0.0017056 0.0225114 0.000096 0.002661 0.00059373 0.001958 129.5 685.15
0.0016544 0.0224736 0.000144 0.003702 0.00087336 0.002937 127.25 686.15
0.0016032 0.0224358 0.000192 0.004743 0.00115299 0.003916 125 683.15
0.001552 0.022398 0.00024 0.005784 0.00143262 0.004895 122.75 685.15
0.0015008 0.0223602 0.000288 0.006825 0.00171225 0.005874 120.5 685.15
0.0014496 0.0223224 0.000336 0.007866 0.00199188 0.006853 118.25 686.15
0.0013984 0.0222846 0.000384 0.008907 0.00227151 0.007832 116 685.15
0.0013472 0.0222468 0.000432 0.009948 0.00255114 0.008811 113.75 685.15
0.001296 0.022209 0.00048 0.010989 0.00283077 0.00979 111.5 684.15
0.0012448 0.0221712 0.000528 0.01203 0.0031104 0.010769 109.25 686.15
0.0011936 0.0221334 0.000576 0.013071 0.00339003 0.011748 107 684.15
0.0011424 0.0220956 0.000624 0.014112 0.00366966 0.012727 104.75 684.15
0.0010912 0.0220578 0.000672 0.015153 0.00394929 0.013706 102.5 685.15
0.00104 0.02202 0.00072 0.016194 0.00422892 0.014685 100.25 686.15
0.0009888 0.0219822 0.000768 0.017235 0.00450855 0.015664 98 685.15
0.0009376 0.0219444 0.000816 0.018276 0.00478818 0.016643 95.75 684.15
0.0008864 0.0219066 0.000864 0.019317 0.00506781 0.017622 93.5 685.15
0.0008352 0.0218688 0.000912 0.020358 0.00534744 0.018601 91.25 683.15
0.000784 0.021831 0.00096 0.021399 0.00562707 0.01958 89 686.15
0.0007328 0.0217932 0.001008 0.02244 0.0059067 0.020559 86.75 685.15
0.0006816 0.0217554 0.001056 0.023481 0.00618633 0.021538 84.5 685.15
0.0006304 0.0217176 0.001104 0.024522 0.00646596 0.022517 82.25 685.15
0.0005792 0.0216798 0.001152 0.025563 0.00674559 0.023496 80 686.15
0.000528 0.021642 0.0012 0.026604 0.00702522 0.024475 77.75 686.15
0.0004768 0.0216042 0.001248 0.027645 0.00730485 0.025454 75.5 683.15
0.0004256 0.0215664 0.001296 0.028686 0.00758448 0.026433 73.25 686.15
0.0003744 0.0215286 0.001344 0.029727 0.00786411 0.027412 71 686.15
0.0003232 0.0214908 0.001392 0.030768 0.00814374 0.028391 68.75 683.15
0.000272 0.021453 0.00144 0.031809 0.00842337 0.02937 66.5 685.15
];
par0=[0.0001;2250;0.2;0.00001;0.00001;0.1;1;0.001;1;1];
[par,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,par0,W,var);
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(par)
fprintf(1, '\t\tpar(%d) = %8.5f\n', k1, par(k1))
end
Wv = linspace(min(W), max(W));
Cfit = kinetics(par, Wv);
figure(1)
plot(W, var, 'p')
hold on
plot(Wv,Cfit);
hold off
grid
xlabel('Mass, kg')
ylabel('All, K')
I write t instead of w.
I got complex double, and i got only good agreement for 8th equation.
Please, try my MATLAB code and if you know what is problem, i will appreciate help.
Thanks in advance.
9 Comments
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_702436
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_702436
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_702711
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_702711
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_702827
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_702827
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_702936
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_702936
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_703102
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_703102
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_703138
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_703138
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_703325
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_703325
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_703372
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_703372
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_847216
Direct link to this comment
https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_847216
Sign in to comment.