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!

 Risposta accettata

After going through the codes you provided, I can suggest you to try some workarounds and see if the code gets optimized as you expect.
  1. Increase MaxFunctionEvaluations and MaxIterations.
  2. Add FunctionTolerance and ConstraintTolerance.
  3. Use GlobalSearch for potentially better global optimization results.
  4. Adjust the optimoptions to include these new parameters.
clear variables
clc
Objective = @MassTransferErrors_Closed_loop;
A = [];
b = [];
Aeq = [];
beq = [];
lb = [0 ; 0];
ub = [1000;1000];
kL = 500;
kH = 500;
p0 = [kL, kH];
% Optimoptions with improved parameters
options = optimoptions('fmincon', ...
'Display', 'iter', ...
'Algorithm', 'sqp-legacy', ...
'StepTolerance', 1e-11, ...
'MaxFunctionEvaluations', 5e4, ... % Increased function evaluations
'FunctionTolerance', 1e-12, ...
'ConstraintTolerance', 1e-6, ...
'MaxIterations', 1e4, ...
'ScaleProblem', true);
% Using GlobalSearch for potentially better global optimization results
gs = GlobalSearch;
problem = createOptimProblem('fmincon', 'x0', p0, 'objective', Objective, ...
'lb', lb, 'ub', ub, 'options', options);
k = run(gs, problem);
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 ethylene 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];
[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 = sum(y(:, 1:7), 2);
molGend = y(end, 1:7);
molLend = y(end, 8:14);
massGend = molGend' .* moleWt;
massLend = molLend' .* moleWt;
TotalProduct = massGend + massLend;
Experiment_i = cell2mat(Experiments(i));
MTE_i = sum(((TotalProduct(2:6) - Experiment_i') ./ Experiment_i').^2);
MTE_j(i) = MTE_i;
end
MTE = sum(MTE_j);
end
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
Q1 = Q0 * 1e-6 / 60; % Q_G ethylene inflow (m3/s)
P1 = 36e5; % ethylene inflow pressure [Pa]
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)]
%% Setpoint
nGin_setpoint = 0.363839693638512;
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 = - F_G_R * x(2) / sum(x(1:7)) - VR * kH * (x(2) / VG - K(2) * x(9) / VL) + wc * k(1) * x(8); % gas phase butene (mol/s)
S3 = - F_G_R * x(3) / sum(x(1:7)) - VR * kH * (x(3) / VG - K(3) * x(10) / VL) + wc * k(2) * x(9); % gas phase hexene (mol/s)
S4 = - F_G_R * x(4) / sum(x(1:7)) - VR * kH * (x(4) / VG - K(4) * x(11) / VL) + wc * k(3) * x(10); % gas phase octene (mol/s)
S5 = - F_G_R * x(5) / sum(x(1:7)) - VR * kH * (x(5) / VG - K(5) * x(12) / VL) + wc * k(4) * x(11); % gas phase decene (mol/s)
S6 = - F_G_R * x(6) / sum(x(1:7)) - VR * kH * (x(6) / VG - K(6) * x(13) / VL) + wc * k(5) * x(12); % gas phase dodecene (mol/s)
S7 = - F_G_R * x(7) / sum(x(1:7)) - VR * kH * (x(7) / VG - K(7) * x(14) / VL) + wc * k(6) * x(13); % gas phase undecene (mol/s)
S8 = VR * kL * (x(1) / VG - K(1) * x(8) / VL) - wc * k(1) * x(8); % liquid phase ethylene (mol/s)
S9 = VR * kH * (x(2) / VG - K(2) * x(9) / VL) - wc * k(2) * x(9); % liquid phase butene (mol/s)
S10 = VR * kH * (x(3) / VG - K(3) * x(10) / VL) - wc * k(3) * x(10); % liquid phase hexene (mol/s)
S11 = VR * kH * (x(4) / VG - K(4) * x(11) / VL) - wc * k(4) * x(11); % liquid phase octene (mol/s)
S12 = VR * kH * (x(5) / VG - K(5) * x(12) / VL) - wc * k(5) * x(12); % liquid phase decene (mol/s)
S13 = VR * kH * (x(6) / VG - K(6) * x(13) / VL) - wc * k(6) * x(13); % liquid phase dodecene (mol/s)
S14 = VR * kH * (x(7) / VG - K(7) * x(14) / VL) - wc * k(7) * x(14); % liquid phase undecene (mol/s)
S15 = e_mol_gases; % integral of error (mol/s)
S = [S1; S2; S3; S4; S5; S6; S7; S8; S9; S10; S11; S12; S13; S14; S15]; % the dynamic set
end
I hope this will help you.

10 Commenti

naiva saeedia
naiva saeedia il 21 Mag 2024
Modificato: naiva saeedia il 21 Mag 2024
Hi. Thankyou so much for your answer. I have run your code but I kept get this message:
Until it says:
So it gave me kL=997.1827 and kH=534.6976..I tested these values in my code the results are better than my first attempt but they also gave me big error. What's the problem in first pic?? It seems like my initial guess stopped the optimization..Am I correct??
Manikanta Aditya
Manikanta Aditya il 21 Mag 2024
Modificato: Manikanta Aditya il 21 Mag 2024
Good to know that the results are better than your first attempt.
Yes, you are correct. The message indicates that the optimization process stopped because the initial point is already a local minimum that satisfies the constraints. This means the optimizer found that the objective function does not decrease in feasible directions from the initial guess, so it concluded that it had already found a solution.
Hope this helps.
Thanks for your answer but I get this message for every initial value that I gave to kL and kH. Is this even normal??
It's unusual for an optimization problem to consistently identify the initial point as a local minimum for a variety of initial guesses, especially in a complex system like the one you're working with. This suggests there might be an issue with the objective function or the setup of the optimization problem.
Check the Objective Function:
  1. Ensure the objective function MassTransferErrors_Closed_loop is correctly implemented and sensitive to changes in kL and kH.
  2. Visualize the objective function with respect to kL and kH to understand its behavior. You can plot the function for a range of kL and kH values.
Verify that the model parameters (kL, kH) significantly impact the outputs of the simulation. If changes in these parameters do not lead to meaningful differences in the objective function, the optimizer will struggle to find a minimum.
Thanks a lot for your help. I guess my problem is that kL and kH doesn't change the system behaviour dramatically as you said. I'm going to check my model and do some adjustment to solve this problem.
No mention!
Hi again. I have run my function 1000 times for random kL and kH between 0 to 1000. Here is the plot of the objective function in each run:
As you can see I have a lot of oscillations but majority of them are not big enough. What do you think about this plot? Do you think I should change my model completely?
Most of these oscillations appear to be relatively small in magnitude. This could suggest that while there is variability in the objective function outcomes, there may be a level of consistency or stability in the results as well.
I will recommend you to try out few things before changing the model completely:
  1. Parameter Tuning: The oscillations could be due to the range of kL and kH values you’re using. Try to experiment with different ranges or distributions of these parameters.
  2. Convergence Criteria: Check your model’s convergence criteria. If the model stops too early, it might not reach the optimal solution, leading to oscillations.
  3. Model Complexity: If your model is too complex, it might overfit to the training data and perform poorly on unseen data. Try simplifying your model or using regularization techniques to prevent overfitting.
Thanks! Sorry for my late responce.
Your welcome!

Accedi per commentare.

Più risposte (0)

Community Treasure Hunt

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

Start Hunting!

Translated by