*--------------------------------------------------------------*
* Test Problem 7 from Chapter 8, section 5.9                   *
* Steady-State CSTR Model 
*--------------------------------------------------------------*

Sets
    m  number of data sets      /1*10/
    n  number of variables      /1*5/
    p  number of parameters     /1*2/;
 
Parameters 
   ze(m,n)   observed data values
   std(n)    standard deviations;

Scalars
   tau       residence time          /100/
   dHr       heat of reaction        /-4180/
   rho       density                 /1.0/
   Cp        heat capacity           /4.18/
   Tr        reference temperature   /800/;

Variables 
   z(m,n)    fitted data variables
   a(p)      model parameters
   k1(m)     reaction rate
   c         cost function;
   
table ze(m,n)
        1        2        3        4        5
1    0.9871   0.8906   0.1157   547.47   663.48 
2    1.0003   0.8350   0.1380   531.77   676.04 
3    1.0039   0.8255   0.1850   512.21   684.81
4    0.9760   0.8020   0.2005   490.59   695.47
5    1.0129   0.7520   0.2420   464.67   703.69
6    1.0083   0.7193   0.2739   438.47   714.90
7    1.0075   0.6861   0.3215   408.04   726.09
8    0.9994   0.6388   0.3741   375.56   735.44
9    1.0007   0.5970   0.3926   340.26   745.70
10   0.9973   0.5580   0.4703   306.55   753.94;

std('1') = 0.01;
std('2') = 0.01;
std('3') = 0.01;
std('4') = 1.0;
std('5') = 1.0;

Equations
      obj                    objective function
      mass1(m),mass2(m)      mass balance equations 
      energy(m)              energy balance equations
      reaction(m)            reaction rate expression;     

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

mass1(m).. 1/tau*(z(m,'1') - z(m,'2')) - k1(m)*z(m,'2') =e= 0;
mass2(m).. -1/tau*z(m,'3') + k1(m)*z(m,'2') =e= 0;

energy(m) .. 1/tau*(z(m,'4') - z(m,'5')) - dHr/rho/Cp*k1(m)*z(m,'2') =e= 0;

reaction(m) .. -k1(m) + a('1')*exp(-a('2') * (Tr/z(m,'5') - 1)) =e= 0;
    
model
     problem /obj,mass1,mass2,energy,reaction/;

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

a.lo('1') = 0.0001;
a.lo('2') = 5;

a.up('1') = 0.1;
a.up('2') = 15;

z.l(m,n) = uniform(z.lo(m,n),z.up(m,n));
a.l('1') = 0.02;
a.l('2') = 12.5;

solve problem using nlp minimizing c;