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

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.

Answer 1

[Alan Weiss]

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

[Bosnian Kingdom]

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.