Improving optimization results(Fmincon)
Hi guys. I’m trying to solve an optimization problem but the results I’m getting from fmincon() don’t have the accuracy that I’m looking for. I have tried to change the alghorithm and step tolerance but they didn’t affect the results. Is there any way to improve fmincon() results?? Here is my optimization problem code:
clear variables
clc
Objective=@MassTransferErrors_Closed_loop;
A = [];
b = [];
Aeq = [];
beq = [];
lb=[0 ; 0];
ub=[1000;1000];
kL=500;
kH=500;
p0=[kL,kH];
options = optimoptions(‘fmincon’,’Display’,’iter’,’Algorithm’,’sqp-legacy’,’StepTolerance’,1e-11,"MaxFunctionEvaluations",2e3);
nlcon =[];
k = fmincon(Objective, p0, A, b, Aeq, beq, lb, ub, nlcon, options);
disp(k)
function MTE=MassTransferErrors_Closed_loop(p)
kL=p(1);
kH=p(2);
moleWt=[28;56;84;112;140;168;156]; % mole weight C2,C4,…,C12,C11 [g/mol]
%Initial Condition
Q0=[100 200 350 400 400 400 500]; % Q_G etylene inflow (ml/min)
T1=[230 230 230 180 200 230 230]; %T for different cases;
Kc_Total_gases=1;
tauI_Total_gases=1;
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
y0 = [0.258176232100050 0 0 0 0 0 0.105663461538462 0.0159368044506204 0 0 0 0 0 0.234807692307692 0];
%options = odeset(‘RelTol’,1e-5,’AbsTol’,1e-7);
[t,y]= ode23s(@(t,y) Sec_model_fun_for_optimization(t,y,Q0(i),T1(i),kL,kH,Kc_Total_gases,tauI_Total_gases),[0 18000],y0);
a = zeros(numel(t),1);
for q=1:numel(t)
a(q) = sum(y(q,1:7));
end
% Create plots for the gas mols in the reactor
% Total Gas mol in the reactor %[=mol]
%subplot(2,1,1)
%figure(i)
%plot(t,a,’r’,’LineWidth’,1.5)
%legend(‘Total Gas mols’)
%xlabel(‘Time [s]’)
%ylabel(‘n [mol]’)
%title(‘Total Gas mols’)
% for t=100s
%subplot(2,1,2)
%t_new=1:50; %time vector for interpolation for plotting for the first 300 seconds
%aint =interp1(t,a,t_new,’pchip’); %interpolated state matrix for plotting
%plot(t_new,aint,’b’,’LineWidth’,1.5)
%legend(‘Total Gas mols’)
%xlabel(‘Time [s] (first 50s)’)
%ylabel(‘n [mol]’)
%title(‘Total Gas mols’)
%hold off
% moles of products in gas and liquid at end
molGend=y(end,1:7);
molLend=y(end,8:14);
% product masses [g] in gas and liquid at end
massGend=molGend’.*moleWt;
massLend=molLend’.*moleWt;
%Total mass
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))/Experiment_i(1))^2+((TotalProduct(3)-Experiment_i(2))/Experiment_i(2))^2+((TotalProduct(4)-Experiment_i(3))/Experiment_i(3))^2+((TotalProduct(5)-Experiment_i(4))/Experiment_i(4))^2+((TotalProduct(6)-Experiment_i(5))/Experiment_i(5))^2; %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
Also here is my code for the function that I have used in ode23s:
function S = Sec_model_fun_for_optimization(t,x,Q0,T1,kL,kH,Kc_Total_gases,tauI_Total_gases)
%%Dynamic state inputs
n_C2_gas=x(1); %Ethylene mols in gas phase[mol]
n_C4_gas=x(2); %Butene mols in gas phase[mol]
n_C6_gas=x(3); %Hexene mols in gas phase[mol]
n_C8_gas=x(4); %Octene mols in gas phase[mol]
n_C10_gas=x(5); %Decene mols in gas phase[mol]
n_C12_gas=x(6); %Dodecene mols in gas phase[mol]
n_C11_gas=x(7); %Undecene mols in gas phase[mol]
n_C2_liquid=x(8); %Ethylene mols in liquid phase[mol]
n_C4_liquid=x(9); %Butene mols in liquid phase[mol]
n_C6_liquid=x(10); %Hexene mols in liquid phase[mol]
n_C8_liquid=x(11); %Octene mols in liquid phase[mol]
n_C10_liquid=x(12); %Decene mols in liquid phase[mol]
n_C12_liquid=x(13); %Dodecene mols in liquid phase[mol]
n_C11_liquid=x(14); %Undecene mols in liquid phase[mol]
e_integral = x(15); %Error[mol/s]
%% Constants
%Q0=350; % Q_G etylene inflow (ml/min)
Q1=Q0*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=230+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]
F0=0.0179; % gas outflow rate [mmol/s]
F1=F0.*1e-3; % gas outflow rate [mol/s]
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)]
% 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)
%%Setpoint
nGin_setpoint=0.363839693638512;
%Set Point Tracking & Load Rejection
%t1 = 1800; t2 = 3600; t3 = 5400; t4 = 7200; t5 = 9000; t6 = 10800; t7 = 12600; t8 = 14400; t9 = 16200; t10 = 18000;
% Step 1
% +5% change in set point
%if t >= t1 && t <= t2
% nGin_setpoint = 0.616018502797380*1.05;
% Step 2
% -5% change in set point and disturbance
%elseif t >= t3 && t <= t4
% nGin_setpoint = 0.616018502797380*0.95;
% Step 3
% +10% change in setpoint
%elseif t >=t5 && t<= t6
% nGin_setpoint = 0.616018502797380*1.1;
% Step 4
% -10% change in set point
%elseif t >=t7 && t <= t8
% nGin_setpoint = 0.616018502797380*0.9;
% Step 5
% +20% change in set point
%elseif t >= t9 && t <= t10
% nGin_setpoint = 0.616018502797380*1.2;
%end
e_mol_gases = sum(x(1:7)) – nGin_setpoint;
F_G_R = Kc_Total_gases*(e_mol_gases+tauI_Total_gases*e_integral);
%Right-hand side evaluation of the dynamic model (DAE set)
S1 = Q1*C1-F_G_R*x(1)/sum(x(1:7))-VR*kL*(x(1)/VG-K(1)*x(8)/VL); % gas phase ethylene (mol/s)
S2 = Q2-F_G_R*x(2)/sum(x(1:7))-VR*kL*(x(2)/VG-K(2)*x(9)/VL); % gas phase butene (mol/s);
S3 = Q3-F_G_R*x(3)/sum(x(1:7))-VR*kL*(x(3)/VG-K(3)*x(10)/VL); % gas phase hexene (mol/s);
S4= Q4-F_G_R*x(4)/sum(x(1:7))-VR*kH*(x(4)/VG-K(4)*x(11)/VL); % gas phase octene (mol/s);
S5= Q5-F_G_R*x(5)/sum(x(1:7))-VR*kH*(x(5)/VG-K(5)*x(12)/VL); % gas phase decene (mol/s);
S6= Q6-F_G_R*x(6)/sum(x(1:7))-VR*kH*(x(6)/VG-K(6)*x(13)/VL); % gas phase dodecene (mol/s);
S7= Q7-F_G_R*x(7)/sum(x(1:7))-VR*kH*(x(7)/VG-K(7)*x(14)/VL); % gas phase undecane (mol/s) ;
S8= 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);
S9= 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);
S10= 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);
S11= 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);
S12= 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);
S13= 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);
S14= VR*kH*(x(7)/VG-K(7)*x(14)/VL);
S15= sum(x(1:7))-(nGin_setpoint); %Error
S = ([S1; S2; S3; S4; S5; S6; S7; S8; S9; S10; S11; S12; S13; S14; S15]);
end
Also the results are changing dramatically when I change initial values for kL and kH. I know it’s normal since fmincon() doesn’t compute global maximum/minimum but it’s just so weird to get different resluts whenever I change kL and kH values!Hi guys. I’m trying to solve an optimization problem but the results I’m getting from fmincon() don’t have the accuracy that I’m looking for. I have tried to change the alghorithm and step tolerance but they didn’t affect the results. Is there any way to improve fmincon() results?? Here is my optimization problem code:
clear variables
clc
Objective=@MassTransferErrors_Closed_loop;
A = [];
b = [];
Aeq = [];
beq = [];
lb=[0 ; 0];
ub=[1000;1000];
kL=500;
kH=500;
p0=[kL,kH];
options = optimoptions(‘fmincon’,’Display’,’iter’,’Algorithm’,’sqp-legacy’,’StepTolerance’,1e-11,"MaxFunctionEvaluations",2e3);
nlcon =[];
k = fmincon(Objective, p0, A, b, Aeq, beq, lb, ub, nlcon, options);
disp(k)
function MTE=MassTransferErrors_Closed_loop(p)
kL=p(1);
kH=p(2);
moleWt=[28;56;84;112;140;168;156]; % mole weight C2,C4,…,C12,C11 [g/mol]
%Initial Condition
Q0=[100 200 350 400 400 400 500]; % Q_G etylene inflow (ml/min)
T1=[230 230 230 180 200 230 230]; %T for different cases;
Kc_Total_gases=1;
tauI_Total_gases=1;
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
y0 = [0.258176232100050 0 0 0 0 0 0.105663461538462 0.0159368044506204 0 0 0 0 0 0.234807692307692 0];
%options = odeset(‘RelTol’,1e-5,’AbsTol’,1e-7);
[t,y]= ode23s(@(t,y) Sec_model_fun_for_optimization(t,y,Q0(i),T1(i),kL,kH,Kc_Total_gases,tauI_Total_gases),[0 18000],y0);
a = zeros(numel(t),1);
for q=1:numel(t)
a(q) = sum(y(q,1:7));
end
% Create plots for the gas mols in the reactor
% Total Gas mol in the reactor %[=mol]
%subplot(2,1,1)
%figure(i)
%plot(t,a,’r’,’LineWidth’,1.5)
%legend(‘Total Gas mols’)
%xlabel(‘Time [s]’)
%ylabel(‘n [mol]’)
%title(‘Total Gas mols’)
% for t=100s
%subplot(2,1,2)
%t_new=1:50; %time vector for interpolation for plotting for the first 300 seconds
%aint =interp1(t,a,t_new,’pchip’); %interpolated state matrix for plotting
%plot(t_new,aint,’b’,’LineWidth’,1.5)
%legend(‘Total Gas mols’)
%xlabel(‘Time [s] (first 50s)’)
%ylabel(‘n [mol]’)
%title(‘Total Gas mols’)
%hold off
% moles of products in gas and liquid at end
molGend=y(end,1:7);
molLend=y(end,8:14);
% product masses [g] in gas and liquid at end
massGend=molGend’.*moleWt;
massLend=molLend’.*moleWt;
%Total mass
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))/Experiment_i(1))^2+((TotalProduct(3)-Experiment_i(2))/Experiment_i(2))^2+((TotalProduct(4)-Experiment_i(3))/Experiment_i(3))^2+((TotalProduct(5)-Experiment_i(4))/Experiment_i(4))^2+((TotalProduct(6)-Experiment_i(5))/Experiment_i(5))^2; %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
Also here is my code for the function that I have used in ode23s:
function S = Sec_model_fun_for_optimization(t,x,Q0,T1,kL,kH,Kc_Total_gases,tauI_Total_gases)
%%Dynamic state inputs
n_C2_gas=x(1); %Ethylene mols in gas phase[mol]
n_C4_gas=x(2); %Butene mols in gas phase[mol]
n_C6_gas=x(3); %Hexene mols in gas phase[mol]
n_C8_gas=x(4); %Octene mols in gas phase[mol]
n_C10_gas=x(5); %Decene mols in gas phase[mol]
n_C12_gas=x(6); %Dodecene mols in gas phase[mol]
n_C11_gas=x(7); %Undecene mols in gas phase[mol]
n_C2_liquid=x(8); %Ethylene mols in liquid phase[mol]
n_C4_liquid=x(9); %Butene mols in liquid phase[mol]
n_C6_liquid=x(10); %Hexene mols in liquid phase[mol]
n_C8_liquid=x(11); %Octene mols in liquid phase[mol]
n_C10_liquid=x(12); %Decene mols in liquid phase[mol]
n_C12_liquid=x(13); %Dodecene mols in liquid phase[mol]
n_C11_liquid=x(14); %Undecene mols in liquid phase[mol]
e_integral = x(15); %Error[mol/s]
%% Constants
%Q0=350; % Q_G etylene inflow (ml/min)
Q1=Q0*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=230+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]
F0=0.0179; % gas outflow rate [mmol/s]
F1=F0.*1e-3; % gas outflow rate [mol/s]
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)]
% 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)
%%Setpoint
nGin_setpoint=0.363839693638512;
%Set Point Tracking & Load Rejection
%t1 = 1800; t2 = 3600; t3 = 5400; t4 = 7200; t5 = 9000; t6 = 10800; t7 = 12600; t8 = 14400; t9 = 16200; t10 = 18000;
% Step 1
% +5% change in set point
%if t >= t1 && t <= t2
% nGin_setpoint = 0.616018502797380*1.05;
% Step 2
% -5% change in set point and disturbance
%elseif t >= t3 && t <= t4
% nGin_setpoint = 0.616018502797380*0.95;
% Step 3
% +10% change in setpoint
%elseif t >=t5 && t<= t6
% nGin_setpoint = 0.616018502797380*1.1;
% Step 4
% -10% change in set point
%elseif t >=t7 && t <= t8
% nGin_setpoint = 0.616018502797380*0.9;
% Step 5
% +20% change in set point
%elseif t >= t9 && t <= t10
% nGin_setpoint = 0.616018502797380*1.2;
%end
e_mol_gases = sum(x(1:7)) – nGin_setpoint;
F_G_R = Kc_Total_gases*(e_mol_gases+tauI_Total_gases*e_integral);
%Right-hand side evaluation of the dynamic model (DAE set)
S1 = Q1*C1-F_G_R*x(1)/sum(x(1:7))-VR*kL*(x(1)/VG-K(1)*x(8)/VL); % gas phase ethylene (mol/s)
S2 = Q2-F_G_R*x(2)/sum(x(1:7))-VR*kL*(x(2)/VG-K(2)*x(9)/VL); % gas phase butene (mol/s);
S3 = Q3-F_G_R*x(3)/sum(x(1:7))-VR*kL*(x(3)/VG-K(3)*x(10)/VL); % gas phase hexene (mol/s);
S4= Q4-F_G_R*x(4)/sum(x(1:7))-VR*kH*(x(4)/VG-K(4)*x(11)/VL); % gas phase octene (mol/s);
S5= Q5-F_G_R*x(5)/sum(x(1:7))-VR*kH*(x(5)/VG-K(5)*x(12)/VL); % gas phase decene (mol/s);
S6= Q6-F_G_R*x(6)/sum(x(1:7))-VR*kH*(x(6)/VG-K(6)*x(13)/VL); % gas phase dodecene (mol/s);
S7= Q7-F_G_R*x(7)/sum(x(1:7))-VR*kH*(x(7)/VG-K(7)*x(14)/VL); % gas phase undecane (mol/s) ;
S8= 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);
S9= 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);
S10= 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);
S11= 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);
S12= 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);
S13= 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);
S14= VR*kH*(x(7)/VG-K(7)*x(14)/VL);
S15= sum(x(1:7))-(nGin_setpoint); %Error
S = ([S1; S2; S3; S4; S5; S6; S7; S8; S9; S10; S11; S12; S13; S14; S15]);
end
Also the results are changing dramatically when I change initial values for kL and kH. I know it’s normal since fmincon() doesn’t compute global maximum/minimum but it’s just so weird to get different resluts whenever I change kL and kH values! Hi guys. I’m trying to solve an optimization problem but the results I’m getting from fmincon() don’t have the accuracy that I’m looking for. I have tried to change the alghorithm and step tolerance but they didn’t affect the results. Is there any way to improve fmincon() results?? Here is my optimization problem code:
clear variables
clc
Objective=@MassTransferErrors_Closed_loop;
A = [];
b = [];
Aeq = [];
beq = [];
lb=[0 ; 0];
ub=[1000;1000];
kL=500;
kH=500;
p0=[kL,kH];
options = optimoptions(‘fmincon’,’Display’,’iter’,’Algorithm’,’sqp-legacy’,’StepTolerance’,1e-11,"MaxFunctionEvaluations",2e3);
nlcon =[];
k = fmincon(Objective, p0, A, b, Aeq, beq, lb, ub, nlcon, options);
disp(k)
function MTE=MassTransferErrors_Closed_loop(p)
kL=p(1);
kH=p(2);
moleWt=[28;56;84;112;140;168;156]; % mole weight C2,C4,…,C12,C11 [g/mol]
%Initial Condition
Q0=[100 200 350 400 400 400 500]; % Q_G etylene inflow (ml/min)
T1=[230 230 230 180 200 230 230]; %T for different cases;
Kc_Total_gases=1;
tauI_Total_gases=1;
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
y0 = [0.258176232100050 0 0 0 0 0 0.105663461538462 0.0159368044506204 0 0 0 0 0 0.234807692307692 0];
%options = odeset(‘RelTol’,1e-5,’AbsTol’,1e-7);
[t,y]= ode23s(@(t,y) Sec_model_fun_for_optimization(t,y,Q0(i),T1(i),kL,kH,Kc_Total_gases,tauI_Total_gases),[0 18000],y0);
a = zeros(numel(t),1);
for q=1:numel(t)
a(q) = sum(y(q,1:7));
end
% Create plots for the gas mols in the reactor
% Total Gas mol in the reactor %[=mol]
%subplot(2,1,1)
%figure(i)
%plot(t,a,’r’,’LineWidth’,1.5)
%legend(‘Total Gas mols’)
%xlabel(‘Time [s]’)
%ylabel(‘n [mol]’)
%title(‘Total Gas mols’)
% for t=100s
%subplot(2,1,2)
%t_new=1:50; %time vector for interpolation for plotting for the first 300 seconds
%aint =interp1(t,a,t_new,’pchip’); %interpolated state matrix for plotting
%plot(t_new,aint,’b’,’LineWidth’,1.5)
%legend(‘Total Gas mols’)
%xlabel(‘Time [s] (first 50s)’)
%ylabel(‘n [mol]’)
%title(‘Total Gas mols’)
%hold off
% moles of products in gas and liquid at end
molGend=y(end,1:7);
molLend=y(end,8:14);
% product masses [g] in gas and liquid at end
massGend=molGend’.*moleWt;
massLend=molLend’.*moleWt;
%Total mass
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))/Experiment_i(1))^2+((TotalProduct(3)-Experiment_i(2))/Experiment_i(2))^2+((TotalProduct(4)-Experiment_i(3))/Experiment_i(3))^2+((TotalProduct(5)-Experiment_i(4))/Experiment_i(4))^2+((TotalProduct(6)-Experiment_i(5))/Experiment_i(5))^2; %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
Also here is my code for the function that I have used in ode23s:
function S = Sec_model_fun_for_optimization(t,x,Q0,T1,kL,kH,Kc_Total_gases,tauI_Total_gases)
%%Dynamic state inputs
n_C2_gas=x(1); %Ethylene mols in gas phase[mol]
n_C4_gas=x(2); %Butene mols in gas phase[mol]
n_C6_gas=x(3); %Hexene mols in gas phase[mol]
n_C8_gas=x(4); %Octene mols in gas phase[mol]
n_C10_gas=x(5); %Decene mols in gas phase[mol]
n_C12_gas=x(6); %Dodecene mols in gas phase[mol]
n_C11_gas=x(7); %Undecene mols in gas phase[mol]
n_C2_liquid=x(8); %Ethylene mols in liquid phase[mol]
n_C4_liquid=x(9); %Butene mols in liquid phase[mol]
n_C6_liquid=x(10); %Hexene mols in liquid phase[mol]
n_C8_liquid=x(11); %Octene mols in liquid phase[mol]
n_C10_liquid=x(12); %Decene mols in liquid phase[mol]
n_C12_liquid=x(13); %Dodecene mols in liquid phase[mol]
n_C11_liquid=x(14); %Undecene mols in liquid phase[mol]
e_integral = x(15); %Error[mol/s]
%% Constants
%Q0=350; % Q_G etylene inflow (ml/min)
Q1=Q0*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=230+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]
F0=0.0179; % gas outflow rate [mmol/s]
F1=F0.*1e-3; % gas outflow rate [mol/s]
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)]
% 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)
%%Setpoint
nGin_setpoint=0.363839693638512;
%Set Point Tracking & Load Rejection
%t1 = 1800; t2 = 3600; t3 = 5400; t4 = 7200; t5 = 9000; t6 = 10800; t7 = 12600; t8 = 14400; t9 = 16200; t10 = 18000;
% Step 1
% +5% change in set point
%if t >= t1 && t <= t2
% nGin_setpoint = 0.616018502797380*1.05;
% Step 2
% -5% change in set point and disturbance
%elseif t >= t3 && t <= t4
% nGin_setpoint = 0.616018502797380*0.95;
% Step 3
% +10% change in setpoint
%elseif t >=t5 && t<= t6
% nGin_setpoint = 0.616018502797380*1.1;
% Step 4
% -10% change in set point
%elseif t >=t7 && t <= t8
% nGin_setpoint = 0.616018502797380*0.9;
% Step 5
% +20% change in set point
%elseif t >= t9 && t <= t10
% nGin_setpoint = 0.616018502797380*1.2;
%end
e_mol_gases = sum(x(1:7)) – nGin_setpoint;
F_G_R = Kc_Total_gases*(e_mol_gases+tauI_Total_gases*e_integral);
%Right-hand side evaluation of the dynamic model (DAE set)
S1 = Q1*C1-F_G_R*x(1)/sum(x(1:7))-VR*kL*(x(1)/VG-K(1)*x(8)/VL); % gas phase ethylene (mol/s)
S2 = Q2-F_G_R*x(2)/sum(x(1:7))-VR*kL*(x(2)/VG-K(2)*x(9)/VL); % gas phase butene (mol/s);
S3 = Q3-F_G_R*x(3)/sum(x(1:7))-VR*kL*(x(3)/VG-K(3)*x(10)/VL); % gas phase hexene (mol/s);
S4= Q4-F_G_R*x(4)/sum(x(1:7))-VR*kH*(x(4)/VG-K(4)*x(11)/VL); % gas phase octene (mol/s);
S5= Q5-F_G_R*x(5)/sum(x(1:7))-VR*kH*(x(5)/VG-K(5)*x(12)/VL); % gas phase decene (mol/s);
S6= Q6-F_G_R*x(6)/sum(x(1:7))-VR*kH*(x(6)/VG-K(6)*x(13)/VL); % gas phase dodecene (mol/s);
S7= Q7-F_G_R*x(7)/sum(x(1:7))-VR*kH*(x(7)/VG-K(7)*x(14)/VL); % gas phase undecane (mol/s) ;
S8= 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);
S9= 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);
S10= 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);
S11= 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);
S12= 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);
S13= 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);
S14= VR*kH*(x(7)/VG-K(7)*x(14)/VL);
S15= sum(x(1:7))-(nGin_setpoint); %Error
S = ([S1; S2; S3; S4; S5; S6; S7; S8; S9; S10; S11; S12; S13; S14; S15]);
end
Also the results are changing dramatically when I change initial values for kL and kH. I know it’s normal since fmincon() doesn’t compute global maximum/minimum but it’s just so weird to get different resluts whenever I change kL and kH values! optimization, ode MATLAB Answers — New Questions