I am trying to solve non-linear probrem (NLP) using GMAS. However, I am getting output "Infeasible solution".The material balance equation in chemical engineering, which describe in this paper.
$title a membrane separation model
Sets
gas component gas name / Gas1, Gas2, Gas3, Gas4 /
m membrane module / module1, module2 /;
Parameters
q(gas) permeance Q
/ Gas1 89e-4
Gas2 71.2e-4
Gas3 4.45e-4
Gas4 1.78e-4 /
x_init_gas(gas) process feed component
/ Gas1 0.19
Gas2 0.01
Gas3 0.73
Gas4 0.07 /
membrane_area(m) area settings
/ module1 222.92
module2 164.45 /
high_pressure(m) feed side high pressure
/ module1 3.5
module2 3.5 /
low_pressure(m) feed side high pressure
/ module1 0.105
module2 0.105 / ;
Scalar
U fresh feed flow rate /90 /
Rg /8.3144 /
temperature /313 /
f_mh /200 /
f_cp /1 /
etha_cp /0.7 /
f_cc /0.27 /
f_wk /0.1 /
f_mr /90 /
t_m /0.05 /
f_mt /0.05 /
f_sg /35e-9 /
t_wk /300 /
f_hv /43e6 /
gamma_LB lower bounds for pressure ratio /0.01/
gamma_UB upper bounds for pressure ratio /20/;
Variables
B(m, gas)
E(m)
D(m)
T(m, gas)
Tdot(m, gas)
Tdotdot(m, gas)
total_annual_cost ;
T.l(m, gas) = 0.8 ;
Tdot.l(m, gas) = 0.8 ;
Tdotdot.l(m, gas) = 0.8 ;
B.l(m, gas) = 1.5 ;
E.l(m) = 1 ;
D.l(m) = 1e-4 ;
T.lo(m, gas) = 0 ;
Tdot.lo(m, gas) = 0 ;
Tdotdot.lo(m, gas) = 0 ;
D.lo(m) = min(q('Gas1'), q('Gas2'), q('Gas3'), q('Gas4')) ;
T.up(m, gas) = 1-gamma_LB ;
Tdot.up(m, gas) = 1-gamma_LB ;
Tdotdot.up(m, gas) = 1-gamma_LB ;
D.up(m) = max(q('Gas1'), q('Gas2'), q('Gas3'), q('Gas4')) ;
Positive Variable
F_x(m, gas)
F(m)
xF(m, gas)
L_x(m, gas)
L(m)
xR(m, gas)
V_x(m, gas)
V(m)
yP(m, gas)
Lret_x(gas)
Lret
xRet(gas)
Vper_x(gas)
Vper
yPer(gas)
Mix_x(m, gas)
Mix(m)
x_Mix(m, gas);
F_x.l(m, gas) = 5 ;
F.l(m) = 20 ;
xF.l(m, gas) = 0.25 ;
L_x.l(m, gas) = 5 ;
L.l(m) = 20 ;
xR.l(m, gas) = 0.25 ;
V_x.l(m, gas) = 2.5 ;
V.l(m) = 10 ;
yP.l(m, gas) = 0.25 ;
Lret_x.l(gas) = 5 ;
Lret.l = 20 ;
xRet.l(gas) = 0.25 ;
Vper_x.l(gas) = 2.5 ;
Vper.l = 10 ;
yPer.l(gas) = 0.25 ;
Mix_x.l(m, gas) = 10 ;
Mix.l(m) = 30 ;
x_Mix.l(m, gas) = 0.25 ;
F_x.lo(m, gas) = 0 ;
F.lo(m) = 0 ;
xF.lo(m, gas) = 0 ;
L_x.lo(m, gas) = 0 ;
L.lo(m) = 0 ;
xR.lo(m, gas) = 0 ;
V_x.lo(m, gas) = 0 ;
V.lo(m) = 0 ;
yP.lo(m, gas) = 0 ;
Lret_x.lo(gas) = 0 ;
Lret.lo = 0 ;
xRet.lo(gas) = 0 ;
Vper_x.lo(gas) = 0 ;
Vper.lo = 0 ;
yPer.lo(gas) = 0 ;
Mix_x.lo(m, gas) = 0 ;
Mix.lo(m) = 0 ;
x_Mix.lo(m, gas) = 0 ;
F_x.up(m, gas) = 30 ;
F.up(m) = 30 ;
xF.up(m, gas) = 1.0 ;
L_x.up(m, gas) = 30 ;
L.up(m) = 30 ;
xR.up(m, gas) = 1.0 ;
V_x.up(m, gas) = 30 ;
V.up(m) = 30 ;
yP.up(m, gas) = 1.0 ;
Lret_x.up(gas) = 30 ;
Lret.up = 30 ;
xRet.up(gas) = 1.0 ;
Vper_x.up(gas) = 30 ;
Vper.up = 30 ;
yPer.up(gas) = 1.0 ;
Mix_x.up(m, gas) = 30 ;
Mix.up(m) = 40 ;
x_Mix.up(m, gas) = 1.0 ;
Parameter
u_x(gas) feed flow rate for component g
G(m) parameter G ;
u_x(gas) = U * x_init_gas(gas) ;
G(m) = low_pressure(m) / high_pressure(m) ;
Equations
cost define objective function
const_module_1(m, gas)
const_module_2(m)
const_module_3(m)
const_module_4(m, gas)
const_module_5(m)
const_module_6(m, gas)
const_module_7(m)
const_module_8(m, gas)
const_module_9(m, gas)
const_module_10(m, gas)
const_module_11(m, gas)
const_module_12(m, gas)
const_module_13(m, gas)
const_module_14(m)
const_module_15(m)
const_module_16(m)
const_module_17(m)
const_module_18(m)
const_module_19(m)
tight_constraint_1(m, gas)
tight_constraint_2(m, gas)
tight_constraint_3(gas)
tight_constraint_4(m, gas)
tight_constraint_5(m, gas)
tight_constraint_6(m, gas)
const_mixer1_1(gas)
const_mixer1_2
const_mixer1_3(gas)
const_mixer2_1(gas)
const_mixer2_2
const_mixer2_3(gas)
const_mixer_module_1(m, gas)
const_mixer_module_2(m)
const_mixer_module_3(m, gas)
const_mixer_retentate_1(gas)
const_mixer_retentate_2
const_mixer_retentate_3(gas)
const_mixer_permeate_1(gas)
const_mixer_permeate_2
const_mixer_permeate_3(gas)
const_component_mixer(m)
const_component_retentate
const_component_permeate
const_overall_1(gas)
const_overall_2
const_overall_3(gas);
cost .. total_annual_cost =e= (f_cc*(1+f_wk)*f_mh*sum(m, membrane_area(m))+f_cp*(Rg*temperature*sum(gas, V_x("module2", gas))*log(high_pressure("module1")/low_pressure("module2")))/etha_cp
+f_mr / t_m * (sum(m, membrane_area(m)))
+f_mt * f_mh*sum(m, membrane_area(m))+f_cp*(Rg*temperature*sum(gas, V_x("module2", gas))*log(high_pressure("module1")/low_pressure("module2")))/etha_cp
+(f_sg*t_wk)/(f_hv*etha_cp)*(Rg*temperature*sum(gas, V_x("module2", gas))*log(high_pressure("module1")/low_pressure("module2")))*24*3600
+f_sg*t_wk*sum(gas, Vper_x(gas))*yPer("Gas1")/(xRet("Gas1")+1e-3))
/(U*t_wk);
const_module_1(m, gas) .. B(m, gas) =e= 1 + q(gas) * E(m) ;
const_module_2(m) .. low_pressure(m)/high_pressure(m) =e= D(m) * E(m) ;
const_module_3(m) .. D(m) =e= sum(gas, q(gas)*T(m, gas)) ;
const_module_4(m, gas) .. xR(m, gas) =e= T(m, gas) * B(m, gas) ;
const_module_5(m) .. sum(gas, Tdotdot(m, gas)) =e= 1 - low_pressure(m)/high_pressure(m) ;
const_module_6(m, gas) .. Tdot(m, gas) =e= xF(m, gas) - low_pressure(m)/high_pressure(m) * yP(m, gas) ;
const_module_7(m) .. sum(gas, Tdot(m, gas)) =e= 1 - low_pressure(m)/high_pressure(m) ;
const_module_8(m, gas) .. 2*T(m, gas)**0.3275 =e= Tdot(m, gas)**0.3275 + Tdotdot(m, gas)**0.3275 ;
const_module_9(m, gas) .. V_x(m, gas) =e= q(gas) * T(m, gas) * high_pressure(m) * membrane_area(m) ;
const_module_10(m, gas) .. F_x(m, gas) =e= L_x(m, gas) + V_x(m, gas) ;
const_module_11(m, gas) .. F_x(m, gas) =e= xF(m, gas) * F(m) ;
const_module_12(m, gas) .. L_x(m, gas) =e= xR(m, gas) * L(m) ;
const_module_13(m, gas) .. V_x(m, gas) =e= yP(m, gas) * V(m) ;
const_module_14(m) .. F(m) =e= sum(gas, F_x(m, gas)) ;
const_module_15(m) .. L(m) =e= sum(gas, L_x(m, gas)) ;
const_module_16(m) .. V(m) =e= sum(gas, V_x(m, gas)) ;
const_module_17(m) .. sum(gas, xF(m, gas)) =e= 1 ;
const_module_18(m) .. sum(gas, xR(m, gas)) =e= 1 ;
const_module_19(m) .. sum(gas, yP(m, gas)) =e= 1 ;
tight_constraint_1(m, gas) .. xR(m, gas) =l= xF(m, gas) ;
tight_constraint_2(m, gas) .. xF(m, gas) =l= yP(m, gas) ;
tight_constraint_3(gas) .. B("module1", gas) =g= B("module2", gas) ;
tight_constraint_4(m, gas) .. xF(m, gas) - gamma_UB * yP(m, gas) =l= Tdot(m, gas) ;
tight_constraint_5(m, gas) .. Tdot(m, gas) =l= xF(m, gas) - gamma_LB * yP(m, gas) ;
tight_constraint_6(m, gas) .. T(m, gas) =l= (Tdot(m, gas) + Tdotdot(m, gas)) / 2 ;
const_mixer1_1(gas) .. Mix_x("module1", gas) =e= U*x_init_gas(gas) + V_x("module2", gas) ;
const_mixer1_2 .. Mix("module1") =e= U + V("module2") ;
const_mixer1_3(gas) .. Mix("module1")*x_Mix("module1", gas) =e= U*x_init_gas(gas) + V("module2")*yP("module2", gas) ;
const_mixer2_1(gas) .. Mix_x("module2", gas) =e= V_x("module1", gas) ;
const_mixer2_2 .. Mix("module2") =e= V("module1") ;
const_mixer2_3(gas) .. Mix("module2")*x_Mix("module2", gas) =e= V("module1")*yP("module1", gas) ;
const_mixer_module_1(m, gas).. Mix_x(m, gas) =e= F_x(m, gas) ;
const_mixer_module_2(m) .. Mix(m) =e= F(m) ;
const_mixer_module_3(m, gas) .. Mix(m)*x_Mix(m, gas) =e= F(m) * xF(m, gas) ;
const_mixer_retentate_1(gas) .. Lret_x(gas) =e= L_x("module2", gas) ;
const_mixer_retentate_2 .. Lret =e= L("module2") ;
const_mixer_retentate_3(gas) .. Lret * xRet(gas) =e= L("module2") * xR("module2", gas) ;
const_mixer_permeate_1(gas) .. Vper_x(gas) =e= V_x("module1", gas) ;
const_mixer_permeate_2 .. Vper =e= V("module1") ;
const_mixer_permeate_3(gas) .. Vper*yPer(gas) =e= V("module1") * yP("module1", gas) ;
const_component_mixer(m) .. sum(gas, x_Mix(m, gas)) =e= 1 ;
const_component_retentate .. sum(gas, xRet(gas)) =e= 1 ;
const_component_permeate .. sum(gas, yPer(gas)) =e= 1 ;
const_overall_1(gas) .. U * x_init_gas(gas) =e= Lret_x(gas) + Vper_x(gas) ;
const_overall_2 .. U =e= Lret + Vper ;
const_overall_3(gas) .. U * x_init_gas(gas) =e= Lret*xRet(gas) + Vper*yPer(gas) ;
Model transport /all/ ;
Solve transport using nlp minimizing total_annual_cost ;
Display F.m ;
Output describes that this problem is infeasible. I have no idea to solve this error. Could you please help me solve it? Thank you very much!
