*****************************************************************
*** Chapter 14
*** Section: Locating All Homogeneous Azeotropes
*** Test Problem 7
***
*** UNIQUAC equation
*** Ethanol -- Methylcyclopentane -- Benzene -- Hexane
*****************************************************************


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

PARAMETERS
dummy, xc(i), 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);
*
*
sel('1') = 1;
sel('2') = 1;
sel('3') = 1;
sel('4') = 1;
* sel('5') corresponds to temperature and should always be 0
sel('5') = 0;
*
*  P Pressure 
*
P = 760;
*
* Antoine equation parameters
*
aa('1') = 18.9119; bb('1') = 3803.98; cc('1') = -41.68+273.15;
aa('2') = 6.8664*LOG(10); bb('2') = 1188.05*LOG(10); cc('2') = 226.276;
aa('3') = 15.9008; bb('3') = 2788.51; cc('3') = -52.36+273.15;
aa('4') = 6.91058*LOG(10); bb('4') = 1189.64*LOG(10); cc('4') = 226.28;

*
* Uniquac parameters  
*
zz = 10;                                                           
*
t0 = 60+273.15;
*
q('1') = 1.97;
q('2') = 3.01;
q('3') = 2.40;
q('4') = 3.86;
*
rr('1') = 2.11;
rr('2') = 3.97;
rr('3') = 3.19;
rr('4') = 4.50;
*
qp('1') = 0.92;
qp('2') = 3.01;
qp('3') = 2.40;
qp('4') = 3.86;
*
                      a0('1','2')=-196.545; a0('1','3')=-128.88; a0('1','4')=-142.923;
a0('2','1')=1231.703;                       a0('2','3')=160.346; a0('2','4')=-221.756;
a0('3','1')=997.41;   a0('3','2')=-65.2945;                      a0('3','4')=-26.1184;
a0('4','1')=1058.010; a0('4','2')=272.912;  a0('4','3')=125.311;
* 
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.322; z.lo('1') = 0.000001; z.up('1') = 1.0;
z.l('2') = 0.322; z.lo('2') = 0.000001; z.up('2') = 1.0;
z.l('3') = 0.222; z.lo('3') = 0.000001; z.up('3') = 1.0;
z.l('4') = 0.133; z.lo('4') = 0.000001; z.up('4') = 1.0;
z.l('5') = 63.558; z.lo('5') = 40.0; z.up('5') = 90.0;

SOLVE uniquac USING nlp MINIMIZING ub;