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

5 次查看(过去 30 天)
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
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
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

请先登录,再进行评论。

回答(3 个)

Sam Chak
Sam Chak 2024-4-3
编辑:Sam Chak 2024-4-3
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 个评论
Sam Chak
Sam Chak 2024-4-13
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
naiva saeedia 2024-4-13
编辑:naiva saeedia 2024-4-13
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.

请先登录,再进行评论。


Star Strider
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 个评论
naiva saeedia
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
naiva saeedia 2024-4-5
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.

请先登录,再进行评论。


Sam Chak
Sam Chak 2024-4-13
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
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
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 个评论
naiva saeedia
naiva saeedia 2024-4-15
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
Sam Chak 2024-4-15
编辑:Sam Chak 2024-4-15
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 .

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by