*****************************************************************
*** Chapter 6
*** Test Problem 12
***
*** n-Butyl-Acetate and Water -- Gibbs energy minimization(ASOG)
*****************************************************************

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

*****************************************************************
* P = pressure (atm)
* T = temperature (K)
* liqphase(k) = is phase k a liquid phase?
* ntot(i) = number of moles of component i in feed
* v(m,i) = group-component matrix
* am(m,l), an(m,l), a(m,l) = group binary interaction parameters
* nu(i), vhat(m,i) = pure component and calculated parameters

PARAMETERS
P, T, liqphase(k), ntot(i), z(i), v(m,i), am(m,l), an(m,l), nu(i), a(m,l), vhat(m,i),
vs(i), lambda(i,l), vi(i);

P = 1.0;
T = 298.0;

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

ntot('1') = 0.5;
ntot('2') = 0.5;

z('1') = 0.5;
z('2') = 0.5;

v('1','1') = 5; v('1','2') = 0;
v('2','1') = 3; v('2','2') = 0;
v('3','1') = 0; v('3','2') = 1.6;

am('1','1') = 0.0; am('1','2') = -15.2623; am('1','3') = -0.2727;
am('2','1') = -0.3699; am('2','2') = 0.0; am('2','3') = -2.5548;
am('3','1') = 0.5045; am('3','2') = -2.4686; am('3','3') = 0.0;

an('1','1') = 0.0; an('1','2') = 515.0; an('1','3') = -277.3;
an('2','1') = 162.6; an('2','2') = 0.0; an('2','3') = 659.9;
an('3','1') = -2382.3; an('3','2') = 565.7; an('3','3') = 0.0;

nu('1') = 8;
nu('2') = 1;

a(l,m) = EXP(am(l,m)+an(l,m)/T);
vhat(l,i) = SUM(m, v(m,i)*a(l,m));
vs(i) = SUM(l, v(l,i));
lambda(i,l) = EXP(1-LOG(vhat(l,i)/vs(i))+SUM(m, -v(m,i)*a(m,l)/vhat(m,i)));
vi(i) = SUM(l, v(l,i)*lambda(i,l));
vi('1') = 2.009819;
vi('2') = -0.752006;

DISPLAY vhat;
DISPLAY a;
DISPLAY vs;
DISPLAY vi;

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

VARIABLES gfe, n(i,k);

***************************************************************************
* obj = objective function
* molesum(i) = material balance on component i

EQUATIONS
    obj
    molesum(i);


obj.. gfe =e= SUM(k, SUM(i, n(i,k)*(LOG(nu(i)) - vi(i)))
		    +SUM(i, n(i,k)*LOG(n(i,k)/SUM(j,nu(j)*n(j,k))))
		    +SUM(i, vs(i)*n(i,k))*LOG(SUM(j, vs(j)*n(j,k)))
           +SUM(i, SUM(m, n(i,k)*v(m,i)*LOG(n(i,k)/SUM(j,vhat(m,j)*n(j,k))))) )
	     +SUM(k,SUM(i, -vs(i)*n(i,k)*LOG(n(i,k))));

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.4994; n.l('1','2') = 0.0006;
n.l('2','1') = 0.1179; n.l('2','2') = 0.3821;

SOLVE gmin USING nlp MINIMIZING gfe;