# adapted from LOTSIZELIB's pp08a instance
# (C) Joao Pedro PEDROSO

param T;	# number of periods (months)
param NProd;
param NBlck;
set Prod := {1..NProd};	# products
set Blck := {1..NBlck};	# machines (blocks)
set Compat {p in Prod} within {Blck} default {Blck};	# machines that can be used to produce p


param SetUpCost {Prod};

param HoldingCost {Prod, 1..T} default 1.;
param BacklogCost {Prod, 1..T} default 2.;

param Demand {Prod, 1..T} default 0;	# million pieces (for each month)

param L {1..T} >= 0;	# -- total capacity 

param DSum {p in Prod} := sum {t in 1..T} Demand[p,t];
param Cap {p in Prod, b in Compat[p], t in 1..T} :=  min(DSum[p], L[t]);


param FirstInv {Prod} default 0;	# million pieces
param LastInv {Prod} default 0;		# million pieces


var y{p in Prod, b in Compat[p], 1..T} binary; 	# 1 if p is produced in b on period t, 0 oth.
var x{p in Prod, b in Compat[p], 1..T} >=0;	# quantity of p produced on block b, period t
var h{Prod,0..T} >=0;	# inventory, million pieces
var g{Prod,0..T} >=0;	# backlog, million pieces

var fcost;	# fixed, setup cost
var hcost;	# holding cost
var gcost;	# backlog cost

#
# objective
#
minimize TCost:
  fcost + hcost + gcost;

subject to FCost:
  fcost = sum{p in Prod, b in Compat[p], t in 1..T} SetUpCost[p]*y[p,b,t] ;
subject to HCost:
  hcost = sum{p in Prod, t in 1..T} HoldingCost[p,t]*h[p,t];
subject to GCost:
  gcost = sum{p in Prod, t in 1..T} BacklogCost[p,t]*g[p,t] +
	sum{p in Prod} 10000.*g[p,T];	# !!! for avoid final backlog > 0


#
# machine setup
#
subject to ResourceUpperBound {b in Blck, t in 1..T}:
  sum {p in Prod : b in Compat[p]} x[p,b,t] <= L[t];

subject to ProductionSetupConnect { p in Prod, b in Compat[p], t in 1..T }:
  x[p,b,t] <= Cap[p,b,t] * y[p,b,t];


# material conservation
subject to InitialStock{p in Prod}:
  h[p,0] = FirstInv[p];
subject to FinalStock{p in Prod}:
  h[p,T] = LastInv[p];
subject to InitialBacklog{p in Prod}:
  g[p,0] = 0;
### subject to FinalBacklog{p in Prod}:
###   g[p,T] = 0;

subject to FlowConservation{p in Prod,t in 1..T}:
  (h[p,t-1] - g[p,t-1]) + sum{b in Compat[p]} x[p,b,t]  =  
	Demand[p,t] + (h[p,t] - g[p,t])
;

end;