Improving optimization results(Fmincon)

1 次查看（过去 30 天）
naiva saeedia
naiva saeedia 2024-5-21
评论： Manikanta Aditya 2024-5-26
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!

请先登录，再回答此问题。

采纳的回答

Manikanta Aditya
Manikanta Aditya 2024-5-21
Hi @naiva saeedia,
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 个评论
显示 8更早的评论

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Surrogate Optimization 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by