Using fmincon but converging to an infeasible point because the size of the current step is less than the value of the step size tolerance but constraints are not satisfied to

5 vues (au cours des 30 derniers jours)
clc
clear
T = 1000; %K
R = 8.314; %kJ/kmol.K
P = 300; %kpa in reformer
Po = 100; %kpa standard conditions
mair = 9.5891;
mfuel = 3;
species = {'C12H26' 'O2' 'CO' 'H2' 'CH4' 'CO2' 'C' 'H2O'};
G0 = 1e3*[-913.4878 -220.0367 -211.1858 -145.4028 -209.7328 -55.1894 -2.8287 -494720];%kJ/kmol
fun = @(x)G0*x + R*T*(x(1)*log((P/Po)*x(1)/sum(x))+ x(2)*log((P/Po)*x(2)/sum(x))+ ...
x(3)*log((P/Po)*x(3)/sum(x))+ x(4)*log((P/Po)*x(4)/sum(x))+ x(5)*log((P/Po)*x(5)/sum(x))+ ...
x(6)*log((P/Po)*x(6)/sum(x)) + x(7)*log((P/Po)*x(7)/sum(x)) + x(8)*log((P/Po)*x(8)/sum(x)));
G_CPOX = g_H2*(10) + g_CO*(9) + g_H2O + g_CH4 + g_CO2 + g_C - (g_Dod + g_O2*(6));
Unrecognized function or variable 'g_H2'.
A=zeros(8,8);
b=zeros(8,1);
for i=1:1:8
A(i,i)=-1;
end
%{C12H26 O2 - CO H2 CH4 CO2 C H2O} material balance
Aeq = [12 0 1 0 1 1 1 0 % C balance
0 2 1 0 0 2 0 1 % O balance
26 0 0 2 4 0 0 2]; % H balance
n_O2 = (mair*0.21)/32; %kmol/s
n_N2 = (mair*0.79)/28; %kmol/s
n_fuel = mfuel/170.33; %kmol/s
% molar feed of C12H26 and O2 in kmol/s
beq = [12*n_fuel % mol C atoms fed
2*n_O2 % mol O atoms fed
26*n_fuel]; % mol H atoms fed
O_C_ratio = beq(2,1)/beq(1,1);
%Solving the Gibbs minimization problem
% choose x0(1) , x0(2), x0(3), x0(4) close to the inlet conditions and solve
% for x0(4) to x0(8) using the solution set and use these as initial cndtn
% initial guesses of inlet flow rates kmol/s
x0 = [0.0100
0.0550
0.0001
0.0002
0.0493
0.0078
0.0341
0.0001];
options = optimoptions( 'fmincon' , 'Display' , 'iter' , 'Algorithm' , 'interior-point' );
x = fmincon(fun,x0,A,b,Aeq,beq, [], [], [], options);
fprintf('%10.8f O/C ratio \n ', O_C_ratio)
for i=1:numel(x)
fprintf('%5d%10s%15.5f kmol/s \n',i, species{i}, x(i))
end
K = exp(-G_CPOX/R/T); % equilibrium constant at T
fprintf('%10.8f Equilibrium constant \n ', K)
% molar fractions and partial pressures
yj = x/sum(x);
Pj = yj*P;
for i=1:numel(x)
fprintf('%5d%10s%15.4f molar fraction\n',i, species{i}, yj(i))
end
  2 commentaires
Torsten
Torsten le 13 Août 2022
Modifié(e) : Torsten le 13 Août 2022
You didn't supply all constants to test the code (see above).
Your x0 does not satisfy Aeq*x0 = beq. So choose a different x0 in order to start with a feasible solution:
mair = 9.5891;
mfuel = 3;
%{C12H26 O2 - CO H2 CH4 CO2 C H2O} material balance
Aeq = [12 0 1 0 1 1 1 0 % C balance
0 2 1 0 0 2 0 1 % O balance
26 0 0 2 4 0 0 2]; % H balance
n_O2 = (mair*0.21)/32; %kmol/s
n_N2 = (mair*0.79)/28; %kmol/s
n_fuel = mfuel/170.33; %kmol/s
% molar feed of C12H26 and O2 in kmol/s
beq = [12*n_fuel % mol C atoms fed
2*n_O2 % mol O atoms fed
26*n_fuel]; % mol H atoms fed
%Solving the Gibbs minimization problem
% choose x0(1) , x0(2), x0(3), x0(4) close to the inlet conditions and solve
% for x0(4) to x0(8) using the solution set and use these as initial cndtn
% initial guesses of inlet flow rates kmol/s
x0 = [0.0100
0.0550
0.0001
0.0002
0.0493
0.0078
0.0341
0.0001];
Aeq*x0-beq
ans = 3×1
1.0e-03 * -0.0544 -0.0569 -0.1346
And better use lb instead of A and b to define the condition x>=0.

Connectez-vous pour commenter.

Réponses (1)

Torsten
Torsten le 13 Août 2022
Modifié(e) : Torsten le 13 Août 2022
Works for me.
T = 1000; %K
R = 8.314; %kJ/kmol.K
P = 300; %kpa in reformer
Po = 100; %kpa standard conditions
mair = 9.5891;
mfuel = 3;
species = {'C12H26' 'O2' 'CO' 'H2' 'CH4' 'CO2' 'C' 'H2O'};
G0 = 1e3*[-913.4878 -220.0367 -211.1858 -145.4028 -209.7328 -55.1894 -2.8287 -494720];%kJ/kmol
fun = @(x)G0*x + R*T*(x(1)*log((P/Po)*x(1)/sum(x))+ x(2)*log((P/Po)*x(2)/sum(x))+ ...
x(3)*log((P/Po)*x(3)/sum(x))+ x(4)*log((P/Po)*x(4)/sum(x))+ x(5)*log((P/Po)*x(5)/sum(x))+ ...
x(6)*log((P/Po)*x(6)/sum(x)) + x(7)*log((P/Po)*x(7)/sum(x)) + x(8)*log((P/Po)*x(8)/sum(x)));
%G_CPOX = g_H2*(10) + g_CO*(9) + g_H2O + g_CH4 + g_CO2 + g_C - (g_Dod + g_O2*(6));
%A=zeros(8,8);
%b=zeros(8,1);
%for i=1:1:8
% A(i,i)=-1;
%end
lb = zeros(8,1);
ub = Inf*ones(8,1);
%{C12H26 O2 - CO H2 CH4 CO2 C H2O} material balance
Aeq = [12 0 1 0 1 1 1 0 % C balance
0 2 1 0 0 2 0 1 % O balance
26 0 0 2 4 0 0 2]; % H balance
n_O2 = (mair*0.21)/32; %kmol/s
n_N2 = (mair*0.79)/28; %kmol/s
n_fuel = mfuel/170.33; %kmol/s
% molar feed of C12H26 and O2 in kmol/s
beq = [12*n_fuel % mol C atoms fed
2*n_O2 % mol O atoms fed
26*n_fuel]; % mol H atoms fed
O_C_ratio = beq(2,1)/beq(1,1);
%Solving the Gibbs minimization problem
% choose x0(1) , x0(2), x0(3), x0(4) close to the inlet conditions and solve
% for x0(4) to x0(8) using the solution set and use these as initial cndtn
% initial guesses of inlet flow rates kmol/s
x0 = [0.0100
0.0550
0.0001
0.0002
0.0493
0.0078
0.0341
0.0001];
options = optimoptions( 'fmincon' , 'Display' , 'iter' , 'Algorithm' , 'interior-point' );
%x = fmincon(fun,x0,A,b,Aeq,beq, [], [], [], options);
format long g
x = fmincon(fun,x0,[],[],Aeq,beq,lb,ub,[], options)
First-order Norm of Iter F-count f(x) Feasibility optimality step 0 9 -8.202643e+04 1.346e-04 4.945e+08 1 18 -2.220561e+05 1.110e-16 4.945e+08 4.306e-03 2 27 -5.416010e+06 2.220e-16 4.945e+08 6.430e-02 3 36 -2.326957e+07 5.551e-17 3.454e+08 6.760e-02 4 45 -3.531130e+07 1.832e-15 2.436e+08 4.312e-02 5 54 -3.684624e+07 1.610e-15 2.430e+08 7.415e-03 6 63 -6.112847e+07 1.110e-16 2.311e+05 8.892e-02 7 72 -6.226494e+07 1.110e-16 9.094e+04 4.055e-03 8 81 -6.227101e+07 3.331e-16 9.023e+04 9.867e-04 9 90 -6.227126e+07 5.551e-17 6.695e+04 2.603e-03 10 99 -6.227211e+07 5.551e-17 5.161e+04 1.138e-02 11 108 -6.227494e+07 5.551e-17 3.284e+04 4.366e-02 12 117 -6.227754e+07 5.551e-17 4.864e+04 4.919e-02 13 126 -6.227935e+07 5.551e-17 4.300e+04 3.269e-02 14 135 -6.227935e+07 5.551e-17 4.335e+04 7.888e-05 15 144 -6.227936e+07 5.551e-17 4.669e+04 1.622e-04 16 156 -6.227936e+07 5.551e-17 4.747e+04 1.383e-05 17 165 -6.227936e+07 0.000e+00 1.535e+04 2.557e-06 18 174 -6.227936e+07 0.000e+00 1.879e+04 5.572e-06 19 183 -6.227936e+07 2.776e-17 2.920e+03 3.291e-06 20 192 -6.227936e+07 5.551e-17 3.755e+01 2.807e-07 Local minimum found that satisfies the constraints. Optimization completed because the objective function is non-decreasing in feasible directions, to within the value of the optimality tolerance, and constraints are satisfied to within the value of the constraint tolerance.
x = 8×1
2.30395256114726e-12 2.02268717197313e-15 4.04678605770702e-15 0.103101694470416 4.33338630718072e-06 2.02235042697261e-15 0.211350096222627 0.125856937499988
fun(x)
ans =
-62279357.3685501
Aeq*x-beq
ans = 3×1
1.0e+00 * 0 0 5.55111512312578e-17
fprintf('%10.8f O/C ratio \n ', O_C_ratio)
0.59547812 O/C ratio
for i=1:numel(x)
fprintf('%5d%10s%15.5f kmol/s \n',i, species{i}, x(i))
end
1 C12H26 0.00000 kmol/s 2 O2 0.00000 kmol/s 3 CO 0.00000 kmol/s 4 H2 0.10310 kmol/s 5 CH4 0.00000 kmol/s 6 CO2 0.00000 kmol/s 7 C 0.21135 kmol/s 8 H2O 0.12586 kmol/s
%K = exp(-G_CPOX/R/T); % equilibrium constant at T
%fprintf('%10.8f Equilibrium constant \n ', K)
% molar fractions and partial pressures
yj = x/sum(x);
Pj = yj*P;
for i=1:numel(x)
fprintf('%5d%10s%15.4f molar fraction\n',i, species{i}, yj(i))
end
1 C12H26 0.0000 molar fraction 2 O2 0.0000 molar fraction 3 CO 0.0000 molar fraction 4 H2 0.2342 molar fraction 5 CH4 0.0000 molar fraction 6 CO2 0.0000 molar fraction 7 C 0.4800 molar fraction 8 H2O 0.2858 molar fraction

Catégories

En savoir plus sur Develop Apps Using App Designer dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by