(!*******************************************************
* Mosel Example Problems *
* ====================== *
* *
* file els.mos *
* ```````````` *
* Example for the use of the Mosel language *
* (Economic lot sizing, ELS, problem, solved by *
* adding(l,S)-inequalities in one or several rounds *
* looping over the root node or at tree nodes) *
* *
* ELS considers production planning over a horizon *
* of T periods. In period t, t=1,...,T, there is a *
* given demand DEMAND[p,t] that must be satisfied by *
* production produce[p,t] in period t and by *
* inventory carried over from previous periods. *
* A set-up cost SETUPCOST[t] is associated with *
* production in period t, and the total production *
* capacity per period is limited. The unit production *
* cost in period t is PRODCOST[p,t]. There is no *
* inventory or stock-holding cost. *
* The model implements a configurable cutting plane *
* algorithm for this problem. *
* *
* (c) 2008 Fair Isaac Corporation *
* author: S. Heipcke, 2001, rev. 2006, rev. 2010 *
*******************************************************!)
model Els ! Start a new model
uses "mmxprs","mmsystem" ! Load the optimizer library
parameters
ALG = 6 ! Algorithm choice [default settings: 0]
CUTDEPTH = 10 ! Maximum tree depth for cut generation
EPS = 1e-6 ! Zero tolerance
end-parameters
forward procedure lp_cut_gen ! Declare a procedure that is
! defined later
forward public function cb_node:boolean
forward procedure tree_cut_gen
declarations
TIMES = 1..15 ! Range of time
PRODUCTS = 1..4 ! Set of products
DEMAND: array(PRODUCTS,TIMES) of integer ! Demand per period
SETUPCOST: array(TIMES) of integer ! Setup cost per period
PRODCOST: array(PRODUCTS,TIMES) of integer ! Production cost per period
CAP: array(TIMES) of integer ! Production capacity per period
D: array(PRODUCTS,TIMES,TIMES) of integer ! Total demand in periods t1 - t2
produce: array(PRODUCTS,TIMES) of mpvar ! Production in period t
setup: array(PRODUCTS,TIMES) of mpvar ! Setup in period t
solprod: array(PRODUCTS,TIMES) of real ! Sol. values for var.s produce
solsetup: array(PRODUCTS,TIMES) of real ! Sol. values for var.s setup
starttime: real
end-declarations
initializations from "Data/els.dat"
DEMAND SETUPCOST PRODCOST CAP
end-initializations
forall(p in PRODUCTS,s,t in TIMES) D(p,s,t):= sum(k in s..t) DEMAND(p,k)
! Objective: minimize total cost
MinCost:= sum(t in TIMES) (SETUPCOST(t) * sum(p in PRODUCTS) setup(p,t) +
sum(p in PRODUCTS) PRODCOST(p,t) * produce(p,t) )
! Satisfy the total demand
forall(p in PRODUCTS,t in TIMES)
Dem(p,t):= sum(s in 1..t) produce(p,s) >= sum (s in 1..t) DEMAND(p,s)
! If there is production during t then there is a setup in t
forall(p in PRODUCTS, t in TIMES)
ProdSetup(p,t):= produce(p,t) <= D(p,t,getlast(TIMES)) * setup(p,t)
! Capacity limits
forall(t in TIMES) Capacity(t):= sum(p in PRODUCTS) produce(p,t) <= CAP(t)
! Variables setup are 0/1
forall(p in PRODUCTS, t in TIMES) setup(p,t) is_binary
! Uncomment to get detailed MIP output
setparam("XPRS_VERBOSE", true)
! All cost data are integer, we therefore only need to search for integer
! solutions
setparam("XPRS_MIPADDCUTOFF", -0.999)
! Set Mosel comparison tolerance to a sufficiently small value
setparam("ZEROTOL", EPS/100)
writeln("**************ALG=",ALG,"***************")
starttime:=gettime
SEVERALROUNDS:=false; TOPONLY:=false
case ALG of
1: setparam("XPRS_CUTSTRATEGY", 0) ! No cuts
2: setparam("XPRS_PRESOLVE", 0) ! No presolve
3: tree_cut_gen ! User branch-and-cut + automatic cuts
4: do ! User branch-and-cut (several rounds),
tree_cut_gen ! no automatic cuts
setparam("XPRS_CUTSTRATEGY", 0)
SEVERALROUNDS:=true
end-do
5: do ! User cut-and-branch (several rounds)
tree_cut_gen ! + automatic cuts
SEVERALROUNDS:=true
TOPONLY:=true
end-do
6: do ! User branch-and-cut (several rounds)
tree_cut_gen ! + automatic cuts
SEVERALROUNDS:=true
end-do
7: lp_cut_gen ! User cut-and-branch using loop
! around LP solving and save/load basis
end-case
minimize(MinCost) ! Solve the problem
writeln("Time: ", gettime-starttime, "sec, Nodes: ", getparam("XPRS_NODES"),
", Solution: ", getobjval)
writeln(" "*15, "| Prod. Quantity (Demand)")
write("Period Setup | ")
forall(p in PRODUCTS) write(strfmt(p,-8))
write("\n", "-"*50)
forall(t in TIMES) do
write("\n ", strfmt(t,2), strfmt(getsol(sum(p in PRODUCTS) setup(p,t)),8), " |")
forall(p in PRODUCTS)
write(strfmt(produce(p,t).sol,4,0), " (", DEMAND(p,t), ")")
end-do
writeln
!*************************************************************************
! Cut generation loop at the top node:
! solve the LP and save the basis
! get the solution values
! identify and set up violated constraints
! load the modified problem and load the saved basis
!*************************************************************************
procedure lp_cut_gen
declarations
ncut,npass,npcut:integer ! Counters for cuts and passes
objval,ds: real
cut: array(range) of linctr
bas: basis
end-declarations
setparam("XPRS_CUTSTRATEGY", 0) ! Disable automatic cuts (optional)
setparam("XPRS_PRESOLVE", 0) ! Switch presolve off (required)
ncut:=0
npass:=0
repeat
npass+=1
npcut:= 0
minimize(XPRS_LIN+XPRS_PRI, MinCost) ! Solve the LP using primal simplex
savebasis(bas) ! Save the current basis
objval:= getobjval ! Get the objective value
forall(p in PRODUCTS,t in TIMES) do ! Get the solution values
solprod(p,t):=produce(p,t).sol
solsetup(p,t):=setup(p,t).sol
end-do
! Search for violated constraints:
forall(p in PRODUCTS,l in TIMES) do
ds:=0
forall(t in 1..l)
if(solprod(p,t) < D(p,t,l)*solsetup(p,t) + EPS) then ds += solprod(p,t)
else ds += D(p,t,l)*solsetup(p,t)
end-if
! Add the violated inequality: the minimum of the actual production
! prod(p,t) and the maximum potential production D(p,t,l)*setup(p,t)
! in periods 1 to l must at least equal the total demand in periods
! 1 to l.
if(ds < D(p,1,l) - EPS) then
ncut+=1
npcut+=1
cut(ncut):= sum(t in 1..l)
if(solprod(p,t)<(D(p,t,l)*solsetup(p,t))+EPS, produce(p,t),
D(p,t,l)*setup(p,t)) >= D(p,1,l)
end-if
end-do
writeln("Pass ", npass, " (", gettime-starttime, " sec), objective value ",
objval, ", cuts added: ", npcut, " (total ", ncut,")")
if(npcut=0) then
writeln("Optimal integer solution found:")
else
loadprob(MinCost) ! Reload the problem
loadbasis(bas) ! Load the saved basis
end-if
until (npcut<=0 or npass>15)
end-procedure
!*************************************************************************
! Cut generation loop:
! get the solution values
! identify and set up violated constraints
! load the modified problem and load the saved basis
!*************************************************************************
public function cb_node:boolean
declarations
ncut:integer ! Counter for cuts
cut: array(range) of linctr ! Cuts
cutid: array(range) of integer ! Cut type identification
type: array(range) of integer ! Cut constraint type
objval,ds: real
end-declarations
depth:=getparam("XPRS_NODEDEPTH")
if((TOPONLY and depth<1) or (not TOPONLY and depth<=CUTDEPTH)) then
ncut:=0
! Get the solution values
forall(t in TIMES, p in PRODUCTS) do
solprod(p,t):=produce(p,t).sol
solsetup(p,t):=setup(p,t).sol
end-do
! Search for violated constraints
forall(p in PRODUCTS,l in TIMES) do
ds:=0
forall(t in 1..l)
if(solprod(p,t) < D(p,t,l)*solsetup(p,t) + EPS) then ds += solprod(p,t)
else ds += D(p,t,l)*solsetup(p,t)
end-if
! Generate the violated inequality
if(ds < D(p,1,l) - EPS) then
cut(ncut):= sum(t in 1..l)
if(solprod(p,t)<(D(p,t,l)*solsetup(p,t))+EPS, produce(p,t),
D(p,t,l)*setup(p,t)) - D(p,1,l)
cutid(ncut):= 1
type(ncut):= CT_GEQ
ncut+=1
end-if
end-do
returned:=false ! Call this function once per node
! Add cuts to the problem
if(ncut>0) then
addcuts(cutid, type, cut);
if(getparam("XPRS_VERBOSE")=true) then
writeln("Cuts added : ", ncut, " (depth ", depth, ", node ",
getparam("XPRS_NODES"), ", obj. ", getparam("XPRS_LPOBJVAL"), ")")
end-if
if SEVERALROUNDS then
returned:=true ! Repeat until no new cuts generated
end-if
end-if
end-if
end-function
! ****Optimizer settings for using the cut manager****
procedure tree_cut_gen
setparam("XPRS_PRESOLVE", 0) ! Switch presolve off
setparam("XPRS_EXTRAROWS", 5000) ! Reserve extra rows in matrix
setcallback(XPRS_CB_CM, "cb_node") ! Set the cut-manager callback function
end-procedure
end-model
Economic Lot-Sizing (ELS)
=========================
A well-known class of valid inequalities for ELS are the
(l,S)-inequalities. Let D(pq) denote the total demand in periods p
to q and y(t) be a binary variable indicating whether there is any
production in period t. For each period l and each subset of periods S
of 1 to l, the (l,S)-inequality is
sum (t=1:l | t in S) x(t) + sum (t=1:l | t not in S) D(tl) * y(t)
>= D(1l)
It says that actual production x(t) in periods included S plus maximum
potential production D(tl)*y(t) in the remaining periods (those not in
S) must at least equal total demand in periods 1 to l. Note that in
period t at most D(tl) production is required to meet demand up to
period l.
Based on the observation that
sum (t=1:l | t in S) x(t) + sum (t=1:l | t not in S) D(tl) * y(t)
>= sum (t=1:l) min(x(t), D(tl) * y(t))
>= D(1l)
it is easy to develop a separation algorithm and thus a cutting plane
algorithm based on these (l,S)-inequalities.
|
