*****************************************************************
*** Chapter 6
*** Test Problem 3
***
*** Toluene - Water - Aniline -- Gibbs energy minimization(NRTL)
*****************************************************************

SETS
    i     components                /1*3/
    k     phases                    /1*2/
    alias(i,j);

*****************************************************************
* P = pressure (atm)
* T = temperature (K)
* liqphase(k) = is phase k a liquid phase?
* ntot(i) = number of moles of component i in feed
* alpha(i,j) = binary interaction parameter (symmetric)
* tau(i,j) = binary interaction parameter (nonsymmetric)
* g(i,j) = binary interaction parameter (calculated)

PARAMETERS
P, T, liqphase(k), ntot(i), alpha(i,j), tau(i,j), g(i,j);

P = 1.0;
T = 298;

liqphase('1') = 1; liqphase('2') = 1;

ntot('1') = 0.2995;
ntot('2') = 0.1998;
ntot('3') = 0.4994;

alpha('1','1') = 0.0; alpha('1','2') = 0.2485; alpha('1','3') = 0.30;
alpha('2','1') = 0.2485; alpha('2','2') = 0.0; alpha('2','3') = 0.3412;
alpha('3','1') = 0.30; alpha('3','2') = 0.3412; alpha('3','3') = 0.0;

tau('1','1') = 0.0; tau('1','2') = 4.93035; tau('1','3') = 1.59806;
tau('2','1') = 7.77063; tau('2','2') = 0.0; tau('2','3') = 4.18462;
tau('3','1') = 0.03509; tau('3','2') = 1.27932; tau('3','3') = 0.0;

g(i,j) = EXP(-alpha(i,j)*tau(i,j));

***************************************************************************
* gfe = Gibbs free energy
* n(i,k) = number of moles of component i in phase k
* psi(i,k) = substitution variable

VARIABLES gfe, n(i,k), psi(i,k);

***************************************************************************
* obj = objective function
* trans(i,k) = transformation constraint
* molesum(i) = material balance on component i

EQUATIONS
    obj
    trans(i,k)
    molesum(i);


obj.. gfe =e= SUM(k, SUM(i, n(i,k)*(LOG(n(i,k)) - LOG(SUM(j,n(j,k))))))
	     +SUM(k$liqphase(k), SUM(i, n(i,k)*SUM(j, g(i,j)*tau(i,j)*psi(j,k))));

trans(i,k)$liqphase(k).. psi(i,k)*SUM(j, g(j,i)*n(j,k)) - n(i,k) =e= 0.0;

molesum(i).. SUM(k, n(i,k)) =e= ntot(i);

MODEL gmin / all /;

n.lo(i,k) = 0.0000001; n.up(i,k) = ntot(i);
n.l('1','1') = 0.29949; n.l('1','2') = 0.00001;
n.l('2','1') = 0.06551; n.l('2','2') = 0.13429;
n.l('3','1') = 0.49873; n.l('3','2') = 0.00067;
psi.lo(i,k) = 0.0;
psi.l(i,k) = n.l(i,k)/SUM(j, g(j,i)*n.l(j,k));

SOLVE gmin USING nlp MINIMIZING gfe;