How to optimize the run time in my optimization problem.

3 views (last 30 days)
Hi guys. I have an optimization problem. The code is as follows:
Objective=@MassTransferErrors;
kL=0.5;
kH=0.5;
p0=[kL,kH];
A = [];
b = [];
Aeq = [];
beq = [];
lb=[0 ; 0];
ub=[100;100];
k = fmincon(Objective, p0, A, b, Aeq, beq, lb, ub);
disp(k)
function MTE=MassTransferErrors(p)
kL = p(1);
kH = p(2);
%% Constants
tMax=18000; % reaction duration (s)
Q0=[100 200 350 400 400 400 500]; % Q_G etylene inflow (ml/min)
T_C =[230 230 230 180 200 230 230]; %T for different cases
MTE_j=zeros(1,7);
Experiments = {[ 0.2985 0.6498 0.6147 0.43917 0.40398],[0.68662 1.6373 1.4260 1.4437 1.53169],[2.90493 5.68662 5.75704 2.65845 1.00352],[3.50352 11.3908 6.77817 3.46831 2.2007],[4.73592 10.8979 4.48944 3.01056 2.76408],[4.80634 9.45423 6.60211 4.03169 2.83451],[4.41901 10.4754 7.09507 4.13732 2.27113]};
for i=1:7
Q1=Q0(i)*1e-6/60; % Q_G ethylene inflow (m3/s)
Q2=0; % Q_G butene inflow
Q3=0; % Q_G hexene inflow
Q4=0; % Q_G octene inflow
Q5=0; % Q_G decene inflow
Q6=0; % Q_G dodecene inflow
Q7=0; % Q_G undecane inflow
P1=36e5; % ethylene inflow pressure [Pa]
T1=T_C(i)+273.15; % T_Ethylene [K]
T2=230+273.15; % T_ref [K]
R=8.314; % gas constant [J/(mol.K)]
C1=P1/(R*T1); % ethylene inlet gas concentration [mol/m^3]
VR=300e-6; % reactor volume [m^3]
VG=250e-6; % gas volume [m^3]
VL=50e-6; % liquid volume [m^3]
K=[3.24;2.23;1.72;0.2;0.1;0.08;0.09]; % solubility [nondim]
moleWt=[28;56;84;112;140;168;156]; % mole weight C2,C4,...,C12,C11 [g/mol]
wc=(0.3+0.25)*1e-3; % catalyst weight [kg]
kref=[2.224e-4;1.533e-4;3.988e-5;1.914e-7;4.328e-5;...
2.506e-7;4.036e-5;1.062e-6;6.055e-7;]; % rate at Tref=230C [mol/(s.g_cat)]
Eact=[109.5; 15.23; 7.88; 44.45; 9.438; 8.426; 10.91; 12.54; 7.127]; % activation energy [J/mol];
k=kref.*exp(-Eact*(1/T1-1/T2)/R); % rate at T=T2 [mol/(s.g)]
tauOF=5; % outflow time constant (s)
% Specify initial conditions
xinit=zeros(15,1); % initial state vector
xinit(1)=C1*VR; % initial ethylene in gas (mol)
xinit(14)=36.63/156; % initial undecane in liquid (mol)
xinit(7) = xinit(14)*VG*K(7)/VL; % initial undecane in gas (mol)
xinit(8) = xinit(1)*VL/(K(1)*VG); % initial ethylene in liquid (mol)
xinit(15)=Q1*C1; % initial outflow rate (mol/s)
nToti=sum(xinit(1:7)); % initial moles in gas (mol)
dNdt=@(t,x) [Q1*C1-x(15)*x(1)/sum(x(1:7))-VR*kL*(x(1)/VG-K(1)*x(8)/VL); % gas phase ethylene (mol/s)
Q2-x(15)*x(2)/sum(x(1:7))-VR*kL*(x(2)/VG-K(2)*x(9)/VL); % gas phase butene (mol/s)
Q3-x(15)*x(3)/sum(x(1:7))-VR*kL*(x(3)/VG-K(3)*x(10)/VL); % gas phase hexene (mol/s)
Q4-x(15)*x(4)/sum(x(1:7))-VR*kH*(x(4)/VG-K(4)*x(11)/VL); % gas phase octene (mol/s)
Q5-x(15)*x(5)/sum(x(1:7))-VR*kH*(x(5)/VG-K(5)*x(12)/VL); % gas phase decene (mol/s)
Q6-x(15)*x(6)/sum(x(1:7))-VR*kH*(x(6)/VG-K(6)*x(13)/VL); % gas phase dodecene (mol/s)
Q7-x(15)*x(7)/sum(x(1:7))-VR*kH*(x(7)/VG-K(7)*x(14)/VL); % gas phase undecane (mol/s)
VR*kL*(x(1)/VG-K(1)*x(8)/VL)+wc*(-2*k(1)*x(8)^2/VL^2-k(2)*x(8)*x(9)/VL^2-k(3)*x(8)*x(10)/VL^2-k(5)*x(8)*x(11)/VL^2-k(7)*x(8)*x(12)/VL^2);
VR*kL*(x(2)/VG-K(2)*x(9)/VL)+wc*(k(1)*x(8)^2/VL^2-k(2)*x(8)*x(9)/VL^2-2*k(4)*x(9)^2/VL.^2-k(6)*x(9)*x(10)/VL^2-k(8)*x(9)*x(11)/VL^2);
VR*kL*(x(3)/VG-K(3)*x(10)/VL)+wc*(k(2)*x(8)*x(9)/VL^2-k(3)*x(8)*x(10)/VL^2-k(6)*x(9)*x(10)/VL.^2-2*k(9)*x(10)^2/VL^2);
VR*kH*(x(4)/VG-K(4)*x(11)/VL)+wc*(k(3)*x(8)*x(10)/VL^2+k(4)*x(9)^2/VL^2-k(5)*x(8)*x(11)/VL^2-k(8)*x(9)*x(11)/VL^2);
VR*kH*(x(5)/VG-K(5)*x(12)/VL)+wc*(k(5)*x(8)*x(11)/VL^2+k(6)*x(9)*x(10)/VL^2-k(7)*x(8)*x(12)/VL^2);
VR*kH*(x(6)/VG-K(6)*x(13)/VL)+wc*(k(7)*x(8)*x(12)/VL^2+k(8)*x(9)*x(11)/VL^2+k(9)*x(10)^2/VL^2);
VR*kH*(x(7)/VG-K(7)*x(14)/VL);
(sum(x(1:7))-nToti)/tauOF]; % d(outflow rate)/dt (mol/s^2)
[t,x]=ode45(dNdt,[0,tMax],xinit);
molGend=x(end,1:7);
molLend=x(end,8:14);
massGend=molGend'.*moleWt;
massLend=molLend'.*moleWt;
%Total Product
TotalProduct = zeros(1,7);
for j=1:7
TotalProduct(j) = massGend(j) + massLend(j); %Sum of the liquid and gas phase products(g)
end
Experiment_i = cell2mat(Experiments(i)); %Converting Experiments set to matrix
MTE_i = (TotalProduct(2)-Experiment_i(1))^2+(TotalProduct(3)-Experiment_i(2))^2+(TotalProduct(4)-Experiment_i(3))^2+(TotalProduct(5)-Experiment_i(4))^2+(TotalProduct(6)-Experiment_i(5))^5; %Defining an Mass Transfer Error relation
MTE_j(i) = MTE_i; %Defines a Mass Transfer Error Vector(1*7) that contains the error for each case
end
MTE = sum(MTE_j(1:7)); %Objective function(Sum of the all arrays in MTE_j Vector) what I need to minimize is each array that is on the MTE_j Vector but since I can't return a Vector as an objective function I sum all the arrays as my objective function.
end
What I need to minimize in this problem are all the seven values in MTE_j vector but since objective function can't return a vector. I have used the sum of the all values. I doubt this is the correct way to minimize all the seven values in MTE_j vector also my code took an extremely long run time(6 hours last time I have checked with no answers yet) . I know the objective function is complicated but I guess I'm doing something wrong. I also test my objective function with an assumption for kL and kH values. my objective function code seems to work OK. Here's the test code:
%%Test Objective function
kL = 0.1; %Assumption for kL
kH = 0.1; %Assumption for kH
%% Constants
tMax=18000; % reaction duration (s)
Q0=[100 200 350 400 400 400 500]; % Q_G etylene inflow (ml/min)
T_C =[230 230 230 180 200 230 230]; %T for different cases
MTE_j=zeros(1,7);
Experiments = {[ 0.2985 0.6498 0.6147 0.43917 0.40398],[0.68662 1.6373 1.4260 1.4437 1.53169],[2.90493 5.68662 5.75704 2.65845 1.00352],[3.50352 11.3908 6.77817 3.46831 2.2007],[4.73592 10.8979 4.48944 3.01056 2.76408],[4.80634 9.45423 6.60211 4.03169 2.83451],[4.41901 10.4754 7.09507 4.13732 2.27113]};
for i=1:7
Q1=Q0(i)*1e-6/60; % Q_G ethylene inflow (m3/s)
Q2=0; % Q_G butene inflow
Q3=0; % Q_G hexene inflow
Q4=0; % Q_G octene inflow
Q5=0; % Q_G decene inflow
Q6=0; % Q_G dodecene inflow
Q7=0; % Q_G undecane inflow
P1=36e5; % ethylene inflow pressure [Pa]
T1=T_C(i)+273.15; % T_Ethylene [K]
T2=230+273.15; % T_ref [K]
R=8.314; % gas constant [J/(mol.K)]
C1=P1/(R*T1); % ethylene inlet gas concentration [mol/m^3]
VR=300e-6; % reactor volume [m^3]
VG=250e-6; % gas volume [m^3]
VL=50e-6; % liquid volume [m^3]
K=[3.24;2.23;1.72;0.2;0.1;0.08;0.09]; % solubility [nondim]
moleWt=[28;56;84;112;140;168;156]; % mole weight C2,C4,...,C12,C11 [g/mol]
wc=(0.3+0.25)*1e-3; % catalyst weight [kg]
kref=[2.224e-4;1.533e-4;3.988e-5;1.914e-7;4.328e-5;...
2.506e-7;4.036e-5;1.062e-6;6.055e-7;]; % rate at Tref=230C [mol/(s.g_cat)]
Eact=[109.5; 15.23; 7.88; 44.45; 9.438; 8.426; 10.91; 12.54; 7.127]; % activation energy [J/mol];
k=kref.*exp(-Eact*(1/T1-1/T2)/R); % rate at T=T2 [mol/(s.g)]
tauOF=5; % outflow time constant (s)
% Specify initial conditions
xinit=zeros(15,1); % initial state vector
xinit(1)=C1*VR; % initial ethylene in gas (mol)
xinit(14)=36.63/156; % initial undecane in liquid (mol)
xinit(7) = xinit(14)*VG*K(7)/VL; % initial undecane in gas (mol)
xinit(8) = xinit(1)*VL/(K(1)*VG); % initial ethylene in liquid (mol)
xinit(15)=Q1*C1; % initial outflow rate (mol/s)
nToti=sum(xinit(1:7)); % initial moles in gas (mol)
dNdt=@(t,x) [Q1*C1-x(15)*x(1)/sum(x(1:7))-VR*kL*(x(1)/VG-K(1)*x(8)/VL); % gas phase ethylene (mol/s)
Q2-x(15)*x(2)/sum(x(1:7))-VR*kL*(x(2)/VG-K(2)*x(9)/VL); % gas phase butene (mol/s)
Q3-x(15)*x(3)/sum(x(1:7))-VR*kL*(x(3)/VG-K(3)*x(10)/VL); % gas phase hexene (mol/s)
Q4-x(15)*x(4)/sum(x(1:7))-VR*kH*(x(4)/VG-K(4)*x(11)/VL); % gas phase octene (mol/s)
Q5-x(15)*x(5)/sum(x(1:7))-VR*kH*(x(5)/VG-K(5)*x(12)/VL); % gas phase decene (mol/s)
Q6-x(15)*x(6)/sum(x(1:7))-VR*kH*(x(6)/VG-K(6)*x(13)/VL); % gas phase dodecene (mol/s)
Q7-x(15)*x(7)/sum(x(1:7))-VR*kH*(x(7)/VG-K(7)*x(14)/VL); % gas phase undecane (mol/s)
VR*kL*(x(1)/VG-K(1)*x(8)/VL)+wc*(-2*k(1)*x(8)^2/VL^2-k(2)*x(8)*x(9)/VL^2-k(3)*x(8)*x(10)/VL^2-k(5)*x(8)*x(11)/VL^2-k(7)*x(8)*x(12)/VL^2);
VR*kL*(x(2)/VG-K(2)*x(9)/VL)+wc*(k(1)*x(8)^2/VL^2-k(2)*x(8)*x(9)/VL^2-2*k(4)*x(9)^2/VL.^2-k(6)*x(9)*x(10)/VL^2-k(8)*x(9)*x(11)/VL^2);
VR*kL*(x(3)/VG-K(3)*x(10)/VL)+wc*(k(2)*x(8)*x(9)/VL^2-k(3)*x(8)*x(10)/VL^2-k(6)*x(9)*x(10)/VL.^2-2*k(9)*x(10)^2/VL^2);
VR*kH*(x(4)/VG-K(4)*x(11)/VL)+wc*(k(3)*x(8)*x(10)/VL^2+k(4)*x(9)^2/VL^2-k(5)*x(8)*x(11)/VL^2-k(8)*x(9)*x(11)/VL^2);
VR*kH*(x(5)/VG-K(5)*x(12)/VL)+wc*(k(5)*x(8)*x(11)/VL^2+k(6)*x(9)*x(10)/VL^2-k(7)*x(8)*x(12)/VL^2);
VR*kH*(x(6)/VG-K(6)*x(13)/VL)+wc*(k(7)*x(8)*x(12)/VL^2+k(8)*x(9)*x(11)/VL^2+k(9)*x(10)^2/VL^2);
VR*kH*(x(7)/VG-K(7)*x(14)/VL);
(sum(x(1:7))-nToti)/tauOF]; % d(outflow rate)/dt (mol/s^2)
[t,x]=ode45(dNdt,[0,tMax],xinit);
molGend=x(end,1:7);
molLend=x(end,8:14);
massGend=molGend'.*moleWt;
massLend=molLend'.*moleWt;
%Total Product
TotalProduct = zeros(1,7);
for j=1:7
TotalProduct(j) = massGend(j) + massLend(j); %Sum of the liquid and gas phase products(g)
end
Experiment_i = cell2mat(Experiments(i));
MTE_i = (TotalProduct(2)-Experiment_i(1))^2+(TotalProduct(3)-Experiment_i(2))^2+(TotalProduct(4)-Experiment_i(3))^2+(TotalProduct(5)-Experiment_i(4))^2+(TotalProduct(6)-Experiment_i(5))^2;
MTE_j(i) = MTE_i;
end
MTE = sum(MTE_j(1:7));
I just want the kL and kH values that minimize each value on MTE_j vector and I want a code to actually give me these values. My code take an extremely long run time with no answers.

Accepted Answer

Torsten
Torsten on 28 Apr 2024
I don't see anything obviously wrong in your coding - except for the last term in MTE_i which should be
(TotalProduct(6)-Experiment_i(5))^2
instead of
(TotalProduct(6)-Experiment_i(5))^5
Try whether switching to ode15s from ode45 will reduce the runtime.
  12 Comments
Torsten
Torsten on 29 Apr 2024
Load it into the editor and click the "Run" button under "Octave".
But now with 1 min runtime, I'd prefer MATLAB ...
function main
pkg load "optim"
Objective=@MassTransferErrors;
kL=0.5;
kH=0.5;
p0=[kL,kH];
A = [];
b = [];
Aeq = [];
beq = [];
lb=[0 ; 0];
ub=[100;100];
k = fmincon(Objective, p0, A, b, Aeq, beq, lb, ub);
disp(k)
end
function MTE=MassTransferErrors(p)
kL = p(1)
kH = p(2)
%% Constants
tMax=18000; % reaction duration (s)
Q0=[100 200 350 400 400 400 500]; % Q_G etylene inflow (ml/min)
T_C =[230 230 230 180 200 230 230]; %T for different cases
MTE_j=zeros(1,7);
Experiments = {[ 0.2985 0.6498 0.6147 0.43917 0.40398],[0.68662 1.6373 1.4260 1.4437 1.53169],[2.90493 5.68662 5.75704 2.65845 1.00352],[3.50352 11.3908 6.77817 3.46831 2.2007],[4.73592 10.8979 4.48944 3.01056 2.76408],[4.80634 9.45423 6.60211 4.03169 2.83451],[4.41901 10.4754 7.09507 4.13732 2.27113]};
for i=1:7
Q1=Q0(i)*1e-6/60; % Q_G ethylene inflow (m3/s)
Q2=0; % Q_G butene inflow
Q3=0; % Q_G hexene inflow
Q4=0; % Q_G octene inflow
Q5=0; % Q_G decene inflow
Q6=0; % Q_G dodecene inflow
Q7=0; % Q_G undecane inflow
P1=36e5; % ethylene inflow pressure [Pa]
T1=T_C(i)+273.15; % T_Ethylene [K]
T2=230+273.15; % T_ref [K]
R=8.314; % gas constant [J/(mol.K)]
C1=P1/(R*T1); % ethylene inlet gas concentration [mol/m^3]
VR=300e-6; % reactor volume [m^3]
VG=250e-6; % gas volume [m^3]
VL=50e-6; % liquid volume [m^3]
K=[3.24;2.23;1.72;0.2;0.1;0.08;0.09]; % solubility [nondim]
moleWt=[28;56;84;112;140;168;156]; % mole weight C2,C4,...,C12,C11 [g/mol]
wc=(0.3+0.25)*1e-3; % catalyst weight [kg]
kref=[2.224e-4;1.533e-4;3.988e-5;1.914e-7;4.328e-5;...
2.506e-7;4.036e-5;1.062e-6;6.055e-7;]; % rate at Tref=230C [mol/(s.g_cat)]
Eact=[109.5; 15.23; 7.88; 44.45; 9.438; 8.426; 10.91; 12.54; 7.127]; % activation energy [J/mol];
k=kref.*exp(-Eact*(1/T1-1/T2)/R); % rate at T=T2 [mol/(s.g)]
tauOF=5; % outflow time constant (s)
% Specify initial conditions
xinit=zeros(15,1); % initial state vector
xinit(1)=C1*VR; % initial ethylene in gas (mol)
xinit(14)=36.63/156; % initial undecane in liquid (mol)
xinit(7) = xinit(14)*VG*K(7)/VL; % initial undecane in gas (mol)
xinit(8) = xinit(1)*VL/(K(1)*VG); % initial ethylene in liquid (mol)
xinit(15)=Q1*C1; % initial outflow rate (mol/s)
nToti=sum(xinit(1:7)); % initial moles in gas (mol)
dNdt=@(t,x) [Q1*C1-x(15)*x(1)/sum(x(1:7))-VR*kL*(x(1)/VG-K(1)*x(8)/VL); % gas phase ethylene (mol/s)
Q2-x(15)*x(2)/sum(x(1:7))-VR*kL*(x(2)/VG-K(2)*x(9)/VL); % gas phase butene (mol/s)
Q3-x(15)*x(3)/sum(x(1:7))-VR*kL*(x(3)/VG-K(3)*x(10)/VL); % gas phase hexene (mol/s)
Q4-x(15)*x(4)/sum(x(1:7))-VR*kH*(x(4)/VG-K(4)*x(11)/VL); % gas phase octene (mol/s)
Q5-x(15)*x(5)/sum(x(1:7))-VR*kH*(x(5)/VG-K(5)*x(12)/VL); % gas phase decene (mol/s)
Q6-x(15)*x(6)/sum(x(1:7))-VR*kH*(x(6)/VG-K(6)*x(13)/VL); % gas phase dodecene (mol/s)
Q7-x(15)*x(7)/sum(x(1:7))-VR*kH*(x(7)/VG-K(7)*x(14)/VL); % gas phase undecane (mol/s)
VR*kL*(x(1)/VG-K(1)*x(8)/VL)+wc*(-2*k(1)*x(8)^2/VL^2-k(2)*x(8)*x(9)/VL^2-k(3)*x(8)*x(10)/VL^2-k(5)*x(8)*x(11)/VL^2-k(7)*x(8)*x(12)/VL^2);
VR*kL*(x(2)/VG-K(2)*x(9)/VL)+wc*(k(1)*x(8)^2/VL^2-k(2)*x(8)*x(9)/VL^2-2*k(4)*x(9)^2/VL.^2-k(6)*x(9)*x(10)/VL^2-k(8)*x(9)*x(11)/VL^2);
VR*kL*(x(3)/VG-K(3)*x(10)/VL)+wc*(k(2)*x(8)*x(9)/VL^2-k(3)*x(8)*x(10)/VL^2-k(6)*x(9)*x(10)/VL.^2-2*k(9)*x(10)^2/VL^2);
VR*kH*(x(4)/VG-K(4)*x(11)/VL)+wc*(k(3)*x(8)*x(10)/VL^2+k(4)*x(9)^2/VL^2-k(5)*x(8)*x(11)/VL^2-k(8)*x(9)*x(11)/VL^2);
VR*kH*(x(5)/VG-K(5)*x(12)/VL)+wc*(k(5)*x(8)*x(11)/VL^2+k(6)*x(9)*x(10)/VL^2-k(7)*x(8)*x(12)/VL^2);
VR*kH*(x(6)/VG-K(6)*x(13)/VL)+wc*(k(7)*x(8)*x(12)/VL^2+k(8)*x(9)*x(11)/VL^2+k(9)*x(10)^2/VL^2);
VR*kH*(x(7)/VG-K(7)*x(14)/VL);
(sum(x(1:7))-nToti)/tauOF]; % d(outflow rate)/dt (mol/s^2)
[t,x]=ode15s(dNdt,[0,tMax],xinit);
molGend=x(end,1:7);
molLend=x(end,8:14);
massGend=molGend'.*moleWt;
massLend=molLend'.*moleWt;
%Total Product
TotalProduct = zeros(1,7);
for j=1:7
TotalProduct(j) = massGend(j) + massLend(j); %Sum of the liquid and gas phase products(g)
end
Experiment_i = cell2mat(Experiments(i)); %Converting Experiments set to matrix
MTE_i = (TotalProduct(2)-Experiment_i(1))^2+(TotalProduct(3)-Experiment_i(2))^2+(TotalProduct(4)-Experiment_i(3))^2+(TotalProduct(5)-Experiment_i(4))^2+(TotalProduct(6)-Experiment_i(5))^5; %Defining an Mass Transfer Error relation
MTE_j(i) = MTE_i; %Defines a Mass Transfer Error Vector(1*7) that contains the error for each case
end
MTE = sum(MTE_j(1:7)) %Objective function(Sum of the all arrays in MTE_j Vector) what I need to minimize is each array that is on the MTE_j Vector but since I can't return a Vector as an objective function I sum all the arrays as my objective function.
end

Sign in to comment.

More Answers (0)

Categories

Find more on Elements in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!