I want to find two missing parameters in an ODE system of equations using regression/optimization in MATLAB

Question

naiva saeedia 2024-4-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2102151-i-want-to-find-two-missing-parameters-in-an-ode-system-of-equations-using-regression-optimization-in

编辑： Sam Chak 2024-4-15

Hi guys I have this code for simulation of a Reactor and I have a ODE system to solve. The problem is I don't have the values of k_L_a_G_L_Light and k_L_a_G_L_Heavy..I have experimental values for x(end,2),x(end,3),x(end,4),x(end,5) and x(end,6):

x(end,2)=0.07891

x(end,3)=0.14349

x(end,4)=0.063348

x(end,5)=0.02686

x(end,6)=0.01351

Is it possible for MATLAB to use these datas and fit a value for k_L_a_G_L_Light and k_L_a_G_L_Heavy by performing some kind of regression or optimization? Any help would be really appreciate becaue I'm not able to do it based on my knowledge in MATLAB.

clc

clear

close all

T_Ethylene=230+273.15;

T_ref=230+273.15;

R=8.314;

d_T=12.5e-2;

d_I=3.15e-2;

Q_G_in_ml_min=580; %mL/min

Q_G_in=Q_G_in_ml_min.*1.66667e-8; %m3/s

F_G_out_mmol_s=0.354;

F_G_out=F_G_out_mmol_s.*1e-3;

V_R=300e-6;

V_G=250e-6;

V_L=50e-6;

%Calculating D_I_M for all is

T=T_Ethylene;%K

M_B=28;%g/mol

M_1=M_B;

M_2=56;

M_3=84;

M_4=112;

M_5=140;

M_6=168;

M_7=156;

P=35.5292;%atm

v_C=15.9;%cm3

v_H=2.31;%cm3

D_I_M_1=(1e-3.*T.^1.75.*(1./M_1+1./M_B))./(P.*((2*v_C+4*v_H)^(1/3)+(2*v_C+4*v_H)^(1/3)));

D_I_M_2=(1e-3.*T.^1.75.*(1./M_2+1./M_B))./(P.*((4*v_C+8*v_H)^(1/3)+(2*v_C+4*v_H)^(1/3)));

D_I_M_3=(1e-3.*T.^1.75.*(1./M_3+1./M_B))./(P.*((6*v_C+12*v_H)^(1/3)+(2*v_C+4*v_H)^(1/3)));

D_I_M_4=(1e-3.*T.^1.75.*(1./M_4+1./M_B))./(P.*((8*v_C+16*v_H)^(1/3)+(2*v_C+4*v_H)^(1/3)));

D_I_M_5=(1e-3.*T.^1.75.*(1./M_5+1./M_B))./(P.*((10*v_C+20*v_H)^(1/3)+(2*v_C+4*v_H)^(1/3)));

D_I_M_6=(1e-3.*T.^1.75.*(1./M_6+1./M_B))./(P.*((12*v_C+24*v_H)^(1/3)+(2*v_C+4*v_H)^(1/3)));

D_I_M_7=(1e-3.*T.^1.75.*(1./M_7+1./M_B))./(P.*((11*v_C+12*v_H)^(1/3)+(2*v_C+4*v_H)^(1/3)));

D_o_w=0.3173;

N_cd=(4.*Q_G_in.^0.5.*d_T.^0.25)./d_I.^2;

N_I=89.01;

U_G=1;

k_L_a_o_w=3.35*(N_I/N_cd)^1.464*U_G;

Lambda_L=5.515e-4;

Delta_L=4.095;

Lambda_H=2.320e-6;

Delta_H=2.207;

%Calculating k_L_i

k_L_a_G_L_Light=0.055;

k_L_a_G_L_Heavy=0.029;

k1_0=2.224e-4;

k2_0=1.533e-4;

k3_0=3.988e-5;

k4_0=1.914e-7;

k5_0=4.328e-5;

k6_0=2.506e-7;

k7_0=4.036e-5;

k8_0=1.062e-6;

k9_0=6.055e-7;

E1=109.53;

E2=15.229;

E3=7.881;

E4=44.454;

E5=9.438;

E6=8.426;

E7=10.911;

E8=12.537;

E9=7.127;

k1=k1_0.*exp(-E1.*(1/T_Ethylene-1/T_ref)./R);

k2=k2_0.*exp(-E2.*(1/T_Ethylene-1/T_ref)./R);

k3=k3_0.*exp(-E3.*(1/T_Ethylene-1/T_ref)./R);

k4=k4_0.*exp(-E4.*(1/T_Ethylene-1/T_ref)./R);

k5=k5_0.*exp(-E5.*(1/T_Ethylene-1/T_ref)./R);

k6=k6_0.*exp(-E6.*(1/T_Ethylene-1/T_ref)./R);

k7=k7_0.*exp(-E7.*(1/T_Ethylene-1/T_ref)./R);

k8=k8_0.*exp(-E8.*(1/T_Ethylene-1/T_ref)./R);

k9=k9_0.*exp(-E9.*(1/T_Ethylene-1/T_ref)./R);

w_cat=(0.3+0.25)*1e-3;

%Ethylene:

P_Ethylene=36e5;

C_G_in=P_Ethylene/(R*T_Ethylene);

n_initial_ethylene=C_G_in.*V_R;

n_initial_solvent=0.2348;

Q_G_in_C_G_in_C4=0;

Q_G_in_C_G_in_C6=0;

Q_G_in_C_G_in_C8=0;

Q_G_in_C_G_in_C10=0;

Q_G_in_C_G_in_C12=0;

Q_G_in_C_G_in_C11H24=0;

K_L_G_Ethylene=3.24;

K_L_G_Butene=2.23;

K_L_G_Hexene=1.72;

K_L_G_Octene=1.3;

K_L_G_Decene=0.985;

K_L_G_Dodecene=0.751;

K_L_G_Undecane=0.842;

dNdt=@(t,x) [Q_G_in.*C_G_in-F_G_out-V_R.*k_L_a_G_L_Light.*(x(1)./V_G-K_L_G_Ethylene.*x(8)./V_L);Q_G_in_C_G_in_C4-F_G_out-V_R.*k_L_a_G_L_Light.*(x(2)./V_G-K_L_G_Butene.*(x(9)./V_L));Q_G_in_C_G_in_C6-F_G_out-V_R.*k_L_a_G_L_Light.*(x(3)./V_G-K_L_G_Hexene.*(x(10)./V_L));Q_G_in_C_G_in_C8-F_G_out-V_R.*k_L_a_G_L_Heavy.*(x(4)./V_G-K_L_G_Octene.*(x(11)./V_L));Q_G_in_C_G_in_C10-F_G_out-V_R.*k_L_a_G_L_Heavy.*(x(5)./V_G-K_L_G_Decene.*x(12)./V_L);Q_G_in_C_G_in_C12-F_G_out-V_R.*k_L_a_G_L_Heavy.*(x(6)./V_G-K_L_G_Dodecene.*x(13)./V_L);Q_G_in_C_G_in_C11H24-F_G_out-V_R.*k_L_a_G_L_Heavy.*(x(7)./V_G-K_L_G_Undecane.*x(14)./V_L);V_R.*k_L_a_G_L_Light.*(x(1)./V_G-K_L_G_Ethylene.*x(8)./V_L)+w_cat.*(-2.*k1.*(x(8)./V_L).^2-k2.*(x(8)./V_L).*(x(9)./V_L)-k3.*(x(8)./V_L).*(x(10)./V_L)-k5.*(x(8)./V_L).*(x(11)./V_L)-k7.*x(8).*x(12)./(V_L.^2));V_R.*k_L_a_G_L_Light.*(x(2)./V_G-K_L_G_Butene.*x(9)./V_L)+w_cat.*(k1.*x(8).^2./V_L.^2-k2.*x(8).*x(9)./V_L.^2-2.*k4.*x(9).^2./V_L.^2-k6.*x(9).*x(10)./V_L.^2-k8.*x(9).*x(11)./V_L.^2);V_R.*k_L_a_G_L_Light.*(x(3)./V_G-K_L_G_Hexene.*x(10)./V_L)+w_cat.*(k2.*x(8).*x(9)./V_L.^2-k3.*x(8).*x(10)./V_L.^2-k6.*x(9).*x(10)./V_L.^2-2.*k9.*x(10).^2./V_L.^2);V_R.*k_L_a_G_L_Heavy.*(x(4)./V_G-K_L_G_Octene.*x(11)./V_L)+w_cat.*(k3.*x(8).*x(10)./V_L.^2+k4.*x(9).^2./V_L.^2-k5.*x(8).*x(11)./V_L.^2-k8.*x(9).*x(11)./V_L.^2);V_R.*k_L_a_G_L_Heavy.*(x(5)./V_G-K_L_G_Decene.*x(12)./V_L)+w_cat.*(k5.*x(8).*x(11)./V_L.^2+k6.*x(9).*x(10)./V_L.^2-k7.*x(8).*x(12)./V_L.^2);V_R.*k_L_a_G_L_Heavy.*(x(6)./V_G-K_L_G_Dodecene.*x(13)./V_L)+w_cat.*(k7.*x(8).*x(12)./V_L.^2+k8.*x(9).*x(10)./V_L.^2+k9.*x(10).^2./V_L.^2);V_R.*k_L_a_G_L_Heavy.*(x(7)./V_G-K_L_G_Undecane.*x(14)./V_L)];

[t,x]=ode45(dNdt,[0,18000],[n_initial_ethylene;0;0;0;0;0;0;0;0;0;0;0;0;10]);

plot(t, x), grid on, xlabel('t'), title('Moles of Products')

%moles of products

Moles_C2_G=x(end, 1);

Moles_C4_G=x(end,2);

Moles_C6_G=x(end,3);

Moles_C8_G=x(end,4);

Moles_C10_G=x(end,5);

Moles_C12_G=x(end,6);

Moles_C11H24_G=x(end,7);

Moles_C2_L=x(end, 8);

Moles_C4_L=x(end,9);

Moles_C6_L=x(end,10);

Moles_C8_L=x(end,11);

Moles_C10_L=x(end,12);

Moles_C12_L=x(end,13);

Moles_C11H24_L=x(end,14);

2 个评论
显示无隐藏无

Alan Stevens 2024-4-3

Are you sure your equations are correct? You get a negative value for Moles_C11H24_G. Can you have a negative number of moles?

naiva saeedia 2024-4-3

Thankyou so much for pointing out this problem..I have changed the initial conditions and now it fixed.. Equations are correct but my initial condition for C11H24_G was wrong

请先登录，再进行评论。

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

Answer 1

Sam Chak 2024-4-3

编辑：Sam Chak 2024-4-3

在 MATLAB Online 中打开

Hi @naiva saeedia

The system identification task at hand seems to be quite challenging, as you only have two parameters to manipulate in order to match the experimental final values of the system states {

to

} at 18000 seconds. This becomes even more complex due to the constraints imposed by the high-order dynamics and fixed initial values.

Update: The parameters have been updated to k_L_a_G_L_Light=0.000438931054106292 and k_L_a_G_L_Heavy=0.189834285630455. However, despite these adjustments, some states (moles of products) still exhibit negative values.

Have you ever conducted a stability analysis on the mathematical model of the system? Are the experimental data measured at steady-state?

T_Ethylene=230+273.15;

T_ref=230+273.15;

R=8.314;

d_T=12.5e-2;

d_I=3.15e-2;

Q_G_in_ml_min=580; %mL/min

Q_G_in=Q_G_in_ml_min.*1.66667e-8; %m3/s

F_G_out_mmol_s=0.354;

F_G_out=F_G_out_mmol_s.*1e-3;

V_R=300e-6;

V_G=250e-6;

V_L=50e-6;

%Calculating D_I_M for all is

T=T_Ethylene;%K

M_B=28;%g/mol

M_1=M_B;

M_2=56;

M_3=84;

M_4=112;

M_5=140;

M_6=168;

M_7=156;

P=35.5292;%atm

v_C=15.9;%cm3

v_H=2.31;%cm3

D_I_M_1=(1e-3.*T.^1.75.*(1./M_1+1./M_B))./(P.*((2*v_C+4*v_H)^(1/3)+(2*v_C+4*v_H)^(1/3)));

D_I_M_2=(1e-3.*T.^1.75.*(1./M_2+1./M_B))./(P.*((4*v_C+8*v_H)^(1/3)+(2*v_C+4*v_H)^(1/3)));

D_I_M_3=(1e-3.*T.^1.75.*(1./M_3+1./M_B))./(P.*((6*v_C+12*v_H)^(1/3)+(2*v_C+4*v_H)^(1/3)));

D_I_M_4=(1e-3.*T.^1.75.*(1./M_4+1./M_B))./(P.*((8*v_C+16*v_H)^(1/3)+(2*v_C+4*v_H)^(1/3)));

D_I_M_5=(1e-3.*T.^1.75.*(1./M_5+1./M_B))./(P.*((10*v_C+20*v_H)^(1/3)+(2*v_C+4*v_H)^(1/3)));

D_I_M_6=(1e-3.*T.^1.75.*(1./M_6+1./M_B))./(P.*((12*v_C+24*v_H)^(1/3)+(2*v_C+4*v_H)^(1/3)));

D_I_M_7=(1e-3.*T.^1.75.*(1./M_7+1./M_B))./(P.*((11*v_C+12*v_H)^(1/3)+(2*v_C+4*v_H)^(1/3)));

D_o_w=0.3173;

N_cd=(4.*Q_G_in.^0.5.*d_T.^0.25)./d_I.^2;

N_I=89.01;

U_G=1;

k_L_a_o_w=3.35*(N_I/N_cd)^1.464*U_G;

Lambda_L=5.515e-4;

Delta_L=4.095;

Lambda_H=2.320e-6;

Delta_H=2.207;

%Calculating k_L_i

k_L_a_G_L_Light=0.000438931054106292;

k_L_a_G_L_Heavy=0.189834285630455;

k1_0=2.224e-4;

k2_0=1.533e-4;

k3_0=3.988e-5;

k4_0=1.914e-7;

k5_0=4.328e-5;

k6_0=2.506e-7;

k7_0=4.036e-5;

k8_0=1.062e-6;

k9_0=6.055e-7;

E1=109.53;

E2=15.229;

E3=7.881;

E4=44.454;

E5=9.438;

E6=8.426;

E7=10.911;

E8=12.537;

E9=7.127;

k1=k1_0.*exp(-E1.*(1/T_Ethylene-1/T_ref)./R);

k2=k2_0.*exp(-E2.*(1/T_Ethylene-1/T_ref)./R);

k3=k3_0.*exp(-E3.*(1/T_Ethylene-1/T_ref)./R);

k4=k4_0.*exp(-E4.*(1/T_Ethylene-1/T_ref)./R);

k5=k5_0.*exp(-E5.*(1/T_Ethylene-1/T_ref)./R);

k6=k6_0.*exp(-E6.*(1/T_Ethylene-1/T_ref)./R);

k7=k7_0.*exp(-E7.*(1/T_Ethylene-1/T_ref)./R);

k8=k8_0.*exp(-E8.*(1/T_Ethylene-1/T_ref)./R);

k9=k9_0.*exp(-E9.*(1/T_Ethylene-1/T_ref)./R);

w_cat=(0.3+0.25)*1e-3;

%Ethylene:

P_Ethylene=36e5;

C_G_in=P_Ethylene/(R*T_Ethylene);

n_initial_ethylene=C_G_in.*V_R;

n_initial_solvent=0.2348;

Q_G_in_C_G_in_C4=0;

Q_G_in_C_G_in_C6=0;

Q_G_in_C_G_in_C8=0;

Q_G_in_C_G_in_C10=0;

Q_G_in_C_G_in_C12=0;

Q_G_in_C_G_in_C11H24=0;

K_L_G_Ethylene=3.24;

K_L_G_Butene=2.23;

K_L_G_Hexene=1.72;

K_L_G_Octene=1.3;

K_L_G_Decene=0.985;

K_L_G_Dodecene=0.751;

K_L_G_Undecane=0.842;

dNdt=@(t,x) [Q_G_in.*C_G_in-F_G_out-V_R.*k_L_a_G_L_Light.*(x(1)./V_G-K_L_G_Ethylene.*x(8)./V_L);Q_G_in_C_G_in_C4-F_G_out-V_R.*k_L_a_G_L_Light.*(x(2)./V_G-K_L_G_Butene.*(x(9)./V_L));Q_G_in_C_G_in_C6-F_G_out-V_R.*k_L_a_G_L_Light.*(x(3)./V_G-K_L_G_Hexene.*(x(10)./V_L));Q_G_in_C_G_in_C8-F_G_out-V_R.*k_L_a_G_L_Heavy.*(x(4)./V_G-K_L_G_Octene.*(x(11)./V_L));Q_G_in_C_G_in_C10-F_G_out-V_R.*k_L_a_G_L_Heavy.*(x(5)./V_G-K_L_G_Decene.*x(12)./V_L);Q_G_in_C_G_in_C12-F_G_out-V_R.*k_L_a_G_L_Heavy.*(x(6)./V_G-K_L_G_Dodecene.*x(13)./V_L);Q_G_in_C_G_in_C11H24-F_G_out-V_R.*k_L_a_G_L_Heavy.*(x(7)./V_G-K_L_G_Undecane.*x(14)./V_L);V_R.*k_L_a_G_L_Light.*(x(1)./V_G-K_L_G_Ethylene.*x(8)./V_L)+w_cat.*(-2.*k1.*(x(8)./V_L).^2-k2.*(x(8)./V_L).*(x(9)./V_L)-k3.*(x(8)./V_L).*(x(10)./V_L)-k5.*(x(8)./V_L).*(x(11)./V_L)-k7.*x(8).*x(12)./(V_L.^2));V_R.*k_L_a_G_L_Light.*(x(2)./V_G-K_L_G_Butene.*x(9)./V_L)+w_cat.*(k1.*x(8).^2./V_L.^2-k2.*x(8).*x(9)./V_L.^2-2.*k4.*x(9).^2./V_L.^2-k6.*x(9).*x(10)./V_L.^2-k8.*x(9).*x(11)./V_L.^2);V_R.*k_L_a_G_L_Light.*(x(3)./V_G-K_L_G_Hexene.*x(10)./V_L)+w_cat.*(k2.*x(8).*x(9)./V_L.^2-k3.*x(8).*x(10)./V_L.^2-k6.*x(9).*x(10)./V_L.^2-2.*k9.*x(10).^2./V_L.^2);V_R.*k_L_a_G_L_Heavy.*(x(4)./V_G-K_L_G_Octene.*x(11)./V_L)+w_cat.*(k3.*x(8).*x(10)./V_L.^2+k4.*x(9).^2./V_L.^2-k5.*x(8).*x(11)./V_L.^2-k8.*x(9).*x(11)./V_L.^2);V_R.*k_L_a_G_L_Heavy.*(x(5)./V_G-K_L_G_Decene.*x(12)./V_L)+w_cat.*(k5.*x(8).*x(11)./V_L.^2+k6.*x(9).*x(10)./V_L.^2-k7.*x(8).*x(12)./V_L.^2);V_R.*k_L_a_G_L_Heavy.*(x(6)./V_G-K_L_G_Dodecene.*x(13)./V_L)+w_cat.*(k7.*x(8).*x(12)./V_L.^2+k8.*x(9).*x(10)./V_L.^2+k9.*x(10).^2./V_L.^2);V_R.*k_L_a_G_L_Heavy.*(x(7)./V_G-K_L_G_Undecane.*x(14)./V_L)];

[t,x]=ode45(dNdt,[0,18000],[n_initial_ethylene;0;0;0;0;0;0;0;0;0;0;0;0;10]);

%% plot results

plot(t, x), grid on, xlabel('t'), title('Moles of Products')

%% moles of products (simulated final values)

Moles_C2_G = x(end,1)

Moles_C2_G = 15.2246

Moles_C4_G = x(end,2)

Moles_C4_G = -0.5886

Moles_C6_G = x(end,3)

Moles_C6_G = -0.4970

Moles_C8_G = x(end,4)

Moles_C8_G = 0.0805

Moles_C10_G = x(end,5)

Moles_C10_G = 0.0352

Moles_C12_G = x(end,6)

Moles_C12_G = 0.5796

Moles_C11H24_G = x(end,7)

Moles_C11H24_G = 2.9316

Moles_C2_L = x(end,8)

Moles_C2_L = 0.0064

Moles_C4_L = x(end,9)

Moles_C4_L = 0.0075

Moles_C6_L = x(end,10)

Moles_C6_L = 0.0204

Moles_C8_L = x(end,11)

Moles_C8_L = 0.0126

Moles_C10_L = x(end,12)

Moles_C10_L = 0.0075

Moles_C12_L = x(end,13)

Moles_C12_L = 0.1549

%% Measured experimental final values (data)

x2exp = 0.07891;

x3exp = 0.14349;

x4exp = 0.063348;

x5exp = 0.02686;

x6exp = 0.01351;

%% errors

x2err = x(end,2) - x2exp

x2err = -0.6675

x3err = x(end,3) - x3exp

x3err = -0.6405

x4err = x(end,4) - x4exp

x4err = 0.0172

x5err = x(end,5) - x5exp

x5err = 0.0083

x6err = x(end,6) - x6exp

x6err = 0.5660

29 个评论
显示 27更早的评论隐藏 27更早的评论

naiva saeedia 2024-4-5

编辑：naiva saeedia 2024-4-5

在 MATLAB Online 中打开

@Torsten

Did I do it correctly? and what does that message mean?

clear
clc
f=@(x) [0.2278011427*x(7) - 0.9590428109*x(14);-0.0003540000000 - 0.004800000000*x(2) + 0.05352000000*x(9);-0.0003540000000 - 0.004800000000*x(3) + 0.04128000000*x(10);-0.0003540000000 - 0.2278011427*x(4) + 1.480707428*x(11);-0.0003540000000 - 0.2278011427*x(5) + 1.121920628*x(12);-0.0003540000000 - 0.2278011427*x(6) + 0.8553932909*x(13);-0.0003540000000 - 0.2278011427*x(7) + 0.9590428109*x(14);0.007965028561 - 0.004800000000*x(1) + 0.07776000000*x(8);0.2278011427*x(5) - 1.121920628*x(12) + 9.521600000*x(8)*x(11) + 0.05513200000*x(9)*x(10) - 8.879200000*x(8)*x(12);0.2278011427*x(6) - 0.8553932909*x(13) + 8.879200000*x(8)*x(12) + 0.2336400000*x(9)*x(10) + 0.1332100000*x(10)^2;0.004800000000*x(3) - 0.04128000000*x(10) + 33.72600000*x(8)*x(9) - 8.773600000*x(8)*x(10) - 0.05513200000*x(9)*x(10) - 0.2664200000*x(10)^2;0.2278011427*x(4) - 1.480707428*x(11) + 8.773600000*x(8)*x(10) + 0.04210800000*x(9)^2 - 9.521600000*x(8)*x(11) - 0.2336400000*x(9)*x(11);0.004800000000*x(1) - 0.07776000000*x(8) - 97.85600000*x(8)^2 - 33.72600000*x(8)*x(9) - 8.773600000*x(8)*x(10) - 9.521600000*x(8)*x(11) - 8.879200000*x(8)*x(12);0.004800000000*x(2) - 0.05352000000*x(9) + 48.92800000*x(8)^2 - 33.72600000*x(8)*x(9) - 0.08421600000*x(9)^2 - 0.05513200000*x(9)*x(10) - 0.2336400000*x(9)*x(11)];
lb=[0,0,0,0,0,0,0,0,0,0,0,0,0,0];
ub=[1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000];
x0=[0,0,0,0,0,0,0,0,0,0,0,0,0,0];
x=lsqnonlin(f,x0,lb,ub)
Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.
x = 1x14
    0.8317    0.0000    0.0000    0.0020    0.0005    0.0000    0.0000    0.0058    0.0038    0.0047    0.0005    0.0003    0.0002    0.0002
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>

naiva saeedia 2024-4-6

在 MATLAB Online 中打开

@Torsten

Hi again. Thx for your answer. I don't think I can possibly alter my model and exculde the components which don't have initial conditions but I can run the model with approxiamte values of k_L_a_G_L_Light and k_L_a_G_L_Heavy to find initial conditions for the rest of the concentrations. As you can see now I have all the concentrations at t=t_end(These concentraions(except x(2) to x(6)..I have these based on experimental datas) was computed at k_L_a_G_L_Light=0.055 and k_L_a_G_L_Heavy=0.029)

x(1)=0.2244

x(2)=0.07891

x(3)=0.14349

x(4)=0.063348

x(5)=0.02686

x(6)=0.01351

x(7)=2.9313

x(8)=0.006401

x(9)=0.007486

x(10)=0.020358

x(11)=0.012628

x(12)=0.007462

x(13)=0.57007

x(14)=0.69673

Now with these conditions I used Isqnonlin and I assumed that k_L_a_G_L_Light and k_L_a_G_L_Heavy are x(15) and x(16) in the equations. As you can see I have found k_L_a_G_L_Light and k_L_a_G_L_Heavy values. Is this correct in your opinion?

clc
clear
close all
T_Ethylene=230+273.15;
T_ref=230+273.15;
R=8.314;
d_T=12.5e-2;
d_I=3.15e-2;
Q_G_in_ml_min=580; %mL/min
Q_G_in=Q_G_in_ml_min.*1.66667e-8; %m3/s
F_G_out_mmol_s=0.354;
F_G_out=F_G_out_mmol_s.*1e-3;
V_R=300e-6;
V_G=250e-6;
V_L=50e-6;
%Calculating D_I_M for all is
T=T_Ethylene;%K
M_B=28;%g/mol
M_1=M_B;
M_2=56;
M_3=84;
M_4=112;
M_5=140;
M_6=168;
M_7=156;
P=35.5292;%atm
v_C=15.9;%cm3
v_H=2.31;%cm3
D_I_M_1=(1e-3.*T.^1.75.*(1./M_1+1./M_B))./(P.*((2*v_C+4*v_H)^(1/3)+(2*v_C+4*v_H)^(1/3)));
D_I_M_2=(1e-3.*T.^1.75.*(1./M_2+1./M_B))./(P.*((4*v_C+8*v_H)^(1/3)+(2*v_C+4*v_H)^(1/3)));
D_I_M_3=(1e-3.*T.^1.75.*(1./M_3+1./M_B))./(P.*((6*v_C+12*v_H)^(1/3)+(2*v_C+4*v_H)^(1/3)));
D_I_M_4=(1e-3.*T.^1.75.*(1./M_4+1./M_B))./(P.*((8*v_C+16*v_H)^(1/3)+(2*v_C+4*v_H)^(1/3)));
D_I_M_5=(1e-3.*T.^1.75.*(1./M_5+1./M_B))./(P.*((10*v_C+20*v_H)^(1/3)+(2*v_C+4*v_H)^(1/3)));
D_I_M_6=(1e-3.*T.^1.75.*(1./M_6+1./M_B))./(P.*((12*v_C+24*v_H)^(1/3)+(2*v_C+4*v_H)^(1/3)));
D_I_M_7=(1e-3.*T.^1.75.*(1./M_7+1./M_B))./(P.*((11*v_C+12*v_H)^(1/3)+(2*v_C+4*v_H)^(1/3)));
D_o_w=0.3173;
N_cd=(4.*Q_G_in.^0.5.*d_T.^0.25)./d_I.^2;
N_I=89.01;
U_G=1;
k_L_a_o_w=3.35*(N_I/N_cd)^1.464*U_G;
Lambda_L=5.515e-4;
Delta_L=4.095;
Lambda_H=2.320e-6;
Delta_H=2.207;
%Calculating k_L_i
k1_0=2.224e-4;
k2_0=1.533e-4;
k3_0=3.988e-5;
k4_0=1.914e-7;
k5_0=4.328e-5;
k6_0=2.506e-7;
k7_0=4.036e-5;
k8_0=1.062e-6;
k9_0=6.055e-7;
E1=109.53;
E2=15.229;
E3=7.881;
E4=44.454;
E5=9.438;
E6=8.426;
E7=10.911;
E8=12.537;
E9=7.127;
k1=k1_0.*exp(-E1.*(1/T_Ethylene-1/T_ref)./R);
k2=k2_0.*exp(-E2.*(1/T_Ethylene-1/T_ref)./R);
k3=k3_0.*exp(-E3.*(1/T_Ethylene-1/T_ref)./R);
k4=k4_0.*exp(-E4.*(1/T_Ethylene-1/T_ref)./R);
k5=k5_0.*exp(-E5.*(1/T_Ethylene-1/T_ref)./R);
k6=k6_0.*exp(-E6.*(1/T_Ethylene-1/T_ref)./R);
k7=k7_0.*exp(-E7.*(1/T_Ethylene-1/T_ref)./R);
k8=k8_0.*exp(-E8.*(1/T_Ethylene-1/T_ref)./R);
k9=k9_0.*exp(-E9.*(1/T_Ethylene-1/T_ref)./R);
w_cat=(0.3+0.25)*1e-3;
%Ethylene:
P_Ethylene=36e5;
C_G_in=P_Ethylene/(R*T_Ethylene);
n_initial_ethylene=C_G_in.*V_R;
n_initial_solvent=0.2348;
Q_G_in_C_G_in_C4=0;
Q_G_in_C_G_in_C6=0;
Q_G_in_C_G_in_C8=0;
Q_G_in_C_G_in_C10=0;
Q_G_in_C_G_in_C12=0;
Q_G_in_C_G_in_C11H24=0;
K_L_G_Ethylene=3.24;
K_L_G_Butene=2.23;
K_L_G_Hexene=1.72;
K_L_G_Octene=1.3;
K_L_G_Decene=0.985;
K_L_G_Dodecene=0.751;
K_L_G_Undecane=0.842;
fun=@(x) [Q_G_in.*C_G_in-F_G_out-V_R.*x(15).*(x(1)./V_G-K_L_G_Ethylene.*x(8)./V_L);Q_G_in_C_G_in_C4-F_G_out-V_R.*x(15).*(x(2)./V_G-K_L_G_Butene.*(x(9)./V_L));Q_G_in_C_G_in_C6-F_G_out-V_R.*x(15).*(x(3)./V_G-K_L_G_Hexene.*(x(10)./V_L));Q_G_in_C_G_in_C8-F_G_out-V_R.*x(16).*(x(4)./V_G-K_L_G_Octene.*(x(11)./V_L));Q_G_in_C_G_in_C10-F_G_out-V_R.*x(16).*(x(5)./V_G-K_L_G_Decene.*x(12)./V_L);Q_G_in_C_G_in_C12-F_G_out-V_R.*x(16).*(x(6)./V_G-K_L_G_Dodecene.*x(13)./V_L);Q_G_in_C_G_in_C11H24-F_G_out-V_R.*x(16).*(x(7)./V_G-K_L_G_Undecane.*x(14)./V_L);V_R.*x(15).*(x(1)./V_G-K_L_G_Ethylene.*x(8)./V_L)+w_cat.*(-2.*k1.*(x(8)./V_L).^2-k2.*(x(8)./V_L).*(x(9)./V_L)-k3.*(x(8)./V_L).*(x(10)./V_L)-k5.*(x(8)./V_L).*(x(11)./V_L)-k7.*x(8).*x(12)./(V_L.^2));V_R.*x(15).*(x(2)./V_G-K_L_G_Butene.*x(9)./V_L)+w_cat.*(k1.*x(8).^2./V_L.^2-k2.*x(8).*x(9)./V_L.^2-2.*k4.*x(9).^2./V_L.^2-k6.*x(9).*x(10)./V_L.^2-k8.*x(9).*x(11)./V_L.^2);V_R.*x(15).*(x(3)./V_G-K_L_G_Hexene.*x(10)./V_L)+w_cat.*(k2.*x(8).*x(9)./V_L.^2-k3.*x(8).*x(10)./V_L.^2-k6.*x(9).*x(10)./V_L.^2-2.*k9.*x(10).^2./V_L.^2);V_R.*x(16).*(x(4)./V_G-K_L_G_Octene.*x(11)./V_L)+w_cat.*(k3.*x(8).*x(10)./V_L.^2+k4.*x(9).^2./V_L.^2-k5.*x(8).*x(11)./V_L.^2-k8.*x(9).*x(11)./V_L.^2);V_R.*x(16).*(x(5)./V_G-K_L_G_Decene.*x(12)./V_L)+w_cat.*(k5.*x(8).*x(11)./V_L.^2+k6.*x(9).*x(10)./V_L.^2-k7.*x(8).*x(12)./V_L.^2);V_R.*x(16).*(x(6)./V_G-K_L_G_Dodecene.*x(13)./V_L)+w_cat.*(k7.*x(8).*x(12)./V_L.^2+k8.*x(9).*x(10)./V_L.^2+k9.*x(10).^2./V_L.^2);V_R.*x(16).*(x(7)./V_G-K_L_G_Undecane.*x(14)./V_L)];
x0=[0.2244, 0.07891, 0.14349, 0.063348, 0.02686, 0.01351, 2.9313, 0.006401, 0.007486, 0.020358, 0.012628, 0.007462, 0.57007, 0.69673, 0.055, 0.029 ];
lb=[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0];
ub=[1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000,1000];
x=lsqnonlin(fun,x0,lb,ub)
Warning: Trust-region-reflective algorithm requires at least as many equations as variables; using Levenberg-Marquardt algorithm instead.
Local minimum found.

Optimization completed because the size of the gradient is less than
1e-4 times the value of the function tolerance.
x = 1x16
    0.2276    0.0783    0.1671    0.0490         0    0.1525    2.9307    0.0065    0.0075    0.0201    0.0118    0.0057    0.0482    0.6992    0.0541    0.0113
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>

Sam Chak 2024-4-6

Hi @naiva saeedia

Have you confirmed that the super long ODEs defined in a single-line 'fun' are identical to the ODEs presented in Woo's paper?

fun=@(x) [Q_G_in.*C_G_in-F_G_out-V_R.*x(15).*(x(1)./V_G-K_L_G_Ethylene.*x(8)./V_L);Q_G_in_C_G_in_C4-F_G_out-V_R.*x(15).*(x(2)./V_G-K_L_G_Butene.*(x(9)./V_L));Q_G_in_C_G_in_C6-F_G_out-V_R.*x(15).*(x(3)./V_G-K_L_G_Hexene.*(x(10)./V_L));Q_G_in_C_G_in_C8-F_G_out-V_R.*x(16).*(x(4)./V_G-K_L_G_Octene.*(x(11)./V_L));Q_G_in_C_G_in_C10-F_G_out-V_R.*x(16).*(x(5)./V_G-K_L_G_Decene.*x(12)./V_L);Q_G_in_C_G_in_C12-F_G_out-V_R.*x(16).*(x(6)./V_G-K_L_G_Dodecene.*x(13)./V_L);Q_G_in_C_G_in_C11H24-F_G_out-V_R.*x(16).*(x(7)./V_G-K_L_G_Undecane.*x(14)./V_L);V_R.*x(15).*(x(1)./V_G-K_L_G_Ethylene.*x(8)./V_L)+w_cat.*(-2.*k1.*(x(8)./V_L).^2-k2.*(x(8)./V_L).*(x(9)./V_L)-k3.*(x(8)./V_L).*(x(10)./V_L)-k5.*(x(8)./V_L).*(x(11)./V_L)-k7.*x(8).*x(12)./(V_L.^2));V_R.*x(15).*(x(2)./V_G-K_L_G_Butene.*x(9)./V_L)+w_cat.*(k1.*x(8).^2./V_L.^2-k2.*x(8).*x(9)./V_L.^2-2.*k4.*x(9).^2./V_L.^2-k6.*x(9).*x(10)./V_L.^2-k8.*x(9).*x(11)./V_L.^2);V_R.*x(15).*(x(3)./V_G-K_L_G_Hexene.*x(10)./V_L)+w_cat.*(k2.*x(8).*x(9)./V_L.^2-k3.*x(8).*x(10)./V_L.^2-k6.*x(9).*x(10)./V_L.^2-2.*k9.*x(10).^2./V_L.^2);V_R.*x(16).*(x(4)./V_G-K_L_G_Octene.*x(11)./V_L)+w_cat.*(k3.*x(8).*x(10)./V_L.^2+k4.*x(9).^2./V_L.^2-k5.*x(8).*x(11)./V_L.^2-k8.*x(9).*x(11)./V_L.^2);V_R.*x(16).*(x(5)./V_G-K_L_G_Decene.*x(12)./V_L)+w_cat.*(k5.*x(8).*x(11)./V_L.^2+k6.*x(9).*x(10)./V_L.^2-k7.*x(8).*x(12)./V_L.^2);V_R.*x(16).*(x(6)./V_G-K_L_G_Dodecene.*x(13)./V_L)+w_cat.*(k7.*x(8).*x(12)./V_L.^2+k8.*x(9).*x(10)./V_L.^2+k9.*x(10).^2./V_L.^2);V_R.*x(16).*(x(7)./V_G-K_L_G_Undecane.*x(14)./V_L)];

Additionally, I noticed that there are fourteen initial values provided, but only five final values are measured. Did you utilize im2graph on Figure 2 to estimate both the fourteen initial values and the five final values?

Sam Chak 2024-4-11

编辑：Sam Chak 2024-4-11

在 MATLAB Online 中打开

Hi @naiva saeedia

In theory, is the given chemical system expected to be thermodynamically stable? I am not an expert in this field, so I cannot provide a definitive answer.

However, from a mathematical perspective, we can assess whether the chemical reactions described by the differential equations tend towards equilibrium. By examining the eigenvalues, we can determine if the equilibrium is stable or unstable. If we are unable to find an equilibrium by solving the algebraic equations, it can generally be inferred that the system is unstable.

As for how to write these equations in MATLAB to improve the readability, you can consider writing as follows. I only demo for the first seven ODEs.

Update: The code snippet has been revised to include the remaining seven ODEs. Please note that there are some constants in the code that are not being utilized. This is how the code originally appeared.

%% Parameters to be determined
kL  = param(1);         % k_L_a_G_L_Light, 1st parameter 
kH  = param(2);         % k_L_a_G_L_Heavy, 2nd parameter 
%% Constants
Q0  = 580;              % Q_G_in_ml_min
Q1  = Q0*1.66667e-8;    % Q_G_in (m3/s)
Q2  = 0;                % Q_G_in_C_G_in_C4
Q3  = 0;                % Q_G_in_C_G_in_C6
Q4  = 0;                % Q_G_in_C_G_in_C8
Q5  = 0;                % Q_G_in_C_G_in_C10
Q6  = 0;                % Q_G_in_C_G_in_C12
Q7  = 0;                % Q_G_in_C_G_in_C11H24
P1  = 36e5;             % P_Ethylene
T1  = 230 + 273.15;     % T_Ethylene
T2  = 230 + 273.15;     % T_ref
R   = 8.314;            % Absolute Gas constant in the Universe
C1  = P1/(R*T1);        % C_G_in
F0  = 0.354;            % F_G_out_mmol_s
F1  = F0*1e-3;          % F_G_out
VR  = 300e-6;           % V_R
VG  = 250e-6;           % V_G
VL  = 50e-6;            % V_L
K1  = 3.24;             % K_L_G_Ethylene
K2  = 2.23;             % K_L_G_Butene
K3  = 1.72;             % K_L_G_Hexene
K4  = 1.3;              % K_L_G_Octene
K5  = 0.985;            % K_L_G_Decene
K6  = 0.751;            % K_L_G_Dodecene
K7  = 0.842;            % K_L_G_Undecane
wc  = (0.3 + 0.25)*1e-3;    % w_cat
k10 = 2.224e-4;
k20 = 1.533e-4;
k30 = 3.988e-5;
k40 = 1.914e-7;
k50 = 4.328e-5;
k60 = 2.506e-7;
k70 = 4.036e-5;
k80 = 1.062e-6;
k90 = 6.055e-7;
E1  = 109.53;
E2  =  15.229;
E3  =   7.881;
E4  =  44.454;
E5  =   9.438;
E6  =   8.426;
E7  =  10.911;
E8  =  12.537;
E9  =   7.127;
k1  = k10*exp(- E1*(1/T1 - 1/T2)/R);
k2  = k20*exp(- E2*(1/T1 - 1/T2)/R);
k3  = k30*exp(- E3*(1/T1 - 1/T2)/R);
k4  = k40*exp(- E4*(1/T1 - 1/T2)/R);
k5  = k50*exp(- E5*(1/T1 - 1/T2)/R);
k6  = k60*exp(- E6*(1/T1 - 1/T2)/R);
k7  = k70*exp(- E7*(1/T1 - 1/T2)/R);
k8  = k80*exp(- E8*(1/T1 - 1/T2)/R);
k9  = k90*exp(- E9*(1/T1 - 1/T2)/R);
%% ODEs
dN1 = Q1*C1 - F1 - VR*kL*(x(1)/VG - K1*x(8) /VL);
dN2 = Q2    - F1 - VR*kL*(x(2)/VG - K2*x(9) /VL);
dN3 = Q3    - F1 - VR*kL*(x(3)/VG - K3*x(10)/VL);
dN4 = Q4    - F1 - VR*kH*(x(4)/VG - K4*x(11)/VL);
dN5 = Q5    - F1 - VR*kH*(x(5)/VG - K5*x(12)/VL);
dN6 = Q6    - F1 - VR*kH*(x(6)/VG - K6*x(13)/VL);
dN7 = Q7    - F1 - VR*kH*(x(7)/VG - K7*x(14)/VL);
dN8 = VR*kL*(x(1)/VG - K1*x(8) /VL) + wc*(-2*k1*x(8)^2    /VL^2 - k2*x(8)*x(9) /VL^2 -   k3*x(8)*x(10)/VL^2 -   k5*x(8)*x(11)/VL^2 - k7*x(8)*x(12)/VL^2);
dN9 = VR*kL*(x(2)/VG - K2*x(9) /VL) + wc*(   k1*x(8)^2    /VL^2 - k2*x(8)*x(9) /VL^2 - 2*k4*x(9)^2    /VL^2 -   k6*x(9)*x(10)/VL^2 - k8*x(9)*x(11)/VL^2);
dN10= VR*kL*(x(3)/VG - K3*x(10)/VL) + wc*(   k2*x(8)*x(9) /VL^2 - k3*x(8)*x(10)/VL^2 -   k6*x(9)*x(10)/VL^2 - 2*k9*x(10)^2   /VL^2);
dN11= VR*kH*(x(4)/VG - K4*x(11)/VL) + wc*(   k3*x(8)*x(10)/VL^2 + k4*x(9)^2    /VL^2 -   k5*x(8)*x(11)/VL^2 -   k8*x(9)*x(11)/VL^2);
dN12= VR*kH*(x(5)/VG - K5*x(12)/VL) + wc*(   k5*x(8)*x(11)/VL^2 + k6*x(9)*x(10)/VL^2 -   k7*x(8)*x(12)/VL^2);
dN13= VR*kH*(x(6)/VG - K6*x(13)/VL) + wc*(   k7*x(8)*x(12)/VL^2 + k8*x(9)*x(10)/VL^2 +   k9*x(10)^2   /VL^2);
dN14= VR*kH*(x(7)/VG - K7*x(14)/VL);
%% Group ODEs as a column vector
dN  = [dN1;
       dN2;
       dN3;
       dN4;
       dN5;
       dN6;
       dN7;
       dN8;
       dN9;
       dN10;
       dN11;
       dN12;
       dN13;
       dN14];

Sam Chak 2024-4-11

在 MATLAB Online 中打开

Hi @naiva saeedia

I haven't conducted exhaustive testing on the system, but I have observed that when the "external mass flows" 'Q1' and 'F1' are set to zero, and both parameters

and

are equal to or greater than zero, the system tends to converge towards equilibrium points. The specific equilibrium points depend on the actual values assigned to

and

.

It is possible to perform an exhaustive search of different combinations of

and

within a specified range to determine the desired equilibrium. However, it's important to note that this approach only works if the mass transfer system is stable. Introduction of non-zero values for 'Q1' and 'F1' seems to destabilize the system.

tspan = [0 1e4];

x0 = [ones(7, 1); zeros(7, 1)];

[t, x] = ode45(@massTransfer, tspan, x0);

plot(t, x), grid on

xlabel('t'), ylabel('\bf{x}\rm(t)'), title('System Response')

function dNdt = massTransfer(t, x)

%% Parameters to be determined

kL = 1; % k_L_a_G_L_Light, 1st parameter

kH = 1; % k_L_a_G_L_Heavy, 2nd parameter

%% Constants

Q0 = 580; % Q_G_in_ml_min

% Q1 = Q0*1.66667e-8; % Q_G_in (m3/s)

Q1 = 0; % No "enthalpy"

Q2 = 0; % Q_G_in_C_G_in_C4

Q3 = 0; % Q_G_in_C_G_in_C6

Q4 = 0; % Q_G_in_C_G_in_C8

Q5 = 0; % Q_G_in_C_G_in_C10

Q6 = 0; % Q_G_in_C_G_in_C12

Q7 = 0; % Q_G_in_C_G_in_C11H24

P1 = 36e5; % P_Ethylene

T1 = 230 + 273.15; % T_Ethylene

T2 = 230 + 273.15; % T_ref

R = 8.314; % Absolute Gas constant in the Universe

C1 = P1/(R*T1); % C_G_in

F0 = 0.354; % F_G_out_mmol_s

% F1 = F0*1e-3; % F_G_out

F1 = 0; % No "enthalpy"

VR = 300e-6; % V_R

VG = 250e-6; % V_G

VL = 50e-6; % V_L

K1 = 3.24; % K_L_G_Ethylene

K2 = 2.23; % K_L_G_Butene

K3 = 1.72; % K_L_G_Hexene

K4 = 1.3; % K_L_G_Octene

K5 = 0.985; % K_L_G_Decene

K6 = 0.751; % K_L_G_Dodecene

K7 = 0.842; % K_L_G_Undecane

wc = (0.3 + 0.25)*1e-3; % w_cat

k10 = 2.224e-4;

k20 = 1.533e-4;

k30 = 3.988e-5;

k40 = 1.914e-7;

k50 = 4.328e-5;

k60 = 2.506e-7;

k70 = 4.036e-5;

k80 = 1.062e-6;

k90 = 6.055e-7;

E1 = 109.53;

E2 = 15.229;

E3 = 7.881;

E4 = 44.454;

E5 = 9.438;

E6 = 8.426;

E7 = 10.911;

E8 = 12.537;

E9 = 7.127;

k1 = k10*exp(- E1*(1/T1 - 1/T2)/R);

k2 = k20*exp(- E2*(1/T1 - 1/T2)/R);

k3 = k30*exp(- E3*(1/T1 - 1/T2)/R);

k4 = k40*exp(- E4*(1/T1 - 1/T2)/R);

k5 = k50*exp(- E5*(1/T1 - 1/T2)/R);

k6 = k60*exp(- E6*(1/T1 - 1/T2)/R);

k7 = k70*exp(- E7*(1/T1 - 1/T2)/R);

k8 = k80*exp(- E8*(1/T1 - 1/T2)/R);

k9 = k90*exp(- E9*(1/T1 - 1/T2)/R);

%% ODEs

dN1 = Q1*C1 - F1 - VR*kL*(x(1)/VG - K1*x(8) /VL);

dN2 = Q2 - F1 - VR*kL*(x(2)/VG - K2*x(9) /VL);

dN3 = Q3 - F1 - VR*kL*(x(3)/VG - K3*x(10)/VL);

dN4 = Q4 - F1 - VR*kH*(x(4)/VG - K4*x(11)/VL);

dN5 = Q5 - F1 - VR*kH*(x(5)/VG - K5*x(12)/VL);

dN6 = Q6 - F1 - VR*kH*(x(6)/VG - K6*x(13)/VL);

dN7 = Q7 - F1 - VR*kH*(x(7)/VG - K7*x(14)/VL);

dN8 = VR*kL*(x(1)/VG - K1*x(8) /VL) + wc*(-2*k1*x(8)^2 /VL^2 - k2*x(8)*x(9) /VL^2 - k3*x(8)*x(10)/VL^2 - k5*x(8)*x(11)/VL^2 - k7*x(8)*x(12)/VL^2);

dN9 = VR*kL*(x(2)/VG - K2*x(9) /VL) + wc*( k1*x(8)^2 /VL^2 - k2*x(8)*x(9) /VL^2 - 2*k4*x(9)^2 /VL^2 - k6*x(9)*x(10)/VL^2 - k8*x(9)*x(11)/VL^2);

dN10= VR*kL*(x(3)/VG - K3*x(10)/VL) + wc*( k2*x(8)*x(9) /VL^2 - k3*x(8)*x(10)/VL^2 - k6*x(9)*x(10)/VL^2 - 2*k9*x(10)^2 /VL^2);

dN11= VR*kH*(x(4)/VG - K4*x(11)/VL) + wc*( k3*x(8)*x(10)/VL^2 + k4*x(9)^2 /VL^2 - k5*x(8)*x(11)/VL^2 - k8*x(9)*x(11)/VL^2);

dN12= VR*kH*(x(5)/VG - K5*x(12)/VL) + wc*( k5*x(8)*x(11)/VL^2 + k6*x(9)*x(10)/VL^2 - k7*x(8)*x(12)/VL^2);

dN13= VR*kH*(x(6)/VG - K6*x(13)/VL) + wc*( k7*x(8)*x(12)/VL^2 + k8*x(9)*x(10)/VL^2 + k9*x(10)^2 /VL^2);

dN14= VR*kH*(x(7)/VG - K7*x(14)/VL);

%% Group ODEs as a column vector

dNdt= [dN1;

dN2;

dN3;

dN4;

dN5;

dN6;

dN7;

dN8;

dN9;

dN10;

dN11;

dN12;

dN13;

dN14];

end

naiva saeedia 2024-4-11

在 MATLAB Online 中打开

Hi @Sam Chak

Q1 is the volumetric flow rate of the ethylene that enters the reactor when you set its value to zero that means we don't enter any ethylene to the reactor and because it's a semi-batch reactor it doesn't make sense. F1 is the molar flow of the ethylene that exits the reactor if we set its value to zero that means the reactor is batch not semi-batch but what you did is giving me some insights maybe reactor is a batch reactor after all and the paper is wrong about the reactor type..I have a question. Can we solve an optimization problem with this batch reactor to find optimal values for kL and kH? I want these constraints for the optimization problem:

kL>0

kH>0

Mass_C2_L>0

Mass_C4_L>0

Mass_C6_L>0

Mass_C8_L>0

Mass_C10_L>0

Mass_C12_L>0

Mass_C11H24_L>0

Mass_C2_G>0

Mass_C4_G>0

Mass_C6_G>0

Mass_C8_G>0

Mass_C10_G>0

Mass_C12_G>0

Mass_C11H24_G>0

Also:

Mass_C4_G-4.41901<0.01

Mass_C6_G-10.4754<0.01

Mass_C8_G-7.09507<0.01

Mass_C10_G-4.13732<0.01

Mass_C12_G-2.27113<0.01

453.15=<T1<=503.15

I also added these parameters to your code. Can you please help me to define this optimization problem in Matlab?

tspan   = [0 1e4];
x0      = [ones(7, 1); zeros(7, 1)];
[t, x]  = ode45(@massTransfer, tspan, x0);
plot(t, x), grid on
xlabel('t'), ylabel('\bf{x}\rm(t)'), title('System Response')
Moles_C2_G=x(end, 1);
Moles_C4_G=x(end,2);
Moles_C6_G=x(end,3);
Moles_C8_G=x(end,4);
Moles_C10_G=x(end,5);
Moles_C12_G=x(end,6);
Moles_C11H24_G=x(end,7);
Moles_C2_L=x(end, 8);
Moles_C4_L=x(end,9);
Moles_C6_L=x(end,10);
Moles_C8_L=x(end,11);
Moles_C10_L=x(end,12);
Moles_C12_L=x(end,13);
Moles_C11H24_L=x(end,14);
%Grams of products
Mass_C2_G=Moles_C2_G.*28;
Mass_C4_G=Moles_C4_G.*56;
Mass_C6_G=Moles_C6_G.*73;
Mass_C8_G=Moles_C8_G.*112;
Mass_C10_G=Moles_C10_G.*154;
Mass_C12_G=Moles_C12_G.*168;
Mass_C11H24_G=Moles_C11H24_G.*156;
Mass_C2_L=Moles_C2_L.*28;
Mass_C4_L=Moles_C4_L.*56;
Mass_C6_L=Moles_C6_L.*73;
Mass_C8_L=Moles_C8_L.*112;
Mass_C10_L=Moles_C10_L.*154;
Mass_C12_L=Moles_C12_L.*168;
Mass_C11H24_L=Moles_C11H24_L.*156;

Sam Chak 2024-4-11

在 MATLAB Online 中打开

Hi @naiva saeedia

To perform optimization, we need a real-valued cost function that varies with the parameters

and

. If this function is convex and has a single minimum point within the specified search region, MATLAB Optimization Toolbox offers various tools that can be utilized.

The optimal values for

and

can be loosely interpreted as the most favorable combination of values to achieve a quantifiable objective. The specific objective in this case is not known to me, but you likely have a better understanding.

For instance, in your Mass Transfer system, we could inquire, "What are the best values for

and

that allow the selected masses to satisfy the given constraints (

) in the fastest manner?" If the cost function J in the image below represents the normalized time taken to satisfy the constraint, then we can conclude that

and

are the best values.

[X, Y] = meshgrid(0:100);

Z = 0.5*(((X - 50).^2 + (Y - 50).^2)/5000) + 0.5;

mesh(Z)

xlabel({'$k_{\mathrm{Light}}$'}, 'Interpreter', 'LaTeX', 'fontsize', 16)

ylabel({'$k_{\mathrm{Heavy}}$'}, 'Interpreter', 'LaTeX', 'fontsize', 16)

zlabel({'Cost function, $\mathcal{J}$'}, 'Interpreter', 'LaTeX', 'fontsize', 16)

Sam Chak 2024-4-13

Hi @naiva saeedia,

Have you been able to identify the optimal combinations of

and

that result in the desired steady-state masses precisely at 18000 seconds? It's important to remember that the optimization approach is only meaningful if the mass transfer system is in equilibrium (stable).

My point is that even if there is a specific set of

and

values that meet the objective requirements and yield the desired masses at 18000 seconds, the outcome may be meaningless if the masses continue to deviate significantly over time for

seconds. The persistent deviation serves as an indication that the system may be unstable.

naiva saeedia 2024-4-13

编辑：naiva saeedia 2024-4-13

@Sam Chak

Hi. yes I have found the values for kL and kH that give me the minimum error in mass of the products but the thing is as you said considering F0 and Q0 makes problem unstable and the total mass changes over the time. Also in some cases(different F0 and Q0) my Total_Mass_Growth over the time at 18000s is negative and I don't know how this can be happen.

请先登录，再进行评论。

Answer 2

Star Strider 2024-4-3

I do not completely understand your problem. If y9ou have values for all the variables as functions of time from the beginning to the end of the reaction, you can use the techniques in Parameter Estimation for a System of Differential Equations to estimate the parameters.

If you are fitting parameters with only two data points (the beginning and the end values for the reactions), it is only possible to fit two parameters (slope and intercept, for example) to them. Just use polyfit (and optionally polyval) for that.

4 个评论
显示 2更早的评论隐藏 2更早的评论

naiva saeedia 2024-4-5

@Star Strider but k_L_a_G_L_Light and k_L_a_G_L_Heavy are not kinetic constants. They are the mass transfer coefficient..But I guess I understood your point. I can do a linear regression as you said and then use the function to compute k_L_a_G_L_Light and k_L_a_G_L_Heavy from differential equations that I have.

naiva saeedia 2024-4-5

@Star Strider

but again there is one problem I just have the values at the start for x1,x7,x8,x9,x10,x11,x12,x13,x14 and I don't know thier final value at the end of the time so I guess I can't estimate these two parameters with these datas.

请先登录，再进行评论。

Answer 3

Sam Chak 2024-4-13

在 MATLAB Online 中打开

Hi @naiva saeedia

I have prepared a simple demonstration below to illustrate that the fmincon optimization solver can be utilized to find the optimal combinations of

and

, provided that the mass transfer system is stable and a solution set exists. Although the volumetric inflow rate

and molar outflow

of ethylene are now non-zero, it is crucial to ensure the conservation of flow.

In the event that the measured readings are indeed from a stable mass transfer system, but the mathematical model described in the paper generates unstable responses, we can only conclude that the math model is inaccurate and fails to properly capture certain flow dynamics. It is also possible that human errors occurred during the display or printing of the math model in the paper.

Regardless, it is not your fault because the theoretical concept has demonstrated the possibility of finding the optimal combinations of

and

, provided that the following conditions are satisfied:

the mass transfer system is asymptotically stable,
the math model is relatively accurate, and
a solution set exists.

Please let me know if you are satisfied with these explanations.

obj = @costfun;

lb = [0, 0.5];

ub = [0.5, 1];

k0 = (lb + ub)/2;

k = fmincon(obj, k0, [], [], [], [], lb, ub)

Local minimum possible. Constraints satisfied. fmincon stopped because the size of the current step is less than the value of the step size tolerance and constraints are satisfied to within the value of the constraint tolerance.

k = 1x2

0.2306 0.7502

tspan = [0 1.8e4];

x0 = [ones(7, 1); zeros(7, 1)];

[t, x] = ode45(@(t, x) massTransfer(t, x, k), tspan, x0);

plot(t, x), grid on

xlabel('t'), ylabel('\bf{x}\rm(t)'), title('System Response')

function J = costfun(param)

tspan = [0 1.8e4];

x0 = [ones(7, 1); zeros(7, 1)];

[t, x] = ode45(@(t, x) massTransfer(t, x, param), tspan, x0);

Moles_C4_G = x(end,2);

Moles_C6_G = x(end,3);

Moles_C8_G = x(end,4);

Moles_C10_G = x(end,5);

Moles_C12_G = x(end,6);

J1 = (Moles_C4_G - 0)^2;

J2 = (Moles_C6_G - 0.0168)^2;

J3 = (Moles_C8_G - 0.3376)^2;

J4 = (Moles_C10_G - 0.8776)^2;

J5 = (Moles_C12_G - 1.6810)^2;

J = J1 + J2 + J3 + J4 + J5;

end

function dNdt = massTransfer(t, x, k)

%% Parameters to be determined

kL = k(1); % k_L_a_G_L_Light, 1st parameter

kH = k(2); % k_L_a_G_L_Heavy, 2nd parameter

%% Constants

Q0 = 580; % Q_G_in_ml_min

Q1 = Q0*1.66667e-8; % Q_G_in (m3/s)

P1 = 36e5; % P_Ethylene

T1 = 230 + 273.15; % T_Ethylene

T2 = 230 + 273.15; % T_ref

R = 8.314; % Absolute Gas constant in the Universe

C1 = P1/(R*T1); % C_G_in

F0 = Q1*C1/1e-3; % F_G_out_mmol_s (original value 0.354)

F1 = F0*1e-3; % F_G_out (conservation of flow, inflow Q1*C1 = outflow F1)

Q2 = F1; % Q_G_in_C_G_in_C4

Q3 = F1; % Q_G_in_C_G_in_C6

Q4 = F1; % Q_G_in_C_G_in_C8

Q5 = F1; % Q_G_in_C_G_in_C10

Q6 = F1; % Q_G_in_C_G_in_C12

Q7 = F1; % Q_G_in_C_G_in_C11H24

VR = 300e-6; % V_R

VG = 250e-6; % V_G

VL = 50e-6; % V_L

K1 = 3.24; % K_L_G_Ethylene

K2 = 2.23; % K_L_G_Butene

K3 = 1.72; % K_L_G_Hexene

K4 = 1.3; % K_L_G_Octene

K5 = 0.985; % K_L_G_Decene

K6 = 0.751; % K_L_G_Dodecene

K7 = 0.842; % K_L_G_Undecane

wc = (0.3+0.25)*1e-3; % w_cat

k10 = 2.224e-4;

k20 = 1.533e-4;

k30 = 3.988e-5;

k40 = 1.914e-7;

k50 = 4.328e-5;

k60 = 2.506e-7;

k70 = 4.036e-5;

k80 = 1.062e-6;

k90 = 6.055e-7;

E1 = 109.53;

E2 = 15.229;

E3 = 7.881;

E4 = 44.454;

E5 = 9.438;

E6 = 8.426;

E7 = 10.911;

E8 = 12.537;

E9 = 7.127;

k1 = k10*exp(- E1*(1/T1 - 1/T2)/R);

k2 = k20*exp(- E2*(1/T1 - 1/T2)/R);

k3 = k30*exp(- E3*(1/T1 - 1/T2)/R);

k4 = k40*exp(- E4*(1/T1 - 1/T2)/R);

k5 = k50*exp(- E5*(1/T1 - 1/T2)/R);

k6 = k60*exp(- E6*(1/T1 - 1/T2)/R);

k7 = k70*exp(- E7*(1/T1 - 1/T2)/R);

k8 = k80*exp(- E8*(1/T1 - 1/T2)/R);

k9 = k90*exp(- E9*(1/T1 - 1/T2)/R);

%% ODEs

dN1 = Q1*C1 - F1 - VR*kL*(x(1)/VG - K1*x(8) /VL);

dN2 = Q2 - F1 - VR*kL*(x(2)/VG - K2*x(9) /VL);

dN3 = Q3 - F1 - VR*kL*(x(3)/VG - K3*x(10)/VL);

dN4 = Q4 - F1 - VR*kH*(x(4)/VG - K4*x(11)/VL);

dN5 = Q5 - F1 - VR*kH*(x(5)/VG - K5*x(12)/VL);

dN6 = Q6 - F1 - VR*kH*(x(6)/VG - K6*x(13)/VL);

dN7 = Q7 - F1 - VR*kH*(x(7)/VG - K7*x(14)/VL);

dN8 = VR*kL*(x(1)/VG - K1*x(8) /VL) + wc*(-2*k1*x(8)^2 /VL^2 - k2*x(8)*x(9) /VL^2 - k3*x(8)*x(10)/VL^2 - k5*x(8)*x(11)/VL^2 - k7*x(8)*x(12)/VL^2);

dN9 = VR*kL*(x(2)/VG - K2*x(9) /VL) + wc*( k1*x(8)^2 /VL^2 - k2*x(8)*x(9) /VL^2 - 2*k4*x(9)^2 /VL^2 - k6*x(9)*x(10)/VL^2 - k8*x(9)*x(11)/VL^2);

dN10= VR*kL*(x(3)/VG - K3*x(10)/VL) + wc*( k2*x(8)*x(9) /VL^2 - k3*x(8)*x(10)/VL^2 - k6*x(9)*x(10)/VL^2 - 2*k9*x(10)^2 /VL^2);

dN11= VR*kH*(x(4)/VG - K4*x(11)/VL) + wc*( k3*x(8)*x(10)/VL^2 + k4*x(9)^2 /VL^2 - k5*x(8)*x(11)/VL^2 - k8*x(9)*x(11)/VL^2);

dN12= VR*kH*(x(5)/VG - K5*x(12)/VL) + wc*( k5*x(8)*x(11)/VL^2 + k6*x(9)*x(10)/VL^2 - k7*x(8)*x(12)/VL^2);

dN13= VR*kH*(x(6)/VG - K6*x(13)/VL) + wc*( k7*x(8)*x(12)/VL^2 + k8*x(9)*x(10)/VL^2 + k9*x(10)^2 /VL^2);

dN14= VR*kH*(x(7)/VG - K7*x(14)/VL);

%% Group ODEs as a column vector

dNdt= [dN1;

dN2;

dN3;

dN4;

dN5;

dN6;

dN7;

dN8;

dN9;

dN10;

dN11;

dN12;

dN13;

dN14];

end

4 个评论
显示 2更早的评论隐藏 2更早的评论

naiva saeedia 2024-4-15

@Sam Chak

Hi and yes I can provide my derivation of equations. I share my notes with you. Feel free to ask question about them and thx for giving your time for solving my problem. I really appreciate that. Also the pictures sort based on the number.

Sam Chak 2024-4-15

编辑：Sam Chak 2024-4-15

Hi @naiva saeedia

After comparing your model and Woo's model, here is my analysis. Please note that my expertise is limited to High School Chemistry, and I have no knowledge of PhD-level chemical reactions. Therefore, it is essential for you to verify the analysis by running simulations using educated guesses for the parameters and testing the outcomes.

Gas Phase dynamics

Since

and

, then

Solving

for

The equilibrium for

is found to be a function of

Liquid Phase dynamics

Since

and

, then

Solving

for

If the product terms in R are neglected, then

The equilibrium for

is found to be a function of

Solving simultaneous equations

by substitutions

to obtain the approximate number of moles in both gas and liquid at steady state:

Conditions:

By the way, I noticed that you recently asked a similar question regarding the optimization problem, and @Torsten's fmincon solution was accepted. It's worth mentioning that the authors, Woo et al., claimed to have utilized the lsqcurvefit solver to estimate the kinetic parameters

and

.

请先登录，再进行评论。

I want to find two missing parameters in an ODE system of equations using regression/optimization in MATLAB

2 个评论
显示无隐藏无

回答（3 个）

29 个评论
显示 27更早的评论隐藏 27更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

I want to find two missing parameters in an ODE system of equations using regression/optimization in MATLAB

2 个评论 显示 无隐藏 无

回答（3 个）

29 个评论 显示 27更早的评论隐藏 27更早的评论

4 个评论 显示 2更早的评论隐藏 2更早的评论

4 个评论 显示 2更早的评论隐藏 2更早的评论

另请参阅

类别

标签

Community Treasure Hunt

2 个评论
显示无隐藏无

29 个评论
显示 27更早的评论隐藏 27更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论