* chemeq.mod    OLR2-MY-38-50
* Original AMPL coding by Elena Bobrovnikova (summer 1996 at Bell Labs).

* Chemical equilibrium problem

* Ref.: D. M. Himmelblau, Applied Nonlinear Programming,
* McGraw-Hill Book Company, New York, 1972.  Problem 6.

* Number of variables:  38 (40 before presolve)
* Number of constraints:  50 (56 before presolve)
* Objective partially separable
* Objective convex
* Linear constraints

* Since constraints 14 and 16 imply x[1,5]=x[2,5]=x[3,5]=x[1,7]=x[2,7]=0
* these variables were removed from the original formulation of the problem

Set J  /1*18/;
Set K  /1*7 /;
Set I  /1*16/;

Alias(j,j1);
Set PAIRS(J,K);

Parameter b(I)    ;  b(I)=0;
Parameter c(J,K)  ;
Parameter E(J,K,I);
Parameter xlb     ;  xlb= 0;   xlb= 1.5e-6;

Table LINK(J,K)
              1      2     3     4     5     6     7
       1      1      1     1     1           1      
       2      1      1     1     1           1      
       3      1      1     1     1                  
       4      1      1     1                        
       5             1     1                        
       6             1     1                        
       7             1     1                        
       8             1     1                        
       9             1     1                        
      10             1     1                        
      11             1     1                        
      12             1     1                        
      13             1     1                        
      14                   1                        
      15                   1                        
      16                   1                        
      17                   1                        
      18                   1                          ;

PAIRS(J,K)=yes$(LINK(J,K)=1);

parameter      b(I);  b(i)=0;
      b('1') =  0.6529581;
      b('2') =  0.281941 ;
      b('3') =  3.705233 ;
      b('4') = 47.00022  ;
      b('5') = 47.02972  ;
      b('6') =  0.08005  ;
      b('7') =  0.08813  ;
      b('8') =  0.04829  ;
      b('9') =  0.0155   ;
     b('10') =  0.0211275;
     b('11') =  0.0022725;

table   c(J,k)
              1        2        3         4        5         6        7
        1            -10.94    10.45   -15.639                        
        2    -7.69                                         11.959     
        3   -11.52             -0.5     21.81                         
        4   -36.6                                                     
        5                                                             
        6                                                             
        7                       2.2435                                
        8              2.5966                                         
        9            -39.39   -39.39                                  
       10            -21.35   -21.49                                  
       11            -32.84   -32.84                                  
       12              6.26     6.12                                  
       13                                                             
       14                                                             
       15                      -1.9028                                
       16                      -2.8889                                
       17                      -3.3622                                
       18                      -7.4854                                   ;

parameter   E(J,K,I);  E(J,K,I)=0;

        E('1','1','1')   =  1  ;
        E('2','1','2')   =  1  ;
        E('3','1','3')   =  1  ;
        E('4','1','4')   =  1  ;
        E('4','1','5')   =  1  ;
        E('1','2','1')   =  1  ;
        E('2','2','2')   =  1  ;
        E('3','2','3')   =  1  ;
        E('4','2','4')   =  1  ;
        E('4','2','12')  =  1  ;
        E('5','2','5')   =  1  ;
        E('5','2','12')  = -1  ;
        E('6','2','6')   =  1  ;
        E('6','2','12')  = -1  ;
        E('7','2','7')   =  1  ;
        E('7','2','12')  =  1  ;
        E('8','2','8')   =  1  ;
        E('8','2','12')  =  1  ;
        E('9','2','4')   =  1  ;
        E('9','2','5')   =  1  ;
        E('10','2','2')  =  1  ;
        E('10','2','5')  =  1  ;
        E('10','2','12') = -1  ;
        E('11','2','2')  =  1  ;
        E('11','2','4')  =  1  ;
        E('11','2','5')  =  1  ;
        E('12','2','2')  =  1  ;
        E('12','2','4')  = -1  ;
        E('12','2','5')  =  1  ;
        E('12','2','12') = -2  ;
        E('13','2','9')  =  1  ;
        E('13','2','12') = -1  ;
        E('1','3','1')   =  1  ;
        E('2','3','2')   =  1  ;
        E('3','3','3')   =  1  ;
        E('4','3','4')   =  1  ;
        E('5','3','5')   =  1  ;
        E('6','3','6')   =  1  ;
        E('7','3','7')   =  1  ;
        E('8','3','8')   =  1  ;
        E('9','3','4')   =  1  ;
        E('9','3','5')   =  1  ;
        E('10','3','2')  =  1  ;
        E('10','3','5')  =  1  ;
        E('11','3','2')  =  1  ;
        E('11','3','4')  = -1  ;
        E('11','3','5')  =  1  ;
        E('13','3','10') =  1  ;
        E('14','3','11') =  1  ;
        E('14','3','12') = -4  ;
        E('15','3','1')  =  1  ;
        E('15','3','11') =  1  ;
        E('15','3','12') = -3  ;
        E('15','3','13') = -1  ;
        E('16','3','1')  =  2  ;
        E('16','3','11') =  1  ;
        E('16','3','12') = -2  ;
        E('16','3','13') = -2  ;
        E('17','3','1')  =  3  ;
        E('17','3','11') =  1  ;
        E('17','3','12') = -1  ;
        E('17','3','13') = -3  ;
        E('18','3','1')  =  4  ;
        E('18','3','11') =  1  ;
        E('18','3','13') = -4  ;
        E('1','4','4')   =  1  ;
        E('1','4','13')  =  1  ;
        E('2','4','13')  =  1  ;
        E('3','4','4')   = -1  ;
        E('3','4','13')  =  1  ;
        E('3','4','15')  = -4  ;
        E('1','6','15')  =  1  ;
        E('2','6','2')   =  1  ;
        E('2','6','4')   = -1  ;
        E('2','6','15')  =  1  ;

Variable  x(J,K),
          energy;

Equation  h(i)  ,
          defobj;

defobj.. energy=e= sum{(j,k)$PAIRS(j,k),(x[j,k]*(c[j,k]  +
                   log(x[j,k]/sum{(j1)$PAIRS(j1,k),x[j1,k]})))} ;

h(i)..   sum{(j,k)$PAIRS(j,k),(E[j,k,i]*x[j,k])}-b[i]=e=0;

 x.lo(J,K)=xlb;

**********************************************************
model chemed /all/;
solve chemed using nlp minimaze energy;

file demo /output/;  put demo;
Put "energy = ", energy.l:20:10;  Put /;
Loop(j,put j.tl:3 Loop(k, Put x.l[j,k]:11:5);Put /;);