(!*****************************************************
* File: scheduling.mos                                *
* Reading data from: data.dat                         *
* Ref: Ch 8- Planning and Scheduling in Supply chains *
*     by Mike Pinedo                                  *
* (c) 2004 Dash Optimization Inc                      *
*     author: Nitin Verma                             *
*     e-mail: nitin.verma@dashoptimization.com        *
*****************************************************!)
model "Supply Chain Scheduling"
uses "mmxprs" ! Use Xpress-Optimizer

parameters
	DatFile='data.dat'
end-parameters

declarations
	Weeks = 1..4 !indexed by t
	ProductFamily = 1..2 !indexed by j
	Factories = 1..2 !indexed by i
	Stages = 2..3 !indexed by l
	
	D: array(Stages,ProductFamily, Weeks) of real
	c_p, p_: array(Factories, ProductFamily) of real
	c_m20, c_m03: array(Factories) of real
	MaxHrPerWeek,h,c_m023,si,bigM,tau:real
	w__,w___: array(ProductFamily) of real
	UB,LB:dynamic array(Factories, Stages, ProductFamily) of real
	
	x: array(Factories, ProductFamily, Weeks) of mpvar
	y: array(Factories, Stages, ProductFamily, Weeks) of mpvar
	z,v2,delta_q2,delta_v2,delta_v3: array(ProductFamily,Weeks) of mpvar
	q2,v3: array(ProductFamily, {0}+Weeks) of mpvar
	delta_y: dynamic array (Factories, Stages, ProductFamily,Weeks) of mpvar
	scale=1
end-declarations

initializations from DatFile
	D c_p p_ h c_m20 c_m03 c_m023 w__ w___ tau si bigM UB LB
end-initializations

! Objective
Cost:=
	sum(t in Weeks, j in ProductFamily, i in Factories) c_p(i,j)*x(i,j,t) + 
	sum(t in Weeks, j in ProductFamily, i in Factories) c_m20(i)*y(i,2,j,t) + 
	sum(t in Weeks, j in ProductFamily, i in Factories) c_m03(i)*y(i,3,j,t) + 
	sum(t in Weeks, j in ProductFamily) c_m023*z(j,t) + 
	sum(t in Weeks, j in ProductFamily) h*q2(j,t) + 
	sum(t in 1..3, j in ProductFamily) w__(j)*v2(j,t) + 
	sum(t in 1..3, j in ProductFamily) w___(j)*v3(j,t) + 
	sum(j in ProductFamily) si*v2(j,4) + 
	sum(j in ProductFamily) si*v3(j,4)

! Weekly production capacity constraints       
forall(i in Factories, t in Weeks) 
	ProductCapacity(i,t):=sum(j in ProductFamily) p_(i,j)*x(i,j,t) <= 24*7*tau/scale

forall(i in Factories,l in Stages, j in ProductFamily | exists(UB(i,l,j)), t in Weeks) 
	TransPortUpperBound(i,l,j,t):=y(i,l,j,t) <= UB(i,l,j)/scale

forall(i in Factories,l in Stages, j in ProductFamily | exists(LB(i,l,j)), t in Weeks) do
	create(delta_y(i,l,j,t))
	delta_y(i,l,j,t) is_binary
	TransPortLowerBound(i,l,j,t):=y(i,l,j,t) >= LB(i,l,j)*delta_y(i,l,j,t)/scale
	TransPortUpperBound(i,l,j,t):=y(i,l,j,t) <= if(exists(UB(i,l,j)),UB(i,l,j),bigM)*delta_y(i,l,j,t)/scale
end-do
 
! Transportation constraints
forall(t in Weeks, j in ProductFamily, i in Factories) 
	Transport(t,j,i):=sum(l in Stages) y(i,l,j,t) = x(i,j,t)

forall(t in 1..3, j in ProductFamily)
	TransportLimitPlant(t,j):=sum(i in Factories) y(i,3,j,t) + z(j,t) <= D(3,j,t+1)/scale + v3(j,t)

forall(t in Weeks, j in ProductFamily)
	TransportLimitDC(t,j):=z(j,t) <= q2(j,t-1)+ if( t=1, 0, y(1,2,j,t-1) + y(2,2,j,t-1))
 
! Storage constraints
forall(j in ProductFamily, t in Weeks) do
	StorageLimitDC(j,t):=q2(j,t) <= bigM*delta_q2(j,t)/scale
	StorageMaxDC1(j,t):=q2(j,t) - q2(j,t-1) + 
		if (t=1, 0, - y(1,2,j,t-1) - y(2,2,j,t-1) + v2(j,t-1)) +D(2,j,t)/scale + z(j,t) <= bigM*(1-delta_q2(j,t))/scale
	StorageMaxDC2(j,t):=q2(j,t) - q2(j,t-1) + 
		if (t=1, 0, - y(1,2,j,t-1) - y(2,2,j,t-1) + v2(j,t-1)) +D(2,j,t)/scale + z(j,t) >=0
 end-do

forall(j in ProductFamily, t in Weeks) do
	TardyLimitDC(j,t):=v2(j,t) <= bigM*delta_v2(j,t)/scale
	TardyMaxDC1(j,t):=v2(j,t) - D(2,j,t)/scale + q2(j,t-1) +
		if (t=1, 0, - v2(j,t-1) - z(j,t) + y(1,2,j,t-1) + y(2,2,j,t-1)) <= bigM*(1-delta_v2(j,t))/scale
	TardyMaxDC2(j,t):=v2(j,t) - D(2,j,t)/scale + q2(j,t-1) +
		if (t=1, 0, - v2(j,t-1) - z(j,t) + y(1,2,j,t-1) + y(2,2,j,t-1))>=0
end-do

forall(j in ProductFamily, t in Weeks) do
	if t=1 then
		TardyLimitCust(j,t):=v3(j,t) = maxlist(0,D(3,j,t))/scale
	else
		TardyLimitCust(j,t):=v3(j,t) <= bigM*delta_v3(j,t)/scale
		TardyMaxCust1(j,t):=v3(j,t) - D(3,j,t)/scale- v3(j,t-1) + z(j,t-1) + y(1,3,j,t-1) + y(2,3,j,t-1) <= bigM*(1-delta_v3(j,t))/scale
		TardyMaxCust2(j,t):=v3(j,t) - D(3,j,t)/scale- v3(j,t-1) + z(j,t-1) + y(1,3,j,t-1) + y(2,3,j,t-1) >=0
	end-if
end-do

! Boundary conditions
forall(j in ProductFamily) do
	InitTardy(j):=v3(j,0) = 0
	InitTransport(j):=q2(j,0) = 0
end-do

! Binary variables
forall(t in Weeks, j in ProductFamily) do
	delta_q2(j,t) is_binary
	delta_v2(j,t) is_binary
	delta_v3(j,t) is_binary
end-do

! Solve the problem
setparam('xprs_loadnames',true)
minimize(Cost)

! Solution printing
writeln("Total cost: ", getobjval)
forall(i in Factories, j in ProductFamily, t in Weeks)
	writeln("x(",i,",",j,",",t,") = ", getsol(x(i,j,t)))
forall(i in Factories, l in Stages, j in ProductFamily, t in Weeks)
	writeln("y(",i,",",l,",",j,",",t,") = ", getsol(y(i,l,j,t)))

end-model