sets
  i 'factories' /f1*f3/
  j 'distribution centers' /d1*d5/
  s 'individual scenarios' /lo,mid,hi/
;

parameter capacity(i) /f1 500, f2 450, f3 650/;

table demand(j,s) 'possible outcomes for demand'
      lo  mid   hi
  d1 150  160  170
  d2 100  120  135
  d3 250  270  300
  d4 300  325  350
  d5 600  700  800
;

table prob(j,s) 'probabilities for table demand'
      lo  mid   hi
  d1 .25  .50  .25
  d2 .25  .50  .25
  d3 .25  .50  .25
  d4 .30  .40  .30
  d5 .30  .40  .30
;

loop(j, abort$(abs(sum(s,prob(j,s))-1)>0.001) "probabilities don't add up");

*
* set up joint probabilities
* total scenarios =  3**5 = 243
*
set js 'scenarios for joint probabilities' /js1*js243/;
parameter jdemand(j,js) 'joint distribution: outcomes';
parameter jprob(js)   'joint distribution: probabilities';

alias (s,s1,s2,s3,s4,s5);
set current(js);
current('js1') = yes;
loop((s1,s2,s3,s4,s5),
   jdemand('d1',current) = demand('d1',s1);
   jdemand('d2',current) = demand('d2',s2);
   jdemand('d3',current) = demand('d3',s3);
   jdemand('d4',current) = demand('d4',s4);
   jdemand('d5',current) = demand('d5',s5);
   jprob(current) = prob('d1',s1)*prob('d2',s2)*prob('d3',s3)*prob('d4',s4)*prob('d5',s5);
   current(js) = current(js-1);
);
display jdemand,jprob;

table transcost(i,j) 'unit transportation cost'
       d1    d2    d3    d4    d5
  f1   2.49  5.21  3.76  4.85  2.07
  f2   1.46  2.54  1.83  1.86  4.76
  f3   3.26  3.08  2.60  3.76  4.45
;

scalar prodcost  'unit production cost' /14/;
scalar price     'sales price' /24/;
scalar wastecost 'cost of removal of overstocked products' /4/;


*-----------------------------------------------------------------------
* Form the Benders master problem
*-----------------------------------------------------------------------

set
   iter 'max Benders iterations' /iter1*iter25/
   dyniter(iter) 'dynamic subset'
;

positive variables
   ship(i,j)        'shipments'
   product(i)       'production'
   slackproduct(i)  'slack'
   received(j)      'quantity sent to market'
;
free variables
   zmaster          'objective variable of master problem'
   theta            'extra term in master obj'
;

equations
   masterobj        'master objective function'
   production(i)    'calculate production in each factory'
   receive(j)       'calculate quantity to be send to markets'
   prodcap(i)       'production capacity'
   optcut(iter)     'Benders optimality cuts'
;

parameter
   cutconst(iter)   'constants in optimality cuts'
   cutcoeff(iter,j) 'coefficients in optimality cuts'
;

masterobj..
    zmaster =e=  sum((i,j), transcost(i,j)*ship(i,j))
                 + sum(i,prodcost*product(i)) + theta;

receive(j)..       received(j) =e= sum(i, ship(i,j));
production(i)..    product(i) =e= sum(j, ship(i,j));
prodcap(i)..       product(i) + slackproduct(i) =e= capacity(i);
optcut(dyniter)..  theta =g= cutconst(dyniter) +
                              sum(j, cutcoeff(dyniter,j)*received(j));

model masterproblem /masterobj, receive, production, prodcap, optcut/;

*-----------------------------------------------------------------------
* Form the Benders' subproblem
* Notice in equation selling we use the level value received.l, i.e.
* this is a constant
*-----------------------------------------------------------------------

positive variables
   sales(j)         'sales (actually sold)'
   waste(j)         'overstocked products'
   slacksales(j)    'slack'
;
free variables
   zsub             'objective variable of sub problem'
;

equations
   subobj           'subproblem objective function'
   selling(j)       'part of received is sold'
   selmax(j)        'upperbound on sales'
;

parameter demnd(j) 'demand';


subobj..
   zsub =e= -sum(j, price*sales(j)) + sum(j, wastecost*waste(j));

selling(j)..  sales(j) + waste(j) =e= received.l(j);
selmax(j)..   sales(j) + slacksales(j) =e= demnd(j);


model subproblem /subobj,selling,selmax/;

*-------------------------------------------------------------------
* Dual subproblem
*-------------------------------------------------------------------

variables
   pi_selling(j)    'dual of equation selling'
   pi_selmax(j)     'dual of equation selmax'
;

equations
   dualobj
   dualcon(j)
;

dualobj..
   zsub =e= sum(j, received.l(j)*pi_selling(j)) + sum(j, demnd(j)*pi_selmax(j));

dualcon(j).. pi_selling(j) + pi_selmax(j) =l= -price;

pi_selling.up(j) = wastecost;
pi_selmax.up(j) = 0;

model dualsubproblem /dualobj,dualcon/;



*-------------------------------------------------------------------
* Benders algorithm
*-------------------------------------------------------------------


*
* step 1: solve master without cuts
*
dyniter(iter) = NO;
cutconst(iter) = 0;
cutcoeff(iter,j) = 0;
theta.fx = 0;
solve masterproblem minimizing zmaster using lp;



*
* repair bounds
*
theta.lo = -INF;
theta.up = INF;

scalar lowerbound /-INF/;
scalar upperbound /INF/;
parameter objsub(js);
scalar objmaster;
objmaster = zmaster.l;


option limrow = 0;
option limcol = 0;
*subproblem.solprint = 0;
masterproblem.solprint = 0;

parameter pselling(j),pselmax(j);

loop(iter,

*
* solve subproblems
*
   dyniter(iter) = yes;
   loop(js,
       loop(j,
           if (received.l(j)>jdemand(j,js),
               pselling(j) = wastecost;
               pselmax(j) = -price-wastecost;
           else
               pselling(j) = -price;
               pselmax(j) = 0;
           );
       );
       objsub(js) = sum(j, pselling(j)*received.l(j)) + sum(j, pselmax(j)*jdemand(j,js));
       cutconst(iter) = cutconst(iter) - jprob(js)*sum(j,(-pselmax(j))*jdemand(j,js));
       cutcoeff(iter,j) = cutcoeff(iter,j) - jprob(js)*(-pselling(j));
   );

   upperbound = min(upperbound, objmaster + sum(js, jprob(js)*objsub(js)));

*
* convergence test
*
   display lowerbound,upperbound;
   abort$( (upperbound-lowerbound) < 0.0001*(1+abs(lowerbound)) ) "Converged";

*
* solve masterproblem
*
   solve masterproblem minimizing zmaster using lp;
   display ship.l;

   lowerbound = zmaster.l;
   objmaster = zmaster.l-theta.l;

);