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


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

PARAMETERS
sel(i),
dummy, xc(i),
aa(i), bb(i), cc(i), P,
zz, q(i), rr(i), ll(i), br(g), bq(g), a(g,l), v(i,g),
qv(g,i), phi(g,m), alph(g,i), beta(g,i,j), delt(i),
zi(i), templow, temphigh, t0,
zlbd(i), zubd(i);
*
* select components that participate
*
sel('1') = 1;
sel('2') = 1;
sel('3') = 0;
sel('4') = 0;
* sel('5') corresponds to temperature and should always be 0
sel('5') = 0;
*
*  P Pressure 
*
P = 760;
t0 = 60+273.15;
*
* 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;                                                           
*
bq('1') = 0.848;
bq('2') = 0.540;
bq('3') = 0.228;
bq('4') = 1.124;
bq('5') = 0.400;
*
br('1') = 0.9011;
br('2') = 0.6744;
br('3') = 0.4469;
br('4') = 1.2044;
br('5') = 0.5313;
*
v(i,g)=0;
v('1','1')=1;                             v('1','4')=1;
v('2','1')=1; v('2','2')=4; v('2','3')=1;
                                                        v('3','5')=6;
v('4','1')=2; v('4','2')=4;
*
a('1','1')=0; a('1','2')=0; a('1','3')=0; a('1','4')=931.2; a('1','5')=32.08;
a('2','1')=0; a('2','2')=0; a('2','3')=0; a('2','4')=931.2; a('2','5')=32.08;
a('3','1')=0; a('3','2')=0; a('3','3')=0; a('3','4')=931.2; a('3','5')=32.08;
a('4','1')=169.7; a('4','2')=169.7; a('4','3')=169.7; a('4','4')=0; a('4','5')=83.5;
a('5','1')=15.26; a('5','2')=15.26; a('5','3')=15.26; a('5','4')=705.9; a('5','5')=0;
*
q(i)=SUM(g, v(i,g)*bq(g));
rr(i)=SUM(g, v(i,g)*br(g));
* 
LOOP(i$(ord(i) le n),
ll(i) = zz/2*(rr(i)-q(i))-(rr(i)-1);
    );
qv(g,i)=0;
LOOP(i$(ord(i) le n),
qv(g,i)=bq(g)*v(i,g);
    );
LOOP((g,m),
phi(g,m)=EXP(-a(g,m)/t0);
     );
LOOP(g,LOOP(i$(ord(i) le n),
alph(g,i)=SUM(m, phi(m,g)*qv(m,i));
     ););
LOOP(g,LOOP(i$(ord(i) le n),LOOP(j$(ord(j) le n),
beta(g,i,j)=SUM(m, qv(m,i)*phi(m,g)*qv(g,j))
     );););
LOOP(i$(ord(i) le n),
delt(i)=aa(i)
       -LOG(P)
       +LOG(rr(i))
       +zz/2*q(i)*LOG(q(i)/rr(i))
       +ll(i)
       +SUM(g, qv(g,i)*( LOG(SUM(m, phi(m,g)*qv(m,i)))
			-LOG(q(i))
			+SUM(m, phi(g,m)*qv(m,i)
			       /SUM(l, phi(l,m)*qv(l,i))) ) );
     );
*
templow =  50;  
temphigh = 80; 
*
* 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, s;
POSITIVE VARIABLE  s;

**************************************************************
* 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= s;

eq1(i)$sel(i)..

(zz/2*q(i)-1)*LOG(SUM(j$sel(j), rr(j)*z(j)))
-(zz/2-1)*q(i)*LOG(SUM(j$sel(j), q(j)*z(j)))
-rr(i)*SUM(j$sel(j), ll(j)*z(j))
      /SUM(k$sel(k), rr(k)*z(k))
-SUM(l, qv(l,i)*LOG( SUM(j$sel(j), alph(l,j)*z(j)) ) )
-SUM(l, SUM(j$sel(j), beta(l,i,j)*z(j))
       /SUM(k$sel(k), alph(l,k)*z(k)) )
- bb(i)/(SUM(j$(ord(j) eq n+1), z(j)) + cc(i))
+ delt(i)
- s =l= 0;

eq2(i)$sel(i)..

-(zz/2*q(i)-1)*LOG(SUM(j$sel(j), rr(j)*z(j)))
+(zz/2-1)*q(i)*LOG(SUM(j$sel(j), q(j)*z(j)))
+rr(i)*SUM(j$sel(j), ll(j)*z(j))
      /SUM(k$sel(k), rr(k)*z(k))
+SUM(l, qv(l,i)*LOG( SUM(j$sel(j), alph(l,j)*z(j)) ) )
+SUM(l, SUM(j$sel(j), beta(l,i,j)*z(j))
       /SUM(k$sel(k), alph(l,k)*z(k)) )
+ bb(i)/(SUM(j$(ord(j) eq n+1), z(j)) + cc(i))
- delt(i)
- s =l= 0;


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

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

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

s.lo = 0.0;
z.l('1') = 0.371; z.lo('1') = 0.000001; z.up('1') = 1.0;
z.l('2') = 0.629; z.lo('2') = 0.000001; z.up('2') = 1.0;
z.l('3') = 0.000001; 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') = 60.632; z.lo('5') = 40.0; z.up('5') = 90.0;

SOLVE unifac USING nlp MINIMIZING ub;