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

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

[Torsten]

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.

[Takudzwa Mashamba]

Thank you for this. Let me do so.Appreciateit!

Answer 1

[Torsten]

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