- Increase MaxFunctionEvaluations and MaxIterations.
- Add FunctionTolerance and ConstraintTolerance.
- Use GlobalSearch for potentially better global optimization results.
- Adjust the optimoptions to include these new parameters.
Improving optimization results(Fmincon)
조회 수: 1 (최근 30일)
이전 댓글 표시
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!
댓글 수: 0
채택된 답변
Manikanta Aditya
2024년 5월 21일
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.
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
추가 답변 (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!