5 views (last 30 days)

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.

Stephane Dauvillier
on 9 May 2019

Sorry, I may have been confusing. The LB and UB will be the lower and upper bound of your parameter. Youtr parameter is a 10 elements vector so LB and UB should have 1à elements (that's why you should have a warning message saying

Warning: Length of lower bounds is < length(x); filling in missing lower bounds with -Inf.

I was confused by the 2 variables named both var.

Knowing that, can you set boundaries or do you need to specify constraint (see the docuement page here) ?

For your problem, the column is the "Least Squares" but depending on your constraint you lay have to change the solver.

Stephane Dauvillier
on 9 May 2019

Hi it's seems that the variable U is the first to become complex (just put a breakpoint on the following line):

U=1.22*((var(8)-Ta)/du)^(1/4);

In fact on first iteration (var(8)-Ta)/du is negative, which explains why you have a complex number for the variable U

You may want to specify boundaries on var to avoid this

Stephane Dauvillier
on 9 May 2019

Hi, You can specify lower and upper bounds for your var.

By the way is it normal that your function DifEq compute the following variable but don't use it :

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);

Stephane Dauvillier
on 9 May 2019

Yes, except you have to modify LB and UB to specify these boundaries. These should be vectors (for instance LB(1) and UB(1) will be the lower and upper boundaries for var(1) and so on).

If one of the element has no lower or upper boundaries, just set it to -inf or +inf

Stephane Dauvillier
on 9 May 2019

Have follow the hyperlink:

Knowing that, can you set boundaries or do you need to specify constraint (see the docuement page here) ?

Your objective function is least square since your want to minimize the square root norm.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 16 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_702499

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_702499

## 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_702796

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_702796

## 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_702867

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_702867

## 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_703090

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_703090

## 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_703109

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_703109

## 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_703319

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_703319

## 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_703360

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_703360

## 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_703384

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/460550-why-i-getting-complex-double-in-fitting-ode-with-experimental-data#comment_703384

Sign in to comment.