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

3 views (last 30 days)
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 Comments

Sign in to comment.

Answers (1)

Torsten
Torsten on 13 Aug 2022
Edited: Torsten on 13 Aug 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

Community Treasure Hunt

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

Start Hunting!