$offsymxref
$offsymlist
options limrow = 0,
        limcol = 0,
        iterlim = 1e9,
        reslim = 3600,
        sysout = off,
        optcr = 0   ;

*option minlp = alphaecp;
*option minlp = baron;
*Option MINLP = CoinBonmin;
*option minlp = dicopt;
option minlp = sbb;

sets i products /A, B, C, D, E/
     j stages   /1*6/
     k binary digit /1*2/;

parameters H        horizon timie (h)
           alpha(j) cost coefficient ($)
           beta(j)  cost coefficient
           VL(j)    lower bound of volume  (L)
           VU(j)    upper bound of volume  (L)
           NU(j)    maximum number of parallel units ;

H = 6000;
alpha(j) = 250;
beta(j) = 0.6;
VL(j) = 300; VU(j) = 3000; NU(j) = 4 ;

parameter Q(i) production rate (kg) /
         A       250000
         B       150000
         C       180000
         D       160000
         E       120000        /;

table S(i,j) size factor for product i in stage j (L per kg)
       1         2       3       4       5       6
A     7.9       2.0     5.2     4.9     6.1     4.2
B     0.7       0.8     0.9     3.4     2.1     2.5
C     0.7       2.6     1.6     3.6     3.2     2.9
D     4.7       2.3     1.6     2.7     1.2     2.5
E     1.2       3.6     2.4     4.5     1.6     2.1

table t(i,j) processing time for product i in stage j (L per kg)
       1          2       3       4       5       6
A     6.4        4.7     8.3     3.9     2.1     1.2
B     6.8        6.4     6.5     4.4     2.3     3.2
C     1.0        6.3     5.4     11.9    5.7     6.2
D     3.2        3.0     3.5     3.3     2.8     3.4
E     2.1        2.5     4.2     3.6     3.7     2.2


positive variables N(j), V(j), B(i), TL(i);
binary variables y(k,j);
variable z;
equations obj, c1, c2, c3, c4;

V.lo(j) = VL(j);   V.up(j) = VU(j);
N.lo(j) = 1;   N.up(j) = NU(j) ;
TL.lo(i) = smax(j, t(i,j)/NU(j) );  TL.up(i) = smax(j, t(i,j));
B.lo(i) = Q(i)/ H * smax(j, t(i,j)/NU(j) );     B.up(i) = min( Q(i), smin(j, VU(j)/S(i,j)) );

N.l(j) = N.lo(j); V.l(j) = V.lo(j); B.l(i) = B.lo(i); TL.l(i) = TL.lo(i);

obj..       z =e= sum(j, alpha(j) * N(j) * ( V(j)**beta(j) )  );

c1(i,j)..   V(j) =g= S(i,j) * B(i);
c2(i,j)..   N(j) * TL(i) =g= t(i,j);
c3..        sum(i, Q(i)*TL(i)/B(i) ) =l= H;
c4(j)..     N(j) =e= 1 + sum(k, power(2, ord(k)-1) * Y(k,j) );

model P1 /all/;
P1.workspace = 1024;
P1.workfactor = 1;

solve P1 using minlp minimizing z;