*****************************************************************
*** Chapter 14
*** Section: Locating All Homogeneous Azeotropes
*** Test Problem 1
***
*** Wilson equation
*** Acetone -- Methyl Acetate -- Methanol
*****************************************************************

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

  
SCALARS 

n         number of components                             /3/

********************************************************************
* 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,
sel(i),
zlbd(i), zubd(i), templow, temphigh;

*
* Select which components participate
*
sel('1') =1;
sel('2') =1;
sel('3') =1;
* sel('4') corresponds to temperature and should always be 0
sel('4') =0;
*
*  P Pressure 
*
P = 1.0*760;
t0 = 51+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') = 7.06524*LOG(10); bb('2') = 1157.630*LOG(10); cc('2') = 219.726;
aa('3') = 8.08097*LOG(10); bb('3') = 1582.271*LOG(10); cc('3') = 239.726;
*
* Wilson binary interaction parameters 
*
a0(i,i) = 1;
a0('1','2') = -118.6763; a0('2','1') = 239.6487;
a0('1','3') = -205.3178; a0('3','1') = 658.0138;
a0('2','3') = -31.1923;  a0('3','2') = 813.1843;   

* Liquid densities
*
v('1') = 74.05;
v('2') = 79.84;
v('3') = 39.41;
*
* calculated binary interaction parameters
*
ll(i,i) = 1.0;
ll('1','2') = 0.480; ll('1','3') = 0.768;
ll('2','1') = 1.550; ll('2','3') = 0.544;
ll('3','1') = 0.566; ll('3','2') = 0.650;

templow = 0;
temphigh = 100;
*
* set variable bounds

LOOP(i$sel(i),           
zlbd(i) = 0.000;
zubd(i) = 1.000;  
    );

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.272; z.lo('1') = 0.000001; z.up('1') = 1.0;
z.l('2') = 0.465; 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') = 54.254; z.lo('4') = 20.0; z.up('4') = 80.0;

SOLVE wilson USING nlp MINIMIZING ub;