*****************************************************************
*** Chapter 14
*** Section: Locating All Homogeneous Azeotropes
*** Test Problem 6
***
*** UNIQUAC equation
*** Ethanol -- Benzene -- Water
*****************************************************************


SET
     i             variables                                 /1*4/ 
     alias(i,j)
     alias(i,l);
  
SCALARS 
n         number of components                             /3/

PARAMETERS
dummy, sel(i),
aa(i), bb(i), cc(i), P,
zz, q(i), rr(i), ll(i), qp(i), tt(i,j), a0(i,j), t0,
zi(i), templow, temphigh,
zlbd(i), zubd(i), zc(i),
sortr(i), sortq(i), sortl(i), sortqp(i), sorttqp(i,j);
*
*
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 = 760;
*
* Antoine equation parameters
*
aa('1') = 18.9119; bb('1') = 3803.98; cc('1') = -41.68+273.15;
aa('2') = 15.9008; bb('2') = 2788.51; cc('2') = -52.36+273.15;
aa('3') = 18.3036; bb('3') = 3816.44; cc('3') = -46.13+273.15; 

*
* Uniquac parameters  
*
zz = 10;                                                           
*
t0 = 60+273.15;
*
q('1') = 1.97;
q('2') = 2.40;
q('3') = 1.40;
*
rr('1') = 2.11;
rr('2') = 3.19;
rr('3') = 0.92;
*
qp('1') = 0.92;
qp('2') = 2.40;
qp('3') = 1.00;
*
a0(i,i) = 0;
                       a0('1','2') = -128.88;   a0('1','3') = -64.56;
a0('2','1') = 997.41;                           a0('2','3') = 2057.42;
a0('3','1') = 380.68;  a0('3','2') = 115.13;                        
* 
LOOP(i$(ord(i) le n),
ll(i) = zz/2*(rr(i)-q(i))-(rr(i)-1);
tt(i,i) = 1.0;
    );
LOOP((i,j)$((ord(i) le n) and (ord(j) le n) and (ord(i) ne ord(j))),
tt(i,j) = EXP(-a0(i,j)/t0);
    );
templow =  0;
temphigh = 100;
*
* set variable bounds

zi(i) = 1;
LOOP(i$(ord(i) eq n),
zi(i) = 0;
    );

LOOP(i$(ord(i) le n),
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)..

(zz/2*q(i)-1)*LOG(SUM(j$sel(j), rr(j)*z(j)))
-zz/2*q(i)*LOG(SUM(j$sel(j), q(j)*z(j)))
-rr(i)*SUM(j$sel(j), ll(j)*z(j))
      /SUM(j$sel(j), rr(j)*z(j))
-qp(i)*LOG(SUM(j$sel(j), tt(j,i)*qp(j)*z(j)))
+qp(i)*LOG(SUM(j$sel(j), qp(j)*z(j)))
-qp(i)*SUM(j$sel(j), tt(i,j)*qp(j)*z(j)/
 SUM(l$sel(l), tt(l,j)*qp(l)*z(l)))
-bb(i)/(SUM(j$(ord(j) eq n+1), z(j)) + cc(i)) +
(aa(i) - LOG(P)+LOG(rr(i))+zz/2*q(i)*LOG(q(i)/rr(i))+ll(i)+qp(i))
- a =l= 0;

eq2(i)$sel(i)..

-(zz/2*q(i)-1)*LOG(SUM(j$sel(j), rr(j)*z(j)))
+zz/2*q(i)*LOG(SUM(j$sel(j), q(j)*z(j)))
+rr(i)*SUM(j$sel(j), ll(j)*z(j))
      /SUM(j$sel(j), rr(j)*z(j))
+qp(i)*LOG(SUM(j$sel(j), tt(j,i)*qp(j)*z(j)))
-qp(i)*LOG(SUM(j$sel(j), qp(j)*z(j)))
+qp(i)*SUM(j$sel(j), tt(i,j)*qp(j)*z(j)/
 SUM(l$sel(l), tt(l,j)*qp(l)*z(l)))
+bb(i)/(SUM(j$(ord(j) eq n+1), z(j)) + cc(i))
- (aa(i) - LOG(P)+LOG(rr(i))+zz/2*q(i)*LOG(q(i)/rr(i))+ll(i)+qp(i))
- a =l= 0;

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

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

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

a.lo = 0.0;
z.l('1') = 0.013; z.lo('1') = 0.000001; z.up('1') = 1.0;
z.l('2') = 0.604; z.lo('2') = 0.000001; z.up('2') = 1.0;
z.l('3') = 0.383; z.lo('3') = 0.000001; z.up('3') = 1.0;
z.l('4') = 61.583; z.lo('4') = 40.0; z.up('4') = 90.0;

SOLVE uniquac USING nlp MINIMIZING ub;