#*************************************************************
# Test Problem 5 from Chapter 11, section 4.6  
# Bifunctional Catalyst Blend Problem 
#*************************************************************

OPTION {{
  FTOL = 1e-6;
  OTOL = 1e-6;
  NLP = "SNOPT";
  RANDOMX;
}}

DECLARATION {{
  INDEX {i};
  SET CP = |0:9|;   # Number of control parameters
  SET J  = |1:7|;   # Number of states
  SET I  = |1:10|;  # Number of reaction rates
  SET K  = |1:4|;   # Number of constants in reaction rate expression

  XVAR {z7f,    # Final Time Value of state 7
        w(CP) # control parameters
       };

  ZVAR {z(J),u,k(I)};
  ISPE {z(J)};
  STP z(J) = {1,0,0,0,0,0,0};

  TIME{0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1};

# Constants in the reaction rate expression  
  PARA c(I,K) = { 0.2918487e-2, -0.8045787e-2, 0.6749947e-2,-0.1416647e-2,
                    0.9509977e1,  -0.3500994e2,  0.4283329e2, -0.1733333e2,
                    0.2682093e2,  -0.9556079e2,  0.1130398e3, -0.4429997e2,
                    0.2087241e3,  -0.7198052e3,  0.8277466e3, -0.3166655e3,
                    0.1350005e1,  -0.6850027e1,  0.1216671e2, -0.6666689e1,
                    0.1921995e-1, -0.7945320e-1, 0.1105666,   -0.5033333e-1,
                    0.1323596,    -0.4696255,    0.5539323,   -0.2166664,
                    0.7339981e1,  -0.2527328e2,  0.2993329e2, -0.1199999e2,
                   -0.3950534,     0.1679353e1, -0.1777829e1,  0.4974987,
                   -0.2504665e-4,  0.1005854e-1,-0.1986696e-1, 0.9833470e-2};

# Reactor Volume
  PARA vr = 2000;

# Variable Bounds
  LBDS w(CP) = < i E CP | 0.6>;
  UBDS w(CP) = < i E CP | 0.9>;
  
  LBDS z7f = 0;
  UBDS z7f = 1;

}}

MODEL {{
  MIN:  -z7f;

  dae1: z'(1) =e= -k(1)*z(1);
  dae2: z'(2) =e=  k(1)*z(1) - (k(2) + k(3))*z(2) + k(4)*z(5); 
  dae3: z'(3) =e=  k(2)*z(2);
  dae4: z'(4) =e= -k(6)*z(4) + k(5)*z(5);
  dae5: z'(5) =e=  k(3)*z(2) + k(6)*z(4) - (k(4) + k(5) + k(8) + k(9))*z(5)
                 + k(7)*z(6) + k(10)*z(7);
  dae6: z'(6) =e=  k(8)*z(5) - k(7)*z(6);
  dae7: z'(7) =e=  k(9)*z(5) - k(10)*z(7);

  reaction(i E I): k(i) =e= vr*(c(i,1) + c(i,2)*u + c(i,3)*u^2 + c(i,4)*u^3);

  control: u =e= interv[i E CP| w(i)];

  endtime[10]: z7f - z(7) =e= 0; 

}}