*****************************************************************
*** Chapter 14
*** Section: Locating All Homogeneous Azeotropes
*** Test Problem 4
***
*** NRTL equation
*** Ethanol -- Methyl Ethyl Ketone -- Water
*****************************************************************


SET
     i             variables                                 /1*4/ 
     alias(i,j)
     alias(i,l);

SCALARS 

n         number of components                             /3/


PARAMETERS
aa(i), bb(i), cc(i), ll(i,j), v(i), P, t0, Rg,
zi(i), templow, temphigh,
dum,
sel(i),
alpha(i,j), a0(i,j), tt(i,j), g(i,j),
zlbd(i), zubd(i), zc(i), sorty(i,j), sortw(i,j);

*
* Select which components participate
*
sel('1') =1;
sel('2') =1;
sel('3') =1;
*
*  P Pressure 
*
P = 1*760;
t0 = 73.39+273.15;
Rg = 1.98721;
*
* Antoine equation parameters
*
aa('1') = 8.11220*LOG(10); bb('1') = 1592.864*LOG(10); cc('1') = 226.184;
aa('2') = 7.06356*LOG(10); bb('2') = 1261.340*LOG(10); cc('2') = 221.969;
aa('3') = 8.07131*LOG(10); bb('3') = 1730.630*LOG(10); cc('3') = 233.426;

*
* NRTL parameters 
*
a0('1','2') =  64.4957;  a0('2','1') = 463.1931;  
a0('1','3') =  -175.0164;  a0('3','1') = 1440.3479; 
a0('2','3') = 674.4614; a0('3','2') = 1809.8868; 

alpha('1','2') = 0.3010;                                
alpha('1','3') = 0.2959;                                
alpha('2','3') = 0.3536;                             

tt(i,i) = 0;
LOOP((i,j)$((ord(i) le n) and (ord(j) le n) and (ord(i) ne ord(j))),
tt(i,j) = a0(i,j)/(Rg*t0);
    );
*
*
g(i,i) = 1;
LOOP((i,j)$((ord(i) le n) and (ord(j) le n) and (ord(i) ne ord(j))),
g(i,j) = EXP(-(alpha(i,j)+alpha(j,i))*tt(i,j));
    );
templow =  0;
temphigh = 100;

* set variable bounds


LOOP(i$sel(i),           
zlbd(i) = 0;
zubd(i) = 1;  
    );

LOOP(i$(ord(i) gt n),
zlbd(i) = templow;  
zubd(i) = temphigh
    );


**************************************************************
* z(i) = mole fractions and temperature
* ub = objective value
* a = slack variable

VARIABLES z(i), ub, a;
POSITIVE VARIABLE  a;

**************************************************************
* obj = objective (minimize maximum equality constraint violation
* eq1(i) = equilibrium constraint 1
* eq2(i) = equilibrium constraint 2
* sumi = mole fractions sum to 1

EQUATIONS
         obj
	 eq1(i)   Equilibrium  
	 eq2(i)   Equilibrium  
	 sumi;


obj..    ub  =e= a;

eq1(i)$sel(i)..

 SUM(j$sel(j), tt(j,i)*g(j,i)*z(j))/SUM(j$sel(j), g(j,i)*z(j))
+SUM(j$sel(j), tt(i,j)*g(i,j)*z(j)/SUM(l$sel(l), g(l,j)*z(l)))
-SUM(j$sel(j), g(i,j)*z(j)*SUM(l$sel(l), tt(l,j)*g(l,j)*z(l))
/SQR(SUM(l$sel(l), g(l,j)*z(l))))
-bb(i)/(SUM(j$(ord(j) gt n), z(j)) + cc(i))
+(aa(i) - LOG(P))
- a =l= 0;

eq2(i)$sel(i)..          

-SUM(j$sel(j), tt(j,i)*g(j,i)*z(j))/SUM(j$sel(j),g(j,i)*z(j))
-SUM(j$sel(j), tt(i,j)*g(i,j)*z(j)/SUM(l$sel(l), g(l,j)*z(l)))
+SUM(j$sel(j), g(i,j)*z(j)*SUM(l$sel(l), tt(l,j)*g(l,j)*z(l))
/SQR(SUM(l$sel(l), g(l,j)*z(l))))
+bb(i)/(SUM(j$(ord(j) gt n), z(j)) + cc(i))
-(aa(i) - LOG(P))
- a =l= 0;

sumi.. SUM(i$sel(i), z(i)) =e= 1;

*-----------------------------------------------------------------------

MODEL nrtl / obj, eq1, eq2, sumi /;

a.lo = 0.0;
z.l('1') = 0.187; z.lo('1') = 0.000001; z.up('1') = 1.0;
z.l('2') = 0.560; z.lo('2') = 0.000001; z.up('2') = 1.0;
z.l('3') = 0.253; z.lo('3') = 0.000001; z.up('3') = 1.0;
z.l('4') = 72.957; z.lo('4') = 40.0; z.up('4') = 90.0;

SOLVE nrtl USING nlp MINIMIZING ub;