Is it posible to optimize kinetic parameter in following ordinary differential equation?

Question

Bosnian Kingdom 2019-4-16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/456703-is-it-posible-to-optimize-kinetic-parameter-in-following-ordinary-differential-equation

编辑： Bosnian Kingdom 2019-4-16

采纳的回答： Alan Weiss

在 MATLAB Online 中打开

I have temperature along reactor length in dependence of tubular reactor volume (experimental data).

I need to optimize kinetic parameters (k1, k2, k3, k4, k5, k6 and k7) in (-rA1), (-rA2) and (-rA3).

rA1 = ((k1*k2 * Ca * (Cb ^ k3)) / (1 +k2 * Ca)). In reaction kinetic (-rA2) and (-rA3) i need to optimize and k4, k5, k6 and k7.

I know U, a, Ta, Fa, Fc, Fb, Fe, Fd, Fe, -deltaHr1, -deltaHr2, -deltaHr3, cpA, cpB, cpC, cpD and cpE.

d(T)/d(V) = (((((U * a * (Ta - T) ))) + ((-rA1) * (-deltaHr1) + (-rA2) * (-deltaHr2) + (-rC3) * (-deltaHr3)))) / (FA * cpA + FB * cpB + FC * cpC + FD * cpD + FE * cpE)

First

I load data T=f(V)

Then I wrote function: function [dTdV] = dTdV(C) where C are kinetic parameters k1= C(1), k2=C(2) etc. in which I defined all known variables and wrote the equation dT / dV = ...

then i made script C=fminsearch(@dTdV,[0.1 0.1 0.1 0.001 0.12 0.20 0.01]);

and then I stucked.

Is it possible to solve this and in what way? Do I have to use Euler's method and how?

I tried to use Euler's method in the script, did I write well?

h=0.0001;
V=0:h:0.001177;
f=zeros(size(V));
T0=431.15;
f(1)=431.15;
n=numel(f);
for i=1:n-1
    f=((U*a.......);
    f(i+1)=f*h+f(i);
    
en

I got this error:

Index exceeds matrix dimensions.
Error in Script (line 64)
    f(i+1)=f*h+f(i);

Thank You in advance.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Alan Weiss 2019-4-16

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/456703-is-it-posible-to-optimize-kinetic-parameter-in-following-ordinary-differential-equation#answer_370878

I think that you would do better to use ode45 to solve your ODE, and use lsqcurvefit to optimize your parameters, as in this example.

But if you really want to do it your way, then I think the error is that you specify f instead of f(i) in this line:

f=((U*a.......);
    % Should be
f(i) = ((U*a.......);

Alan Weiss

MATLAB mathematical toolbox documentation

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Bosnian Kingdom 2019-4-16

编辑：Bosnian Kingdom 2019-4-16

在 MATLAB Online 中打开

Thank you for your answer. My codes are:

BEGIN OF FIRST CODE

function [dTdV] = dTdV(C)

U=200;

a=50;

Ta=578;

P0=101;

P=33.8;

T0=303;

FA=0.001678;

FB=0.01998;

FC=0.0001801;

FD=0.0018;

FE=0.00041;

T=680.15;

Ft=FA+FB+FC+FD+FE;

Ct0=P0/(8.314*T0);

Ca=Ct0 * (FA / Ft) * (T0 / T) * (P / P0);

Cb= Ct0 * (FB / Ft) * (T0 / T) * (P / P0);

Cc= Ct0 * (FC / Ft) * (T0 / T) * (P / P0);

Cd= Ct0 * (FD / Ft) * (T0 / T) * (P / P0);

Ce= Ct0 * (FE / Ft) * (T0 / T) * (P / P0);

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;

cpA = cpA1 + cpA2 * T + cpA3 * T ^ 2 + cpA4 * (T ^ 3);

cpB = cpB1 + cpB2 * T + cpB3 * T ^ 2 + cpB4 * (T ^ 3);

cpC = cpC1 + cpC2 * T + cpC3 * T ^ 2 + cpC4 * (T ^ 3);

cpD = cpD1 + cpD2 * T + cpD3 * T ^ 2 + cpD4 * (T ^ 3);

cpE = cpE1 + cpE2 * T + cpE3 * T ^ 2 + cpE4 * (T ^ 3);

ra1 = ((C(1) * C(2) * Ca * (Cb^C(3))) / (1 + C(2) * Ca))*0.001*3600;

ra2 = (C(4) * (Cb ^ C(3)))*0.001*3600;

rc3 = (C(5) * Cc * ((Cb ^ C(6)) / (Ca ^ C(7))))*0.001*3600;

dHr1 = -1012232 + (7.042 * (T - 298.15) + (0.019872 / 2) * (T ^ 2 - 298.15 ^ 2) - (0.0001475 / 3) * (T ^ 3 - 298.15 ^ 3) + (8.0205e-8 / 4) * (T ^ 4 - 298.15 ^ 4));

dHr2 = -2875895.41 + (45.985 * (T - 298.15) - (0.02312 / 2) * (T ^ 2 - 298.15 ^ 2) - (0.0001625 / 3) * (T ^ 3 - 298.15 ^ 3) + (1.2308e-7 / 4) * (T ^ 4 - 298.15 ^ 4));

dHr3 = -1398585.21 + (39.78 * (T - 298.15) - (0.025124 / 2) * (T ^ 2 - 298.15 ^ 2) - (2.09222 / 3) * (T ^ 3 - 298.15 ^ 3) + (2.22171E-8 / 4) * (T ^ 4 - 298.15 ^ 4));

f=((U*a*(Ta-T)+((-dHr1)*(-((C(1) * C(2) * Ca * (Cb^C(3))) / (1 + C(2) * Ca))*0.001*3600)+(-dHr2)*(-(C(4) * (Cb ^ C(3))))+(-dHr3)*(-(C(5) * Cc * ((Cb ^ C(6)) / (Ca ^ C(7)))))))/(FA*cpA+FB*cpB+FC*cpC+FD*cpD+FE*cpE));

dTdV=(f-T).^2;

dTdV=sum(dTdV);

end

END OF FIRST CODE

BEGIN OD SECOND CODE

load('data.mat');
C=fminsearch(@dTdV,[0.1 0.1 0.1 0.001 0.12 0.20 0.01]);

U=200;

a=50;

Ta=578;

P0=101;

P=33.8;

T0=303;

FA=0.001678;

FB=0.01998;

FC=0.0001801;

FD=0.0018;

FE=0.00041;

T=680.15;

Ft=FA+FB+FC+FD+FE;

Ct0=P0/(8.314*T0);

Ca=Ct0 * (FA / Ft) * (T0 / T) * (P / P0);

Cb= Ct0 * (FB / Ft) * (T0 / T) * (P / P0);

Cc= Ct0 * (FC / Ft) * (T0 / T) * (P / P0);

Cd= Ct0 * (FD / Ft) * (T0 / T) * (P / P0);

Ce= Ct0 * (FE / Ft) * (T0 / T) * (P / P0);

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;

cpA = cpA1 + cpA2 * T + cpA3 * T ^ 2 + cpA4 * (T ^ 3);

cpB = cpB1 + cpB2 * T + cpB3 * T ^ 2 + cpB4 * (T ^ 3);

cpC = cpC1 + cpC2 * T + cpC3 * T ^ 2 + cpC4 * (T ^ 3);

cpD = cpD1 + cpD2 * T + cpD3 * T ^ 2 + cpD4 * (T ^ 3);

cpE = cpE1 + cpE2 * T + cpE3 * T ^ 2 + cpE4 * (T ^ 3);

ra1 = ((C(1) * C(2) * Ca * (Cb^C(3))) / (1 + C(2) * Ca))*0.001*3600;

ra2 = (C(4) * (Cb ^ C(3)))*0.001*3600;

rc3 = (C(5) * Cc * ((Cb ^ C(6)) / (Ca ^ C(7))))*0.001*3600;

dHr1 = -1012232 + (7.042 * (T - 298.15) + (0.019872 / 2) * (T ^ 2 - 298.15 ^ 2) - (0.0001475 / 3) * (T ^ 3 - 298.15 ^ 3) + (8.0205e-8 / 4) * (T ^ 4 - 298.15 ^ 4));

dHr2 = -2875895.41 + (45.985 * (T - 298.15) - (0.02312 / 2) * (T ^ 2 - 298.15 ^ 2) - (0.0001625 / 3) * (T ^ 3 - 298.15 ^ 3) + (1.2308e-7 / 4) * (T ^ 4 - 298.15 ^ 4));

dHr3 = -1398585.21 + (39.78 * (T - 298.15) - (0.025124 / 2) * (T ^ 2 - 298.15 ^ 2) - (2.09222 / 3) * (T ^ 3 - 298.15 ^ 3) + (2.22171E-8 / 4) * (T ^ 4 - 298.15 ^ 4));

h=0.0001;
V=0:h:0.001177;
f=zeros(size(V));
T0=431.15;
f(1)=431.15;
n=numel(f);
for i=1:n-1
    f=((U*a*(Ta-T)+((-dHr1)*(-((C(1) * C(2) * Ca * (Cb^C(3))) / (1 + C(2) * Ca))*0.001*3600)+(-dHr2)*(-(C(4) * (Cb ^ C(3)))*0.001*3600)+(-dHr3)*(-(C(5) * Cc * ((Cb ^ C(6)) / (Ca ^ C(7))))*0.001*3600)))/(FA*cpA+FB*cpB+FC*cpC+FD*cpD+FE*cpE));
    f(i+1)=f*h+f(i);
    
end
load('data.mat');
figure(1);
hold on;
plot(V,T);
fplot(@(V)f*h+f(i), [0 30]);

END OF SECOND CODE

And T=f(V) is:

T V

303 0

680.15 0.000040474

680.15 7.50925E-05

681.15 0.000109711

678.15 0.00014433

680.15 0.000178948

680.15 0.000213567

681.15 0.000248185

680.15 0.000282804

680.15 0.000317422

679.15 0.000317422

681.15 0.000352041

679.15 0.000352041

679.15 0.000386659

680.15 0.000386659

681.15 0.000421278

680.15 0.000455896

679.15 0.000490515

680.15 0.000525133

678.15 0.000559752

681.15 0.00059437

680.15 0.000628989

680.15 0.000663607

680.15 0.000698226

681.15 0.000732844

681.15 0.000767463

678.15 0.000802081

681.15 0.000871318

681.15 0.000940555

678.15 0.001009792

680.15 0.001079029

Thanks in advance.

请先登录，再进行评论。

Is it posible to optimize kinetic parameter in following ordinary differential equation?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Is it posible to optimize kinetic parameter in following ordinary differential equation?

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论