I want to optimize an ODE system of equations with constraints

Hi..I have a set of ODEs..I want to find the optimal value of kL and kH with these constraints:

0=<kL<=100

0=<kH<=100

Mass_C4_G-0.2985<0.1

Mass_C6_G-0.6498<0.1

Mass_C8_G-0.6147<0.1

Mass_C10_G-0.43917<0.1

Mass_C12_G-0.40398<0.1

Also all the values of Moles_Ci_G and Moles_Ci_L should be positive.

Can someone help me with the code?

clc
clear
close all
%%Parameters that I want to change
kL=0.1;      % k_L_a_G_L_Light, 1st parameter 
kH=0.1;   % k_L_a_G_L_Heavy, 2nd parameter 
%%Constants
Q0=100;            % 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_C11
Q7=0;              % Q_G_in_C_G_in_C12

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.0179;         % F_G_out_mmol_s
F1=F0.*1e-3;       % F_G_out
VR=300e-6;         %Volume of the reactor
VG=250e-6;         %Volume of the gas
VL=50e-6;          %Volume of the liquid

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;  %Catalyst weight

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

n1=C1.*VR;  %Ethylene initial condition(moles)
n2=10;      %Solvent initial condition(moles)

dNdt=@(t,x) [Q1.*C1-F1-VR.*kL.*(x(1)./VG-K1.*x(8)./VL);Q2-F1-VR.*kL.*(x(2)./VG-K2.*(x(9)./VL));Q3-F1-VR.*kL.*(x(3)./VG-K3.*(x(10)./VL));Q4-F1-VR.*kH.*(x(4)./VG-K4.*(x(11)./VL));Q5-F1-VR.*kH.*(x(5)./VG-K5.*x(12)./VL);Q6-F1-VR.*kH.*(x(6)./VG-K6.*x(13)./VL);Q7-F1-VR.*kH.*(x(7)./VG-K7.*x(14)./VL);VR.*kL.*(x(1)./VG-K1.*x(8)./VL)+wc.*(-2.*k1.*(x(8)./VL).^2-k2.*(x(8)./VL).*(x(9)./VL)-k3.*(x(8)./VL).*(x(10)./VL)-k5.*(x(8)./VL).*(x(11)./VL)-k7.*x(8).*x(12)./(VL.^2));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);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);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);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);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);VR.*kH.*(x(7)./VG-K7.*x(14)./VL)];
[t,x]=ode45(dNdt,[0,18000],[n1;0;0;0;0;0;0;0;0;0;0;0;0;n2]);
%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);
%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;

[Torsten]

Hi..I have a set of ODEs..I want to find the optimal value of kL and kH with these constraints:

You forgot to specify what you consider as "optimal" in your case.

[naiva saeedia]

@Torsten

Hi..I just want to know if this problem has any answers or not. I don't have a specific value for kL and kH in my mind

Answer 1

[Torsten]

I ran your code with k in the order of 0.1,1 and 10, but it seems the k-values have amost no influence on the result of the integration. Since the concentrations you get are far off the end concentrations you prescribed, "fmincon" tells you that no feasible solution for the k values could be found.

kL = 0.1;
kH = 0.1;
p0 = [kL,kH];
lb = [0;0];
ub = [100;100]; 
obj = @(p)0;
%[c,ceq] = nonlcon(p0)
sol = fmincon(obj,p0,[],[],[],[],lb,ub,@nonlcon);

function [c,ceq] = nonlcon(p)
%%Parameters that I want to change
kL=p(1);      % k_L_a_G_L_Light, 1st parameter 
kH=p(2);   % k_L_a_G_L_Heavy, 2nd parameter 
%%Constants
Q0=100;            % 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_C11
Q7=0;              % Q_G_in_C_G_in_C12
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.0179;         % F_G_out_mmol_s
F1=F0.*1e-3;       % F_G_out
VR=300e-6;         %Volume of the reactor
VG=250e-6;         %Volume of the gas
VL=50e-6;          %Volume of the liquid
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;  %Catalyst weight
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);
n1=C1.*VR;  %Ethylene initial condition(moles)
n2=10;      %Solvent initial condition(moles)
dNdt=@(t,x) [Q1.*C1-F1-VR.*kL.*(x(1)./VG-K1.*x(8)./VL);Q2-F1-VR.*kL.*(x(2)./VG-K2.*(x(9)./VL));Q3-F1-VR.*kL.*(x(3)./VG-K3.*(x(10)./VL));Q4-F1-VR.*kH.*(x(4)./VG-K4.*(x(11)./VL));Q5-F1-VR.*kH.*(x(5)./VG-K5.*x(12)./VL);Q6-F1-VR.*kH.*(x(6)./VG-K6.*x(13)./VL);Q7-F1-VR.*kH.*(x(7)./VG-K7.*x(14)./VL);VR.*kL.*(x(1)./VG-K1.*x(8)./VL)+wc.*(-2.*k1.*(x(8)./VL).^2-k2.*(x(8)./VL).*(x(9)./VL)-k3.*(x(8)./VL).*(x(10)./VL)-k5.*(x(8)./VL).*(x(11)./VL)-k7.*x(8).*x(12)./(VL.^2));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);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);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);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);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);VR.*kH.*(x(7)./VG-K7.*x(14)./VL)];
[t,x]=ode45(dNdt,[0,18000],[n1;0;0;0;0;0;0;0;0;0;0;0;0;n2]);
%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);
%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;
c(1) = Mass_C4_G-0.2985-0.1;
c(2) = -(Mass_C4_G-0.2985)-0.1;
c(3) = Mass_C6_G-0.6498-0.1;
c(4) = -(Mass_C6_G-0.6498)-0.1;
c(5) = Mass_C8_G-0.6147-0.1;
c(6) = -(Mass_C8_G-0.6147)-0.1;
c(7) = Mass_C10_G-0.43917-0.1;
c(8) = -(Mass_C10_G-0.43917)-0.1;
c(9) = Mass_C12_G-0.40398-0.1;
c(10) = -(Mass_C12_G-0.40398)-0.1;
ceq = [];
end

[naiva saeedia]

@Torsten

Hi and thx for simulating this optimization. Can we use a different approach that @William Rose has described? I want to minimize MTE in this function using fmincon() and the accepted range for kL and kH. Could you please help me to do this?

function MTE=MassTransferError(kL,kH)
clc
clear
close all 
%%Constants
Q0=100;            % 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_C11
Q7=0;              % Q_G_in_C_G_in_C12

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.0179;         % F_G_out_mmol_s
F1=F0.*1e-3;       % F_G_out
VR=300e-6;         %Volume of the reactor
VG=250e-6;         %Volume of the gas
VL=50e-6;          %Volume of the liquid

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;  %Catalyst weight

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

n1=C1.*VR;  %Ethylene initial condition(moles)
n2=10;      %Solvent initial condition(moles)

dNdt=@(t,x) [Q1.*C1-F1-VR.*kL.*(x(1)./VG-K1.*x(8)./VL);Q2-F1-VR.*kL.*(x(2)./VG-K2.*(x(9)./VL));Q3-F1-VR.*kL.*(x(3)./VG-K3.*(x(10)./VL));Q4-F1-VR.*kH.*(x(4)./VG-K4.*(x(11)./VL));Q5-F1-VR.*kH.*(x(5)./VG-K5.*x(12)./VL);Q6-F1-VR.*kH.*(x(6)./VG-K6.*x(13)./VL);Q7-F1-VR.*kH.*(x(7)./VG-K7.*x(14)./VL);VR.*kL.*(x(1)./VG-K1.*x(8)./VL)+wc.*(-2.*k1.*(x(8)./VL).^2-k2.*(x(8)./VL).*(x(9)./VL)-k3.*(x(8)./VL).*(x(10)./VL)-k5.*(x(8)./VL).*(x(11)./VL)-k7.*x(8).*x(12)./(VL.^2));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);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);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);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);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);VR.*kH.*(x(7)./VG-K7.*x(14)./VL)];
[t,x]=ode45(dNdt,[0,18000],[n1;0;0;0;0;0;0;0;0;0;0;0;0;n2]);
%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);
%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;
MTE=(abs(Mass_C4_G-0.2985))+(abs(Mass_C6_G-0.6498))+(abs(Mass_C8_G-0.6147))+(abs(Mass_C10_G-0.43917));
end

[Torsten]

function MTE=MassTransferError(p)
  kL = p(1);
  kH = p(2);
...

instead of

function MTE=MassTransferError(kL,kH)
...
    

To make your objective differentiable, better use

MTE = (Mass_C4_G-0.2985)^2+(Mass_C6_G-0.6498)^2+(Mass_C8_G-0.6147)^2+(Mass_C10_G-0.43917)^2;

instead of

MTE=(abs(Mass_C4_G-0.2985))+(abs(Mass_C6_G-0.6498))+(abs(Mass_C8_G-0.6147))+(abs(Mass_C10_G-0.43917));

[naiva saeedia]

@Torsten

Tnq so much.

[naiva saeedia]

@Torsten

Hi again. I have tried to run fmincon with the changes in my function but with this code:

Objective=@MassTransferError;
k0=[1;1];
A = [];
b = [];
Aeq = [];
beq = [];
lb=[0 ; 0];
ub=[1000;1000];
k = fmincon(Objective, k0, A, b, Aeq, beq, lb, ub);
disp(k)

I got this error:

Reference to a cleared variable p.

Error in MassTransferError (line 6)

kL = p(1);

Error in fmincon (line 568)

initVals.f = feval(funfcn{3},X,varargin{:});

Error in Mass_Error_Optimization (line 9)

k = fmincon(Objective, k0, A, b, Aeq, beq, lb, ub);

Caused by:

Failure in initial objective function evaluation. FMINCON cannot continue.

Can you please help me. why is this happening?

my function code:

function MTE=MassTransferError(p)
clc
clear
close all 
%%Constants
kL = p(1);
kH = p(2);
Q0=100;            % 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_C11
Q7=0;              % Q_G_in_C_G_in_C12

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.0179;         % F_G_out_mmol_s
F1=F0.*1e-3;       % F_G_out
VR=300e-6;         %Volume of the reactor
VG=250e-6;         %Volume of the gas
VL=50e-6;          %Volume of the liquid

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;  %Catalyst weight

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

n1=C1.*VR;  %Ethylene initial condition(moles)
n2=10;      %Solvent initial condition(moles)

dNdt=@(t,x) [Q1.*C1-F1-VR.*kL.*(x(1)./VG-K1.*x(8)./VL);Q2-F1-VR.*kL.*(x(2)./VG-K2.*(x(9)./VL));Q3-F1-VR.*kL.*(x(3)./VG-K3.*(x(10)./VL));Q4-F1-VR.*kH.*(x(4)./VG-K4.*(x(11)./VL));Q5-F1-VR.*kH.*(x(5)./VG-K5.*x(12)./VL);Q6-F1-VR.*kH.*(x(6)./VG-K6.*x(13)./VL);Q7-F1-VR.*kH.*(x(7)./VG-K7.*x(14)./VL);VR.*kL.*(x(1)./VG-K1.*x(8)./VL)+wc.*(-2.*k1.*(x(8)./VL).^2-k2.*(x(8)./VL).*(x(9)./VL)-k3.*(x(8)./VL).*(x(10)./VL)-k5.*(x(8)./VL).*(x(11)./VL)-k7.*x(8).*x(12)./(VL.^2));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);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);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);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);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);VR.*kH.*(x(7)./VG-K7.*x(14)./VL)];
[t,x]=ode45(dNdt,[0,18000],[n1;0;0;0;0;0;0;0;0;0;0;0;0;n2]);
%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);
%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;
MTE = (Mass_C4_G-0.2985)^2+(Mass_C6_G-0.6498)^2+(Mass_C8_G-0.6147)^2+(Mass_C10_G-0.43917)^2;
end

Update:I have solved it. I shouldn't have used clear() on my function.

[naiva saeedia]

@Torsten

Hi again. I have another question. I have managed to find the kL and kH values that makes the MTE minimum but the problem is these values are making some of the moles_Ci_G and moles_Ci_L negative and I want fmincon minimize this objective function with all the moles_Ci_G and moles_Ci_L being positive. How can I define this in nonlcon function??

Update:I have wrote this function for my nlcon function. Is this correct?

function [c, ceq] = nlcon(p)
kL = p(1);
kH = p(2);
Q0=100;            % 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_C11
Q7=0;              % Q_G_in_C_G_in_C12

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.0179;         % F_G_out_mmol_s
F1=F0.*1e-3;       % F_G_out
VR=300e-6;         %Volume of the reactor
VG=250e-6;         %Volume of the gas
VL=50e-6;          %Volume of the liquid

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;  %Catalyst weight

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

n1=C1.*VR;  %Ethylene initial condition(moles)
n2=10;      %Solvent initial condition(moles)

dNdt=@(t,x) [Q1.*C1-F1-VR.*kL.*(x(1)./VG-K1.*x(8)./VL);Q2-F1-VR.*kL.*(x(2)./VG-K2.*(x(9)./VL));Q3-F1-VR.*kL.*(x(3)./VG-K3.*(x(10)./VL));Q4-F1-VR.*kH.*(x(4)./VG-K4.*(x(11)./VL));Q5-F1-VR.*kH.*(x(5)./VG-K5.*x(12)./VL);Q6-F1-VR.*kH.*(x(6)./VG-K6.*x(13)./VL);Q7-F1-VR.*kH.*(x(7)./VG-K7.*x(14)./VL);VR.*kL.*(x(1)./VG-K1.*x(8)./VL)+wc.*(-2.*k1.*(x(8)./VL).^2-k2.*(x(8)./VL).*(x(9)./VL)-k3.*(x(8)./VL).*(x(10)./VL)-k5.*(x(8)./VL).*(x(11)./VL)-k7.*x(8).*x(12)./(VL.^2));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);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);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);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);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);VR.*kH.*(x(7)./VG-K7.*x(14)./VL)];
[t,x]=ode45(dNdt,[0,18000],[n1;0;0;0;0;0;0;0;0;0;0;0;0;n2]);
%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);

c(1) = -(Moles_C2_G);
c(2) = -(Moles_C4_G);
c(3) = -(Moles_C6_G);
c(4) = -(Moles_C8_G);
c(5) = -(Moles_C10_G);
c(6) = -(Moles_C12_G);
c(7) = -(Moles_C11H24_G);
c(8) = -(Moles_C2_L);
c(9) = -(Moles_C4_L);
c(10) = -(Moles_C6_L);
c(11) = -(Moles_C8_L);
c(12) = -(Moles_C10_L);
c(13) = -(Moles_C12_L);
c(14) = -(Moles_C11H24_L);

ceq=[];

end

[Torsten]

I have wrote this function for my nlcon function. Is this correct?

Looks correct.

[naiva saeedia]

@Torsten

Thx a lot🙏🏻

Answer 2

[William Rose]

@naiva saeedia,

You stated goals such as

Mass_C4_G-0.2985<0.1

I assume you really would like

abs(Mass_C4_G-0.2985)<0.

and likewise for C6, C8, C10, C12.

I put your code into a while loop. The loop runs until conditions 1 and 2 are met, or after 50 passes, whichever comes first. Inside each loop pass, the following happens:

• Make guesses for 0.1<kL<100 and 0.1<kH<100.
• Run simulation.
• Check condition 1: all mole quantities are positive
• Check condition 2: five inequality constraints are met: abs(Mass...)<0.1.
• Compute massError=sum of absolute deviations of 5 masses from their target values

Results: All 50 guesses for (kL,kH) satisfied condition 1. None of the 50 guesses satisfied condition 2. The histogram of massError, for the 50 trials, is below. The histogram shows that the massError was between 428.4 and 428.75, for all the trials, and was between 428.7 and 428.75 for 27 out of 50 trials. In other words, the massError did not vary much, despite the wide range of random guesses for kL and kH.

So I aded code to save the five individual mass errors on each pass, to determine which mass(es) are causing the (total) massError to be large. Then I re-ran the code with 100 passes, instead of 50. This set of 100 random (kL,kH) pairs is independent of the previous 50 pairs. See attached code.

Results: All 100 guesses for (kL,kH) satisfied condition 1. None of the 100 guesses satisfied condition 2.

I displayed the minimum, mean, and maximum error for all five species: C4_G, C6_G, C8_G, C10_G, C12_G, with this command

disp([min(results(:,6:10));mean(results(:,6:10));max(results(:,6:10))]);

1.6875 5.6660 5.3837 5.6904 410.0246

1.6935 5.6738 5.3963 5.7076 410.2164

1.6946 5.6753 5.3995 5.7110 410.2561

As you can see from the printout above, the errors for the individual species did not change much, and the mass error for C12_G was the largest, by far. This makes me pessimistic that you will find values for kL and kH that satisfy your constraints. Also interesting is that fact that all individual mass errors are positive.

Next steps I recommend: Put the simulation inside a function whose output is massError. Use fmincon() to find the values of kL and kH that minimize massError, within the allowed range for kL and kH. See if the individual constraints on the 5 masses are met at this minimum.

clear
c1=false;   % condition 1: all masses>0 at end
c2=false;   % condition 2: five masses within 0.1 of target values 
i=0;        % pass counter
while ~(c1 & c2) & i<20 
%%Parameters that I want to change
kL=10^(-1+3*rand());             % k_L_a_G_L_Light, 1st parameter 
kH=10^(-1+3*rand());             % k_L_a_G_L_Heavy, 2nd parameter 
%%Constants
Q0=100;            % 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_C11
Q7=0;              % Q_G_in_C_G_in_C12
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.0179;         % F_G_out_mmol_s
F1=F0.*1e-3;       % F_G_out
VR=300e-6;         %Volume of the reactor
VG=250e-6;         %Volume of the gas
VL=50e-6;          %Volume of the liquid
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;  %Catalyst weight
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);
n1=C1.*VR;  %Ethylene initial condition(moles)
n2=10;      %Solvent initial condition(moles)
dNdt=@(t,x) [Q1.*C1-F1-VR.*kL.*(x(1)./VG-K1.*x(8)./VL);Q2-F1-VR.*kL.*(x(2)./VG-K2.*(x(9)./VL));Q3-F1-VR.*kL.*(x(3)./VG-K3.*(x(10)./VL));Q4-F1-VR.*kH.*(x(4)./VG-K4.*(x(11)./VL));Q5-F1-VR.*kH.*(x(5)./VG-K5.*x(12)./VL);Q6-F1-VR.*kH.*(x(6)./VG-K6.*x(13)./VL);Q7-F1-VR.*kH.*(x(7)./VG-K7.*x(14)./VL);VR.*kL.*(x(1)./VG-K1.*x(8)./VL)+wc.*(-2.*k1.*(x(8)./VL).^2-k2.*(x(8)./VL).*(x(9)./VL)-k3.*(x(8)./VL).*(x(10)./VL)-k5.*(x(8)./VL).*(x(11)./VL)-k7.*x(8).*x(12)./(VL.^2));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);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);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);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);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);VR.*kH.*(x(7)./VG-K7.*x(14)./VL)];
[t,x]=ode45(dNdt,[0,18000],[n1;0;0;0;0;0;0;0;0;0;0;0;0;n2]);
%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);
%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;

c1=length(x(end,(x(end,:)<0)))==0;   % condition 1: all masses positive at end
c2=(abs(Mass_C4_G-0.2985)<0.1) &...
   (abs(Mass_C6_G-0.6498)<0.1) &...
   (abs(Mass_C8_G-0.6147)<0.1) &...
   (abs(Mass_C10_G-0.43917)<0.1) &...
   (abs(Mass_C12_G-0.40398)<0.1);  % condition 2: five masses within 0.1 of target values
fprintf(' %i',i)
end         % end while loop
fprintf('\nkL=%.2f, kH=%.2f.\n',kL,kH)

You replied to @Torsten: "I just want to know if this problem has any answers or not. I don't have a specific value for kL and kH in my mind"

[naiva saeedia]

@William Rose

Tnq so much for doing this amazing work. Can we ignore C_12_G constraint since it's causing so much trouble and tbh it's not that important for me in this simulation? Also I haven't use fmincon() a lot. I'm going to try to do as you said and I'm also going to ask from @Torsten to help me with that too.

[naiva saeedia]

@William Rose

Hi again I have tried to use fmincon by myself but my problem is defining the objective function..I defined this as my objective function:

function MTE=MassTransferError(kL,kH)
clc
clear
close all 
%%Constants
Q0=100;            % 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_C11
Q7=0;              % Q_G_in_C_G_in_C12

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.0179;         % F_G_out_mmol_s
F1=F0.*1e-3;       % F_G_out
VR=300e-6;         %Volume of the reactor
VG=250e-6;         %Volume of the gas
VL=50e-6;          %Volume of the liquid

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;  %Catalyst weight

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

n1=C1.*VR;  %Ethylene initial condition(moles)
n2=10;      %Solvent initial condition(moles)

dNdt=@(t,x) [Q1.*C1-F1-VR.*kL.*(x(1)./VG-K1.*x(8)./VL);Q2-F1-VR.*kL.*(x(2)./VG-K2.*(x(9)./VL));Q3-F1-VR.*kL.*(x(3)./VG-K3.*(x(10)./VL));Q4-F1-VR.*kH.*(x(4)./VG-K4.*(x(11)./VL));Q5-F1-VR.*kH.*(x(5)./VG-K5.*x(12)./VL);Q6-F1-VR.*kH.*(x(6)./VG-K6.*x(13)./VL);Q7-F1-VR.*kH.*(x(7)./VG-K7.*x(14)./VL);VR.*kL.*(x(1)./VG-K1.*x(8)./VL)+wc.*(-2.*k1.*(x(8)./VL).^2-k2.*(x(8)./VL).*(x(9)./VL)-k3.*(x(8)./VL).*(x(10)./VL)-k5.*(x(8)./VL).*(x(11)./VL)-k7.*x(8).*x(12)./(VL.^2));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);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);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);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);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);VR.*kH.*(x(7)./VG-K7.*x(14)./VL)];
[t,x]=ode45(dNdt,[0,18000],[n1;0;0;0;0;0;0;0;0;0;0;0;0;n2]);
%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);
%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;
MTE=(abs(Mass_C4_G-0.2985))+(abs(Mass_C6_G-0.6498))+(abs(Mass_C8_G-0.6147))+(abs(Mass_C10_G-0.43917));
end

but my problem is when I ran fmincon() with this code:

Objective=@MassTransferError;
k0=[1;1];
A = [];
b = [];
Aeq = [];
beq = [];
lb=[0 ; 0];
ub=[1000;1000];
k = fmincon(Objective, k0, A, b, Aeq, beq, lb, ub);
disp(k)

I get this error:

Unrecognized function or variable 'kL'.

Caused by:

Failure in initial objective function evaluation. FMINCON cannot continue.

What should I do? Can you please help me?

Update: I have solved the issues. Thank you so much for giving me valuable suggestion.

[William Rose]

@naiva saeedia, you are welcome. @Torsten always gives excellent advice. Including his recommendation to use an error function for fmincon() that involves the squared errors, instead of the sum of the absolute values of the errors. It is easier for fmincon() to minimize an error function that is differentiable.

[naiva saeedia]

@William Rose

Yeah @Torsten is amazing. I have another problem though could you please check my last comment that I wrote for @Torsten? I want to know how can I make moles_Ci_G and moles_Ci_L positive in fmincon.

[William Rose]

@naiva saeedia,

It is interesting that you are getting negative masses for some species. With what values for kL and kH? I never got negative moles for any species, in the 150 simulations I did, using random (kL,kH) pairs in the range 0.1 to 100.

My initial thought is to check the equations for the reaction kinetics. Negative masses suggest an error in the equations. There could be an error in the "raw" chemical reaction equations, or in the code that represents those equations in nfunction dNdt.

You could compute the sum of all masses at each instant, after ode45() has completed, to see if the total mass is constant across time. If your equations include all reactants and products - and I don't know if they do - then the total mass should remain constant over time. If the total mass changes over time, then you know something is wrong. I have used this strategy to check for errors in complex cardiovascular simulations, where the total blood volume in all compartments should remain constant, in a simulated circulation with no leaks.

Do G and L refer to liquid and gas phase?

[naiva saeedia]

@William Rose

Well I didn't get negative values in this case but I have 6 other cases when I should change Q0 and F0 and n2 values. In one of the cases with Q0=350, F0=0.169 and n2 =3.05 I got negative value with kL=0.002158233, kH=0.0034029.(with the constraint function nlcon(p) that I wrote in @Torsten comment ). You think I should compute total mass in different times to see if it remains constant or not? That's actually a great idea. Biut the problem is I have inlet and outlet in my reactor. What you said is true for batch systems where there is no inlet and outlet I think. Yes it's liquid phase and gas phase.

[William Rose]

@naiva saeedia,

I looked at your differential equations. I see that they do not conserve mass, because:

dx1/dt = Q1*C1 - F1 - ...

dx2/dt = Q2 - F1 - ...

...

dx7/dt = Q7 - F1 - ...

where Q1*C1=1.43E-3 represents a constant positive influx of gas-phase ethylene that is not balanced by a corresponding loss anywhere, and -F1=-1.79E-5 is a constant rate of loss of moles that affects all 7 gas phase species, and is not balanced by a gain anywhere. What physical process removes moles of all 7 gas phase species at a constant rate, independent of their concentration? I wonder if the -F1 term is the cause of the negative mass values you have observed. (I have not yet observed negative mass values in the simulation.)

Influx rates Q2 through Q7 are zero, so those terms do not upset the mass balance.

When I look at the differential equations in your code, I also notice:

The differential equations include terms for exchange between the gas and liquid phase for all species. These equations are balanced: x1 with x8, x2 with x9, ..., x7 with x14.

The gas phase species, x1 through x7, and x14 (liquid phase undecane) only change due to gas-liquid exchange, and the influx and efflux terms mentioned above.

Species x8 through x13 (all liquid phase) also participate in catalyzed chemical reactions with one another. I have not carefully checked the differential equations for x8 through x13 to see if everything balances.

[William Rose]

@naiva saeedia,

I added code to compute and plot the total mass of all 14 species as a function of time. I ran the simuation with [kL,kH]=[0.1,0.1] and [1,1] and [10,10] and [99,99]. The total mass versus time plot was identical, to my eye, for all four simulations. The plots show that the total mass increases at an approximately constant rate throughout the simulaiton. The code to compute total mass (after ode45() has run) and make the plot is

% Compute and plot total mass for last simulation pass
massPerMole=[28;56;73;112;154;168;156;28;56;73;112;154;168;156];
massTot=x*massPerMole;
figure;
plot(t,massTot,'-r');
xlabel('Time'); ylabel('Mass'); grid on
title(['Total Mass, kL=',num2str(kL),' kH=',num2str(kH)]);

The plots for kL,kH=0.1,0.1 and for kL,kH=99,99 are below.

The observed average rate of total mass growth is

(massTot(end)-massTot(1))/(t(end)-t(1))=0.0271

The rate of mass growth we expect, due to the Q1*C1 gain in species 1, and the -F1 loss in species 1 through 7 (see my previous comment), is

28*Q1*C1-F1*sum(massPerMole(1:7))=0.0268

The observed and expected rate of mass growth are quite similar, but not identical. The small difference may reflect an imbalance in the differential equations for the reactions.

[naiva saeedia]

@William Rose

Hi. Sorry for late response. I have simulated the equations of this paper. It's for ethylene oligomerization and as you see in the first pic the paper said the GasOutlet(F0) should remain constant and the paper gave the values for each case in the second pic. Also the values of(Q0) aka ethylene feeding rate are given in the third pic and I calculated C1 by assuming ethylene as an ideal gas.

[naiva saeedia]

@William Rose

Thx for your great analysis about over all mass growth. It's really giving me some insights. I have a question due. How can I solve this problem and make these two values Identical. based on the paper I shouldn't change the value of neither Q0 and F0 but since Q0 is ethylene feeding rate I can change that imo. do you thing changing Q0 value is going to solve anything? and in some cases(for example case7) my Total_Mass_Growth and Observed_Mass_Growth is negative. Does this make sense in your opinion?(I have found kL and kH values with my optimization code) I'm going to share the case 7 code for you to see it.

clc
clear
close all
%%Parameters that I want to change
kL=0.601479349419186;      % k_L_a_G_L_Light, 1st parameter 
kH=0.044;   % k_L_a_G_L_Heavy, 2nd parameter 
%%Constants
Q0=400;            % 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_C11
Q7=0;              % Q_G_in_C_G_in_C12

P1=36e5;           % P_Ethylene
T1=180+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.291;          % F_G_out_mmol_s
F1=F0.*1e-3;       % F_G_out
VR=300e-6;         %Volume of the reactor
VG=250e-6;         %Volume of the gas
VL=50e-6;          %Volume of the liquid

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;  %Catalyst weight

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

n1=C1.*VR;  %Ethylene initial condition(moles)
n2=5.24;      %Solvent initial condition(moles)

dNdt=@(t,x) [Q1.*C1-F1-VR.*kL.*(x(1)./VG-K1.*x(8)./VL);Q2-F1-VR.*kL.*(x(2)./VG-K2.*(x(9)./VL));Q3-F1-VR.*kL.*(x(3)./VG-K3.*(x(10)./VL));Q4-F1-VR.*kH.*(x(4)./VG-K4.*(x(11)./VL));Q5-F1-VR.*kH.*(x(5)./VG-K5.*x(12)./VL);Q6-F1-VR.*kH.*(x(6)./VG-K6.*x(13)./VL);Q7-F1-VR.*kH.*(x(7)./VG-K7.*x(14)./VL);VR.*kL.*(x(1)./VG-K1.*x(8)./VL)+wc.*(-2.*k1.*(x(8)./VL).^2-k2.*(x(8)./VL).*(x(9)./VL)-k3.*(x(8)./VL).*(x(10)./VL)-k5.*(x(8)./VL).*(x(11)./VL)-k7.*x(8).*x(12)./(VL.^2));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);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);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);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);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);VR.*kH.*(x(7)./VG-K7.*x(14)./VL)];
[t,x]=ode45(dNdt,[0,18000],[n1;0;0;0;0;0;0;0;0;0;0;0;0;n2]);
%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);
%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;
%Mass Fractions
Total_Mass=Mass_C2_G+Mass_C4_G+Mass_C6_G+Mass_C8_G+Mass_C10_G+Mass_C12_G+Mass_C11H24_G;
C2_Fraction=Mass_C2_G./Total_Mass;
C4_Fraction=Mass_C4_G./Total_Mass;
C6_Fraction=Mass_C6_G./Total_Mass;
C8_Fraction=Mass_C8_G./Total_Mass;
C10_Fraction=Mass_C10_G./Total_Mass;
C12_Fraction=Mass_C12_G./Total_Mass;
C11H24_Fraction=Mass_C11H24_G./Total_Mass;
T1=table(C2_Fraction,C4_Fraction,C6_Fraction,C8_Fraction,C10_Fraction,C12_Fraction,C11H24_Fraction);
T2=table(Mass_C2_G,Mass_C4_G,Mass_C6_G,Mass_C8_G,Mass_C10_G,Mass_C12_G,Mass_C11H24_G);
T3=table(Mass_C2_L,Mass_C4_L,Mass_C6_L,Mass_C8_L,Mass_C10_L,Mass_C12_L,Mass_C11H24_L);
Hydrocarbons={'C2H4';'C4H8';'C6H12';'C8H16';'C10H20';'C12H24';'Solvent'};
LiquidPhase_g=[Mass_C2_L;Mass_C4_L;Mass_C6_L;Mass_C8_L;Mass_C10_L;Mass_C12_L;Mass_C11H24_L];
GasPhase_g=[Mass_C2_G;Mass_C4_G;Mass_C6_G;Mass_C8_G;Mass_C10_G;Mass_C12_G;Mass_C11H24_G];
T=table(GasPhase_g,LiquidPhase_g,'RowNames',Hydrocarbons)

T = 7x2 table
               GasPhase_g    LiquidPhase_g
               __________    _____________

    C2H4         2.8032         0.15844   
    C4H8         4.0551         0.36571   
    C6H12         10.94          1.2755   
    C8H16        6.9811           1.169   
    C10H20       3.3425         0.85102   
    C12H24       108.85          29.271   
    Solvent     0.22044        0.091559   

T4=table(Moles_C2_G,Moles_C4_G,Moles_C6_G,Moles_C8_G,Moles_C10_G,Moles_C12_G,Moles_C11H24_G);
T5=table(Moles_C2_L,Moles_C4_L,Moles_C6_L,Moles_C8_L,Moles_C10_L,Moles_C12_L,Moles_C11H24_L);
% Compute and plot total mass for last simulation pass
massPerMole=[28;56;73;112;154;168;156;28;56;73;112;154;168;156];
massTot=x*massPerMole;
Total_Mass_Growth = (massTot(end)-massTot(1))/(t(end)-t(1));
Observed_Mass_Growth = 28*Q1*C1-F1*sum(massPerMole(1:7));
error = abs(Total_Mass_Growth-Observed_Mass_Growth)/Observed_Mass_Growth

error = -0.0670

T = table(Total_Mass_Growth, Observed_Mass_Growth)

T = 1x2 table
    Total_Mass_Growth    Observed_Mass_Growth
    _________________    ____________________

        -0.036394             -0.039009      

[William Rose]

@naiva saeedia,

Thanks for the paper by Woo et al., 2021. Is it your goal to implement the model of Woo et al.?

You wrote "How can I solve this problem and make these two values [meaning inflow and outflow, I think] identical."

Before I attempt to answer that question, I have a few suggesitons:

Molecular weights for hexene and decene appear to be incorrect. This is a likely cause for why the observed rate of change of mass is not equal to the value expected. Incorrect masses mean that mass is not conserved in reactions involving these species. See comments I added below:

Mass_C6_G=Moles_C6_G.*73;       % should be C6H12 = 84
Mass_C10_G=Moles_C10_G.*154;    % should be C10H20=140

Document your code by including units for various quantities. This will make it easier to understand the code, and easier to notice possible mistakes.

You want to assure that inflow=outflow. Inflow and outflow, in this model, have units of total moles per second. Woo et al. say "The gas outlet flow was calculated so that the total number of moles in the gas phase was maintained as the initial number of moles to fix the gas pressure at a specified value..." You could try to implement equation 3 dynamically in your model. In this case, the outflow rate, Fout, would be another variable in the state vector x. Its units will be moles/s. You can write a differential equation for Fout. Fout is the TOTAL moles out per second. The effect of Fout on the individual gas species is proportional to the amount of each species in the gas, as a fraction of all the species in the gas. Therefore, there will be a term in the differential equation for ethylene (x1), corresponding to the effect of Fout, as follows:

dx(1)/dt = ... -Fout*x1/(sum(x(1:7)) + ...

and so on for the other gas species.

Think about what would be a good differential equation for Fout.

[naiva saeedia]

@William Rose

Hi. Yes my goal is implement woo et al's model in Matlab. Sorry for the mistakes in C10H20 and C6H12 molecular weights and thx for pointing this mistake out. I really appreciate that. What you wrote about Fout is pretty reasonable but I don't really know how can I have a differential equation for Fout. Fout is total moles out per second but the problem is I don't know the rate of total moles in woo et al's model and If I remember correctly the woo et al's model doesn't give specific information for this. Can't I just use the values in table 6 in the paper and add the *x1/(sum(x(1:7)) term to it in my equations?

[William Rose]

@naiva saeedia,

Your code includes

Q6=0;              % Q_G_in_C_G_in_C11
Q7=0;              % Q_G_in_C_G_in_C12

which suggests that x6=undecane and x7=dodecene. But your code also includes

K6=0.751;          % K_L_G_Dodecene
K7=0.842;          % K_L_G_Undecane

and

Moles_C12_G=x(end,6);
Moles_C11H24_G=x(end,7);

which means x6=dodecene and x7=undecane. I recommend you check the code carefully to make sure the equations and constants are correct and consistent for dodecene and undecane.

Woo et al. say "The gas outlet flow was calculated so that the total number of moles in the gas phase was maintained as the initial number of moles to fix the gas pressure at a specified value..." To make an outflow variable that adjusts dynamically to do what Woo et al say, proceed as follows:

Compute the initial number of gas phase moles

nTot0=n1+n2;

Let = gas outflow rate (mol/s). The initial conditon for x15 is that it equals the initial inflow rate:

x150=Q1*C1;  % initial outflow=initial inflow (mol/s)

The differential equation for outflow rate, x15, is

where = time constant for outflow adjustment (s). I think this will work, but you have to try it to be sure. The idea of the equation above is ... (more later)

[naiva saeedia]

Hi @William Rose

Thx for pointing this mistake out. I should change the comment section cause x6 is dodecene and x7 is undecane and, yeah, I'm going to check my equations rn. I have a question, though n2 is not the gas phase moles of undecane. It's the initial condition of undecane that's in the liquid phase, so I guess nTot0 should be equal to n1. I would love to hear the idea behind the equation and thx again for helping me and giving me your time.

[William Rose]

@naiva saeedia,

I have quesitons about the initial conditions and the gas-liquid balance.

n1 = C1*VR; % initial ethylene in gas (mol) 
n2= 10; % initial undecane in liquid (mol)

I think it should be n1=C1*VG and n2=36.63/156 (p.504 says 36.63 g undacane in reactor initially).

I would use different notation, in order to be consistent with other usage in the script.

When I tried introducing an equation for dx15/dt, the integration bew up after about 0.1 second. I am not sure why yet. It is possible that the problem is due to initial conditions that are very far away from equilibrium. When this is the case, there can be inital values for the derivative that are very vaery large, and this can cause instability in the nmerical solution. Therefore it is a good idea to make the initial conditions somewhere close to an equilibrium set of values. I noted above that n2=initial undecane in liquid=10 mol is probaby incorrect. I htink it should be 36.63/156, which means n2 is too high by about a factor of 40. We can also estimate initial values for ethylene in liquid and undecane in gas.

x1i = C1*VG;                % initial ethylene in gas (mol) 
x14i= 36.63/156;            % initial undecane in liquid (mol)
x7i = x14i*VG*K7/VL;        % initial undecane in gas (mol)
x8i = x1i*VL/(K1*VG);       % initial ethylene in liquid (mol)
x15i= Q1*C1;                % initial outflow rate (mol/s)
%...
xinit = [x1i; 0; 0; 0; 0; 0; x7i; x8i; 0; 0; 0; 0; 0; x14i; x15i];
[t,x]=ode45(dNdt,[0,18000],xinit);
dNdt=@(t,x) [Q1*C1-etc; Q2-etc; Q3-etc; etc];  % skipped a lot here    

x7i and x8i are chosen to acheive initial gas-liquid equilibrium if there is no inflow or outflow, or if there is equal inflow and outflow, and if there is no enzymatic consumption of ethylene. (There is never enzymatic production of ethylene, nor is there ever enzymatic consumption or production of undecane.)

x15i (initial outflow rate, mol/s) is chosen to produce net rate of change of moles=0, initially. (Woo et al. say on p.503 that the outflow is chosen to keep a constant number of total moles in the gas phase.) The only inflow is ethylene, at rate Q1*C1 mol/s, at all times.

I highly recommend that you check my equations for x1i, x14i, x7i, x8i, x15i, in case I made a mistake.

[William Rose]

@naiva saeedia, I have more comments and suggestions. Please email me securely by clicking on the envelope icon which appears when you click on the "WR" next to my Matlab answers posts, so we can continue the dicussion offline.

[naiva saeedia]

@William Rose

Hi. About your first question yeah I have changed initial conditions of undecane because x14i(36.63/156) gave me negative values for undecane in gas and liquid phases with this initial condition. I will check your equations for x1i,..,x15i and about n1 you are right since ethylene in the feed is all in gas phase. Sure. I will email you rn and thx again.

[naiva saeedia]

@William Rose

I have another question about time constant for outflow adjustment variable. we should assign something to this variable. right? How do we know the correct value for Tau_OF?