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

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

[Alan Stevens]

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]

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]

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

[naiva saeedia]

Thx a lot for your answer..I have also changed n_initial_solvent value to 10 in my original code cause all of the C11H24_G_moles were negeative with the currnet value of n_Initial_Solvent..Yes It's a challenging task cause I don't have the experimental value for all of the x(1) to x(14)..Can you please check my new code and see if the values of k_L_a_G_L_Light and k_L_a_G_L_Heavy is changing or not? Also with your values I get negative value for moles_C4_G and moles_C6_G and It's not physically OK since moles(or x(end,i)) can't be negative..What you have done really helped me to understand the margines of k_L_a_G_L_Light and k_L_a_G_L_Heavy and I really appreciate that but how you calculate these two values for these two parameters? Can you please share your code? (I have to present the steps for the calculation of k_L_a_G_L to my professor so any further help would be greatly helpful)

[Sam Chak]

Hi @naiva saeedia, Please check my 'Update' message above. I have some queries about the system and the measured data.

[naiva saeedia]

Hi @Sam Chak, Thx for the quick reply. No I haven't conducted a stability analysis on the mathematical model of the system but the equations and parameters are simulated based on this paper at the attachments and I don't think the system itself is unstable. I have these datas for x2,x3,x4,x5,x6 based on Fig2 on the page 9 of this paper. I have used im2graph to collect the experimental values for C4 to C12 and based on this paper. the values of these datas measured at 18000s(5hours) and based on the paper the system is at steady-state condition because the gas outlet(F_G_out) is constant all the time. I hope this helps. Also I have tried to find values for k_L_a_G_L by hand and from what I have seen the x6 condition(x6(18000)=0.01351) is hard to achieve maybe you can ignore this condition and see what is happening?

[Sam Chak]

Hi @naiva saeedia

I have briefly reviewed the article. Since I am not familiar with mass balance modeling and the associated symbols, I cannot immediately determine whether the system described in the article is stable or unstable. However, it is possible to identify the equilibrium point by setting the derivative terms to zero and solving the resulting algebraic equations. Afterwards, the Jacobian matrix can be obtained by evaluating it at the equilibrium point, and the eigenvalues of that matrix can be computed.

In my opinion, the next step you should take (after performing the stability analysis) is to replicate the results presented in the article.

[naiva saeedia]

Hi @Sam Chak

I have tried to find the Eq points by setting the derivative terms to zero and solve the system of the algebric equations. However as you see in my photo the constant term of some equations is equal to zero and as a result I can't possibly solve this system of equations. What does this mean in your opinion? Can you please share your thoughts and knowledge with me? It would be really appericated.

[Torsten]

What's the reason for the parameter k in the last equation ?

However as you see in my photo the constant term of some equations is equal to zero and as a result I can't possibly solve this system of equations.

Why do you think so ?

[naiva saeedia]

@Torsten there is a tiny mistake..I will fix it rn .

[naiva saeedia]

@Torsten

I fixed the Equations sorry for that and thanks for pointing it out. That because in equation 1,9,10,11,12,13,14 we don't have any numerical values on the left side and the equations are equal to zero. If I remembered correctly when this is happening the system of Equations has no answer based on a theory in algebra. I used MAPLE to solve the system of equations and it gave me no answer..I guess it's going to be the same in MATLAB too(I haven't tried it on MATLAB)

[Torsten]

If I remembered correctly when this is happening the system of Equations has no answer based on a theory in algebra.

No, you don't remember correctly. x(7) = x(14) = 0,e.g., is a solution of the first equation.

Check whether MATLAB's "lsqnonlin" together with lower bound 0 for all x returns a solution of your quadratic system.

But in order to solve your problem, your algebraic equations should depend on the two parameters to be fitted, shouldn't they ? And you had to assume that your system is already in equilibrium at t = tend.

[naiva saeedia]

@Torsten

Thx...But how you can use Isqnonlin to solve the quadratic system? Based on what I see on MATLAB HELP it's for non-linear data fitting problems. I'm a little torn rn. Can you please explain what should I do with more details?

[naiva saeedia]

@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]

@Torsten

But in order to solve your problem, your algebraic equations should depend on the two parameters to be fitted, shouldn't they ? And you had to assume that your system is already in equilibrium at t = tend.

Yes and actually this is a great idea to find k_L_a_G_light and k_L_a_G_Heavy..If I understood correctly, your point is using isqnonlin and my equation at t=t_end to find the k_L_a_G_light and k_L_a_G_Heavy right?

[Torsten]

Write your equations as F(x) = 0 and call lsqnonlin:

f = @(x)[x(1)^2+x(2)^2-2;x(1)+x(2)-1];
x0 = [-0.3;0.5];
x = lsqnonlin(f,x0)

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

x = 2x1

-0.3660 1.3660

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

norm(f(x))

ans = 8.1179e-10

[Torsten]

If I understood correctly, your point is using isqnonlin and my equation at t=t_end to find the k_L_a_G_light and k_L_a_G_Heavy right?

Yes.

But note that the equilibrium will depend on your initial conditions for all the concentrations, but you wrote you only have initial conditions for some components of the x-vector. In this case, your attempt to fit the parameters doesn't make sense in my opinion. Or you must try to alter your model and exclude the components for which initial conditions are unknown.

[naiva saeedia]

@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]

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?

[naiva saeedia]

Hi @Sam Chak

Yes the equations are identical to the equations of the paper. The reason why Q_G_in*C_G_in is only presented in the first ODE is because it's only has value for Ethylene(X(1)) for other variables(x(2),...,x(14)) It's zero and I figured this out based on what the paper said about the feed and reactor condition..Also as you can see in the Fig2 legend we have C4,C6,C8,C10,C12 points on different cases(That's why there is so many points on the fig2) also it shows the product distribution and the product is only gas phase..I changed the equations format to simulated in Matlab. X(1)--->C2H4_gas X(2)--->C4H8_gas ... X(7)--->Undecane(solvent)_gas X(8)--->C2H4_Liquid X(9)--->C4H8_Liquid .... X(14)--->Undecane_Liquid So based on this and fig2 we have only C4H8_gas,C6H12_gas,C8H16_gas,C10H20_gas,C12H24_gas(Five different points for case1) I hope it wasn't confusing..My point is The equations are correct and I only have five initial values from fig2 not 14.

[Torsten]

@naiva saeedia

But if you found the initial conditions with assumed parameters for k_L_a_G_L_Light and k_L_a_G_L_Heavy, it doesn't make sense trying to find different parameters for k_L_a_G_L_Light and k_L_a_G_L_Heavy with these initial conditions afterwards, does it ?

[naiva saeedia]

@Torsten

Well technically you are right but I just want to find the order of magnityude of k_L_a_G_L_Light and k_L_a_G_L_Heavy..and I want to present something that MATLAB calculated rather than something that I have found by guess and at this point I ran out of solutions :(...It doesn't make sense but it's better than nothing I guess.

[Sam Chak]

@naiva saeedia

I'm running out of ideas. As a last resort, I suggest going back and reviewing the 14 differential equations presented in elegant mathematical notation. It's possible that you may extract some useful information or express the equilibrium in terms of and .

However, it's important to avoid the cluttered and excessively long equation form, as it becomes challenging for humans to decipher the code and grasp the mathematical concepts at a glance. I certainly cannot do.

[naiva saeedia]

@Sam Chak

Hi and sorry for late response..As you said previously I have computed Eigenvalues of the Jacobian matrix at steady state points. I guess the system of the equations is unstable and I can't figure out why the system is unstable. I just want to know your opinion And how I can wrote these equations in matlab without using long equation forms? Do you have any example of a neat MATLAB code for the system of ODE Equations?

[Sam Chak]

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

[naiva saeedia]

@Sam Chak

Hi again and thx for improving the code readability. I really appreciate that. about your first question yes I know we should solve the algebric equations to find Equilibrium points. I have solved the algebric equations with lsqnonlin. I have also computed Jacobian matrix at these points and calculated eigenvalues from it and I shared the results with you in the Eigenvalues.mat. as you see the eigen values are complex numbers and some the real part of these numbers are positive and some of them are negative and we have imaginary parts for some of the eigenvalues so from the mathematics the system is unstable. In theory after 18000s system should be at steady state. I don't have information about the thermodynamics of this system but from the mass transfer prospective system should be stable.

[Sam Chak]

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]

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]

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)

[naiva saeedia]

@Sam Chak

Thank you so much for explanation.

[Sam Chak]

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]

@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]

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.

[naiva saeedia]

No actually I don't have the end values for all variables..I have all the initial values for x1 to x14 and I have x2,x3,x4,x5,x6 values at the end..And I'm trying to estimate two unknown parameters in the system of the equations..I want to know if it's even possible for MATLAB to do that...What you meant by using polyval is fit a plot to my points? I have these datas for xs:

initial values:

x1=0.0083

x2=0

x3=0

...

x13=0

x14=10

The End values:

x2=0.07891

x3=0.14349

x4=0.063348

x5=0.02686

x6=0.01351

and the values for the rest of xs are unknown..I don't know if i was able to explain it correctly.

[Star Strider]

If you are only fitting two points, you are essentially doing a linear regression. The polyval function makes this easier. However if the beginning point is always zero with respect to both time and concentration, then the kinetic constant between the two points is simply: end_value/elapsed_time. The units are whatever you want them to be.

[naiva saeedia]

@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]

@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]

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

<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

[naiva saeedia]

@Sam Chak

Hi and thanks a lot for defining this code in Matlab. Your opinion about conversion of mass(F1 should be equal to Q1*C1) is really interesting. I have question due this system in mass transfer is a two phase semi-batch system. I have some liquid in the reactor and then I enter the ethylene as gas phase to the reactor. Reactions are happening in the liquid phase and then the products go to the gas phase and exit the reactor. My question is maybe the inlet mass trapped in the liquid phase and that's why F1 is not equal to Q1*C1 in the paper that I have shared with you. Outlet must be equal to Inlet in all cases??

[Sam Chak]

Hi @naiva saeedia

Could you please provide me with your own derivation of the mathematical model for the two-phase semi-batch system? This can be a lengthy task but this will enable me to compare it with Woo's model. If there are any discrepancies, we can investigate what might have gone wrong. If both models are correct, then it suggests that the dynamic model described in the MATLAB code is likely incorrect.

[naiva saeedia]

@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]

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 .