param N;	# number of input items
param K;	# number of machines

param D;	# demand or consumption rate for the end product (units/period)
param r;	# inventory carrying cost rate ($/$/period)

param A{1..N};	# ordering cost for each input item ($/order)
param b{1..N};	# purchase cost for each input item ($/unit)

param S{1..K};	# setup cost per batch for each machine ($/setup)
param c{1..K};	# total production cost of an inventory item after processing on a machine ($/unit)
param P{1..K};	# production rate for each machine (units/period)

var L{1..N} >= 1, := 1, integer;
var R{1..N} >= 1, := 1, integer;

var Q >= 0, := 1;	# production lot size (units)

minimize Total_Cost:
	sum {i in 1..N} D*A[i]/(L[i]+1/R[i]-1)/Q + sum {j in 1..K} D*S[j]/Q + 
	D*Q*r*(sum {i in 1..N} (b[i]/R[i]/P[1] + (L[i]-1)/D) + 
	sum {j in 1..K-1} c[j]*abs(1/P[j] - 1/P[j+1]) + c[K]*(1/D-1/P[K]))/2;

subject to complementarity{i in 1..N}:
	0 <= (L[i]-1) complements (R[i]-1) >= 0;

#data "data/cell-8-4-15.dat";
#let Q := sqrt((2*A[1]+2*sum{j in 1..K}S[j])/(r*b[1]/P[1] + 
#		r*(sum {j in 1..K-1} c[j]*abs(1/P[j] - 1/P[j+1]) + c[K]*(1/D-1/P[K]))));
#display Total_Cost;
#solve;
#display L, R;
#display Q;

#commands cell.cmd;

# Explicit formula for M when N = 1 and M < 1
#display 
#  sqrt( (b[1]*sum {j in 1..K} S[j]) / (A[1]*P[1]*(sum {j in 1..K-1} c[j]*abs(1/P[j] - 1/P[j+1]) + c[K]*(1/D-1/P[K]))));

# Explicit formula for M when N = 1 and M > 1
# display 1/sqrt( (b[1]*sum {j in 1..K} S[j]) / (A[1]*D*(b[1]*(1/P[1]-1/D) + sum {j in 1..K-1} c[j]*abs(1/P[j] - 1/P[j+1]) + c[K]*(1/D-1/P[K]))));