*****************************************************************
*** Chapter 14
*** Section: Locating All Homogeneous Azeotropes
*** Test Problem 3
***
*** Wilson equation
*** Acetone -- Chloroform -- Methanol -- Ethanol -- Benzene
*****************************************************************


SET
     i             variables                                 /1*6/ 
     alias(i,j)
     alias(i,k);

SCALARS 

n         number of components                             /5/

********************************************************************
* P = pressure (mmHg)
* Rg = gas constant
* t0 = reference temperature
* aa(i), bb(i), cc(i) = Antoine parameters
* v(i) = liquid densities
* a0(i,j), ll(i,j) = binary interaction parameters
* zlbd(i), zubd(i), templow, temphigh = bounds on variables

PARAMETERS
aa(i), bb(i), cc(i), a0(i,j), ll(i,j), v(i), P, t0, Rg,
zlbd(i), zubd(i), templow, temphigh, sel(i);

*
* Select which components participate
*
sel('1') =1;
sel('2') =1;
sel('3') =1;
sel('4') =0;
sel('5') =1;
* sel('6') corresponds to temperature and should always be 0
sel('6') = 0;
*
*  P Pressure 
*
P = 1*760;
t0 = 65.09+273.15;
Rg = 1.98721;
*
* Antoine equation parameters
*
aa('1') = 7.11714*LOG(10); bb('1') = 1210.595*LOG(10); cc('1') = 229.664;       
aa('2') = 6.95465*LOG(10); bb('2') = 1170.966*LOG(10); cc('2') = 226.232;       
aa('3') = 8.08097*LOG(10); bb('3') = 1582.271*LOG(10); cc('3') = 239.726;
aa('4') = 8.11220*LOG(10); bb('4') = 1592.864*LOG(10); cc('4') = 226.184;       
aa('5') = 6.87987*LOG(10); bb('5') = 1196.760*LOG(10); cc('5') = 219.161;

*
* Wilson parameters 
*
a0(i,i) = 1.0;
a0('1','2') = 28.8819; a0('2','1') = -484.3856;
a0('1','3') = -161.8813; a0('3','1') = 583.1054; 
a0('1','4') = 180.4365; a0('4','1') = 252.2864; 
a0('1','5') = 543.9352; a0('5','1') = -182.5230;
a0('2','3') = -351.1964; a0('3','2') = 1760.6741;  
a0('2','4') = -315.6205; a0('4','2') = 1434.1695;  
a0('2','5') = -161.8065; a0('5','2') = 49.6010;    
a0('3','4') =  -65.7022; a0('4','3') = 143.6658;  
a0('3','5') =  1717.8985; a0('5','3') = 196.8233;   
a0('4','5') =  1399.9279; a0('5','4') = 207.3433;   

*
*
v('1') = 71.1633;
v('2') = 94.2678;
v('3') = 38.9474;
v('4') = 59.2123;
v('5') = 94.3214;
*
*
ll(i,i) = 1;
LOOP((i,j)$((ord(i) le n) and (ord(j) le n) and (ord(i) ne ord(j))),
ll(i,j) = v(j)/v(i)*EXP(-a0(i,j)/(Rg*t0));
    );

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)..

LOG(SUM(j$sel(j), z(j)*ll(i,j)))
+SUM(j$sel(j), z(j)*ll(j,i) / SUM(k$sel(k), z(k)*ll(j,k)))
+bb(i)/(SUM(j$(ord(j) eq n+1), z(j)) + cc(i))
-aa(i) - 1.0 + LOG(P)
- a =l= 0;

eq2(i)$sel(i)..
-LOG(SUM(j$sel(j), z(j)*ll(i,j)))
-SUM(j$sel(j), z(j)*ll(j,i) / SUM(k$sel(k), z(k)*ll(j,k)))
-bb(i)/(SUM(j$(ord(j) eq n+1), z(j)) + cc(i))
+aa(i) + 1.0 - LOG(P)
- a =l= 0;

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

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

MODEL wilson / all /;

a.lo = 0.0;
z.l('1') = 0.295; z.lo('1') = 0.000001; z.up('1') = 1.0;
z.l('2') = 0.148; z.lo('2') = 0.000001; z.up('2') = 1.0;
z.l('3') = 0.463; z.lo('3') = 0.000001; z.up('3') = 1.0;
z.l('4') = 0.000001; z.lo('4') = 0.000001; z.up('4') = 1.0;
z.l('5') = 0.094; z.lo('5') = 0.000001; z.up('5') = 1.0;
z.l('6') = 57.154; z.lo('6') = 20.0; z.up('6') = 80.0;

SOLVE wilson USING nlp MINIMIZING ub;