$title A Production Planning Example (PRODPLAN, SEQ=356)

$onText
This uncapacitated lot-sizing problem finds the least cost production plan
meeting demand requirements. There are costs given for producing, stocking,
and setting up the machines.

Four solving approaches are presented:
1) Solving the original model as a MIP
2) Solving a tight reformulation as an RMIP
3) Solving a tight reforumulation without stock as an RMIP
4) Solving the original model as an RMIP using a separation algorithm

Pochet, Y, and Wolsey, L A, Production Planning by Mixed Integer
Programming (Springer Series in Operations Research and Financial
Engineering). Springer-Verlag New York, Inc., 2006.

Keywords: mixed integer linear programming, relaxed mixed integer linear
          programming, production planning
$offText

*1) Initial tiny formulation
Set
   t       'time periods'  / t1*t8 /
   ut(t,t) 'upper triangle';

Alias (t,tt,k);

Parameter
   DEMAND(T) 'demand per period'          / (t1,t2) 400, (t3,t4) 800, (t5*t8) 1200 /
   SETUPCOST 'setup cost per period'      / 5000 /
   PRODCOST  'production cost per period' /  100 /
   INVCOST   'production cost per period' /    5 /
   STOCKINI  'production cost per period' /  200 /
   BigM(T)   'max production - BigM';

*We assume that the initial stock is lower equal than the demand in the first period
abort$(Demand('t1') < STOCKINI) 'Initial stock is too large';

ut(k,t) = ord(k) <= ord(t);

BigM(t) = sum(k$(ord(k) >= ord(t)), DEMAND(k) - STOCKINI$(ord(t) = 1));

display ut, BigM;

Variable
   s(t) 'inventory in period t'
   x(t) 'production in period t'
   y(t) 'setup in period t'
   cost;

Binary   Variable y;
Positive Variable s, x;

Equation
   Balance(t)    'stock balance'
   Production(t) 'production set-up'
   Mincost       'objective function';

Mincost..
   cost =e= sum(t, ifthen(ord(t) < card(t),INVCOST,INVCOST/2)*s(t)) + sum(t, SETUPCOST*y(t) + PRODCOST*x(t));

Production(t)..
   x(t) =L= BigM(t)*y(t);

Balance(t)..
   STOCKINI$(ord(t) = 1) + s(t-1) + x(t) =e= DEMAND(t) + s(t);

Model tiny / Mincost, Production, Balance /;

tiny.optCr=0;

solve tiny minimizing cost using mip;

*2) Multi-commodity formulation (tight reformulation)
Variable
   smc(t,t) 'inventory entered in period i for period t'
   xmc(t,t) 'production in period i for demand in t';

Positive Variable smc, xmc;

Equation
   Balancemc(t,t)    'stock balance'
   Productionmc(t,t) 'production set-up'
   Mincostmc         'objective function';

Mincostmc..
   cost =e= sum(ut, PRODCOST*xmc(ut) + INVCOST*smc(ut)) +  sum(t, SETUPCOST*y(t));

Balancemc(ut(k,t))..
   STOCKINI$(ord(t) = 1) + smc(k-1,t) + xmc(k,t) =e= smc(k,t) + diag(k,t)*DEMAND(t);

Productionmc(ut(k,t))..
   xmc(k,t) =l= (DEMAND(t) - STOCKINI$(ord(t) = 1))*y(k);

Model tinymc / Mincostmc, Balancemc, Productionmc /;

solve tinymc minimizing cost using rmip;

*3) Multi-commodity formulation without stock (tight reformulation)
Parameter dist(t,t) 'distance between time periods';
dist(ut(k,t)) = ord(t) - ord(k);

display dist;

Equation
   Demandmcws(t) 'demand satisfaction'
   Mincostmcws   'objective function';

Mincostmcws..
   cost =e= sum(ut, PRODCOST*xmc(ut) + INVCOST*dist(ut)*xmc(ut)) + sum(t, SETUPCOST*y(t));

Demandmcws(t)..
   sum(ut(k,t), xmc(k,t)) =g= DEMAND(t) - STOCKINI$(ord(t) = 1);

Model tinymcws / Mincostmcws, Demandmcws, Productionmc /;

solve tinymcws minimizing cost using rmip;

*4) Separation Algorithm
Set
   j           'iterations'  /j1*j10/
   n(j,t)      'set of cuts'
   Scon(j,t,t) 'set of violated constraints';

n(j,t)      = no;
Scon(j,t,t) = no;

Alias (t,l), (j,jj);

Parameter
   D(t,t)    'accumulated demand'
   left(t,t) 'left side of cut';

D(ut(t,k)) = sum[tt$(ord(tt) <= ord(k) and ord(tt) >= ord(t)), DEMAND(tt)];

Equation cuts(j,t) 'cuts for the RMIP (complete linear description)';

cuts(n(jj,t)).. sum(Scon(jj,t,k), x(k) - D(k,t)*y(k)) =l= s(t);

Model tinycuts / tiny, cuts /;

Scalar
   more    / 1    /
   epsilon / 1e-6 /;

*If STOCKINI < DEMAND(t1) there has to be production in the first period
y.fx('t1') = 1;

loop(j$more,
   solve tinycuts using rmip min cost;
   option limCol = 0, limRow = 0, solPrint = silent;

*  Store the left hand side of potential cuts
   left(ut(tt,l)) = x.l(tt)-d(tt,l)*y.l(tt);

*  Use only those LHS which are greater zero
   Scon(j,l,tt) = left(tt,l) > epsilon;

*  If the sum of those is greater than the inventory level: violation found
*  Add this cut to the model
   n(j,l) = sum[Scon(j,l,tt), left(tt,l)] - epsilon > s.l(l);

*  Proceed if at least one cut was added during this iteration
   more = sum(n(j,l), yes);
);

put_Utility$(not more) 'log' /
  '>>>>Integer solution found. A total of 'sum(n(j,t),1):0:0' cuts were added.';