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)
Afficher commentaires plus anciens
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));
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
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
And better use lb instead of A and b to define the condition x>=0.
Réponses (1)
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)
fun(x)
Aeq*x-beq
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
0 commentaires
Voir également
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!