*--------------------------------------------------------------*
* Test Problem 8 from Chapter 8, section 5.10                  *
* Vapor-Liquid Equilibrium Model                               * 
*--------------------------------------------------------------*

Sets       
    m  number of data sets             /1*5/
    n  number of variables             /1*4/
    i  number of components            /1*2/
    j  number of Antione coefficients  /1*3/
    p  number of parameters            /1*2/;
 
Parameters 
   ze(m,n)   observed data values
   c(i,j)    Antoine coefficients
   std(n)    standard deviations;

Scalars
   Tr        reference temperature   /323.15/;

Variables 
   z(m,n)     fitted data variables
   a(p)       model parameters
   gamma(m,i) activity coefficients
   p0(m,i)    pure component vapor pressure
   cost       cost function;
   
table ze(m,n)
     1      2      3       4
1  0.30   0.591   1.00   483.8
2  0.40   0.602   1.00   493.2
3  0.50   0.612   1.00   499.9
4  0.70   0.657   1.00   501.4
5  0.90   0.814   1.00   469.7;

table c(i,j)
      1        2       3
1  18.5875  3626.55  34.29
2  16.1764  2927.17  50.22;

std('1') = 0.005;
std('2') = 0.015;
std('3') = 0.1/Tr;
std('4') = 0.75;

Equations
      obj                  objective function
      antoine(m,i)         Antoine equation
      equil1(m),equil2(m)  Equilibrium Expressions       
      vl1(m),vl2(m)        Van Laar equations;

obj.. cost =e= sum(m,sum(n,sqr((z(m,n)-ze(m,n))/std(n))));

equil1(m).. gamma(m,'1')*z(m,'1')*p0(m,'1') - z(m,'2')*z(m,'4') =e= 0;
equil2(m).. gamma(m,'2')*(1 - z(m,'1'))*p0(m,'2') - (1 - z(m,'2'))*z(m,'4') =e= 0;

antoine(m,i).. -p0(m,i) + exp( c(i,'1') - c(i,'2')/( Tr*z(m,'3') - c(i,'3'))) =e= 0;

vl1(m).. -log(gamma(m,'1')) + a('1')/z(m,'3')/
             ((1 + a('1')/a('2')*z(m,'1')/(1-z(m,'1')))**2) =e= 0;
vl2(m).. -log(gamma(m,'2')) + a('2')/z(m,'3')/
             ((1 + a('2')/a('1')*(1-z(m,'1'))/z(m,'1'))**2) =e= 0;

model
     problem /obj,antoine,equil1,equil2,vl1,vl2/;

z.lo(m,n) = ze(m,n) - 3*std(n);
z.up(m,n) = ze(m,n) + 3*std(n);

a.lo(p) = 1;
a.up(p) = 2;

z.l(m,n) = uniform(z.lo(m,n),z.up(m,n));
a.l('1') = 1.9;
a.l('2') = 1.6;
gamma.l(m,i) = 1;

solve problem using nlp minimizing cost;