ELS - Solving several model instances in parallel
|
|
|
| Type: | Lot sizing |
| Rating: | 5 (difficult) |
| Description: | In its basic version (els.mos) this model has the following features:
A third implementation (master model: runels.mos or runelst.mos, submodel: elsp.mos; version with SVG graph drawing: runels_graph.mos with elspg.mos) parallelizes the execution of several model instances, showing the following features:
|
| File(s): | runels.mos, runels_graph.mos, runelst.mos, elsp.mos (submodel), elspg.mos (submodel), elsc.mos, elscb.mos, elscb_graph.mos, elsglobal.mos |
| Data file(s): | els.dat, els4.dat, els5.dat |
|
|
|
| runels.mos |
(!*******************************************************
Mosel Example Problems
======================
file runels.mos
```````````````
Run several instances of the model elsp.mos in
parallel and coordinate the solution update.
-- Using the 'bin' IO driver --
*** ATTENTION: This model will return an error if ***
*** no more than one Xpress licence is available. ***
*** This model cannot be run with a Community Licence
for the provided data instance ***
(c) 2008 Fair Isaac Corporation
author: S. Heipcke, Nov. 2004, rev. Jul. 2017
*******************************************************!)
model "Els main"
uses "mmjobs", "mmsystem"
parameters
DATAFILE = "els.dat"
T = 15
P = 4
end-parameters
declarations
RM = 0..5 ! Range of models
TIMES = 1..T ! Time periods
PRODUCTS = 1..P ! Set of products
solprod: array(PRODUCTS,TIMES) of real ! Sol. values for var.s produce
solsetup: array(PRODUCTS,TIMES) of real ! Sol. values for var.s setup
DEMAND: array(PRODUCTS,TIMES) of integer ! Demand per period
modELS: array(RM) of Model ! Models
NEWSOL = 2 ! Identifier for "sol. found" event
Msg: Event ! Messages sent by models
params: text ! Submodel runtime parameters
end-declarations
! Compile, load, and run models M1 and M2
M1:= 1; M2:=3
res:= compile("elsp.mos") ! Compile the submodel
load(modELS(M1), "elsp.bim") ! Load submodels
load(modELS(M2), "elsp.bim")
forall(m in RM) modELS(m).uid:= m
setmodpar(params, "DATAFILE", DATAFILE)
setmodpar(params, "T", T); setmodpar(params, "P", P)
forall(m in [M1,M2]) do
setworkdir(modELS(m), ".")
setmodpar(params, "ALG", m); ! Configure submodels
run(modELS(m), params) ! Start submodel runs
end-do
objval:= MAX_REAL
algsol:= -1; algopt:= -1
repeat
wait ! Wait for the next event
Msg:= getnextevent ! Get the event
if getclass(Msg)=NEWSOL then ! Get the event class
if getvalue(Msg) < objval then ! Value of the event (= obj. value)
algsol:= Msg.fromuid ! ID of model sending the event
objval:= getvalue(Msg)
writeln("Improved solution ", objval, " found by model ", algsol)
forall(m in RM | m <> algsol) send(modELS(m), NEWSOL, objval)
else
writeln("Solution ", getvalue(Msg), " found by model ", Msg.fromuid)
end-if
end-if
until getclass(Msg)=EVENT_END ! A model has finished
algopt:= Msg.fromuid ! Retrieve ID of terminated model
forall(m in RM) stop(modELS(m)) ! Stop all running models
if algsol=-1 then
writeln("No solution available")
exit(1)
else ! Retrieve the best solution from shared memory
initializations from "bin:shmem:sol"+algsol
solprod
solsetup
end-initializations
initializations from DATAFILE
DEMAND
end-initializations
! Solution printing
writeln("Best solution found by model ", algsol)
writeln("Optimality proven by model ", algopt)
writeln("Objective value: ", objval)
write("Period setup ")
forall(p in PRODUCTS) write(strfmt(p,-7))
forall(t in TIMES) do
write("\n ", strfmt(t,2), strfmt(sum(p in PRODUCTS) solsetup(p,t),8), " ")
forall(p in PRODUCTS) write(strfmt(solprod(p,t),3), " (",DEMAND(p,t),")")
end-do
writeln
end-if
! Cleaning up temporary files
fdelete("elsp.bim")
fdelete("shmem:sol"+M1); fdelete("shmem:sol"+M2)
end-model
|
| runels_graph.mos |
(!*******************************************************
Mosel Example Problems
======================
file runels_graph.mos
`````````````````````
Run several instances of the model elspg.mos in
parallel and coordinate the solution update.
-- Using the 'bin' IO driver --
-- Graphical solution output --
*** ATTENTION: This model will return an error if ***
*** no more than one Xpress licence is available. ***
*** This model cannot be run with a Community Licence
for the provided data instance ***
(c) 2008 Fair Isaac Corporation
author: S. Heipcke, Nov. 2004, rev. July 2017
*******************************************************!)
model "Els main"
uses "mmjobs", "mmsystem", "mmsvg"
parameters
DATAFILE = "els5.dat" ! "els5.dat" !"els.dat" !"els4.dat"
T = 45 !45 !15 !60
P = 4
MAXDUR = 30000 ! Timer delay in millisec
end-parameters
declarations
RM = 0..5 ! Range of models
TIMES = 1..T ! Time periods
PRODUCTS = 1..P ! Set of products
solprod: array(PRODUCTS,TIMES) of real ! Sol. values for var.s produce
solsetup: array(PRODUCTS,TIMES) of real ! Sol. values for var.s setup
DEMAND: array(PRODUCTS,TIMES) of integer ! Demand per period
modELS: array(RM) of Model ! Models
NEWSOL = 2 ! Identifier for "sol. found" event
INITVAL = 3 ! "Initial LP value" event class
NEWGAP = 4 ! Identifier for "gapnotify" event
Msg: Event ! Messages sent by models
params: text ! Submodel runtime parameters
initval: array(RM) of real
lastsol: array(RM) of real
lastxb,lastxs,lastys,lastyb: array(RM) of real
end-declarations
! Compile, load, and run selected submodel versions
AlgSel:= [2,3,4] ![0,1,2,3,4,5]
res:= compile("elspg.mos") ! Compile the submodel
forall(a in AlgSel) load(modELS(a), "elspg.bim") ! Load submodels
forall(m in RM) modELS(m).uid:= m
setmodpar(params, "DATAFILE", DATAFILE)
setmodpar(params, "T", T); setmodpar(params, "P", P)
! Start submodel runs
forall(a in AlgSel) do
setmodpar(params, "ALG", a); run(modELS(a), params)
end-do
svgaddgroup("Msg", "", SVG_GREY)
forall(a in AlgSel) svgaddgroup(string(a), "Model with ALG="+a)
xoffset:=-150.0
svgsetgraphlabels("Time (in sec)", "MIP gap")
svgsetreffreq(5) ! High update frequency
svgrefresh
objval:= MAX_REAL
algsol:= -1; algopt:= -1; iffirst:=true; indl:= 0.0
forall(m in AlgSel) do
lastxb(m):=0.0; lastxs(m):=0.0; lastys(m):=0.0; lastyb(m):=0.0
end-do
tid:=settimer(0,MAXDUR,false) ! Start the timer (false=single measure)
repeat
repeat
wait(1) ! Wait for the next event
! Closing the graphical display will interrupt the algorithm
if svgclosing then
writeln("Stopped by closing display window")
break 2
end-if
until not isqueueempty
Msg:= getnextevent ! Get the event
if getclass(Msg)=NEWSOL then ! Get the event class
mid:= Msg.fromuid ! ID of model sending the event
val:= getvalue(Msg)
grtime:=2*gettime
if lastys(mid)<>initval(mid) then
svgaddline(string(mid), lastxs(mid), lastys(mid), grtime, val)
end-if
svgaddpoint(string(mid), grtime, val)
lastxs(mid):=grtime; lastys(mid):=val
svgaddline(string(mid), lastxb(mid), lastyb(mid), grtime, lastyb(mid))
lastxb(mid):=grtime
if getvalue(Msg) < objval then ! Value of the event (= obj. value)
algsol:= mid
objval:= val
writeln("Improved solution ", objval, " found by model ", algsol)
svgaddtext(string(algsol), xoffset, indl, "Improved solution "+ objval)
svgrefresh
indl+=10
forall(m in RM | m <> algsol) send(modELS(m), NEWSOL, objval)
else
writeln("Solution ", val, " found by model ", mid)
svgaddtext(string(mid), xoffset, indl, "Solution "+ val)
svgrefresh
indl+=10
end-if
elif getclass(Msg)=INITVAL then ! Get the event class
mid:= Msg.fromuid ! ID of model sending the event
initval(mid):= getvalue(Msg) ! Value of the event (= bound value)
lastys(mid):=initval(mid); lastyb(mid):=initval(mid)
grtime:=2*gettime
lastxb(mid):=grtime; lastxs(mid):=grtime
if iffirst then
svgsetgraphviewbox(xoffset-10,initval(mid),600,1000)
indl:=initval(mid)*1.01;
iffirst:=false
svgaddtext(string(mid), xoffset, indl, "Initial LP solution "+ getvalue(Msg))
svgrefresh
indl+=10
end-if
writeln("Initial LP solution ", getvalue(Msg), " from model ", mid)
elif getclass(Msg)=NEWGAP then ! Get the event class
mid:= Msg.fromuid ! ID of model sending the event
val:=getvalue(Msg) ! Value of the event (= bound value)
writeln("New best bound ", getvalue(Msg), " from model ", mid)
grtime:=2*gettime
if lastys(mid)<>initval(mid) then
svgaddline(string(mid), lastxs(mid), lastys(mid), grtime, lastys(mid))
end-if
lastxs(mid):=grtime
svgaddline(string(mid), lastxb(mid), lastyb(mid), grtime, val)
lastxb(mid):=grtime; lastyb(mid):=val
svgrefresh
elif getclass(Msg)=EVENT_TIMER then ! Handle timer events
svgaddtext("Msg", xoffset, indl, "Time limit reached")
svgrefresh
writeln("Time limit reached, stopping all submodels.")
break ! Exit from 'repeat' loop
elif getclass(Msg)=EVENT_END then
algopt:= Msg.fromuid ! Retrieve ID of terminated model
end-if
until getclass(Msg)=EVENT_END ! A model has finished
forall(m in RM) stop(modELS(m)) ! Stop all running models
if algsol=-1 then
writeln("No solution available")
exit(1)
else ! Retrieve the best solution from shared memory
initializations from "bin:shmem:sol"+algsol
solprod
solsetup
end-initializations
initializations from DATAFILE
DEMAND
end-initializations
! Solution printing
writeln("Best solution found by model ", algsol)
if algopt<>-1 then writeln("Optimality proven by model ", algopt); end-if
writeln("Objective value: ", objval)
write("Period setup ")
forall(p in PRODUCTS) write(strfmt(p,-7))
forall(t in TIMES) do
write("\n ", strfmt(t,2), strfmt(sum(p in PRODUCTS) solsetup(p,t),8), " ")
forall(p in PRODUCTS) write(strfmt(solprod(p,t),3), " (",DEMAND(p,t),")")
end-do
writeln
end-if
! Cleaning up temporary files
fdelete("elsp.bim")
forall(a in AlgSel) fdelete("shmem:sol"+a)
! svgsave("runels.svg")
svgwaitclose
end-model
|
| runelst.mos |
(!*******************************************************
Mosel Example Problems
======================
file runelst.mos
````````````````
Run several instances of the model elsp.mos in
parallel and coordinate the solution update.
-- Using the 'bin' IO driver --
-- Implements a timer to interrupt submodels --
*** ATTENTION: This model will return an error if ***
*** no more than one Xpress licence is available. ***
*** This model cannot be run with a Community Licence
for the provided data instance ***
(c) 2017 Fair Isaac Corporation
author: S. Heipcke, June 2017
*******************************************************!)
model "Els main"
uses "mmjobs", "mmsystem"
parameters
DATAFILE = "els5.dat" !els !els4
T = 45 !15 !60
P = 4
MAXDUR = 1000 ! Timer delay in millisec
end-parameters
declarations
RM = 0..5 ! Range of models
TIMES = 1..T ! Time periods
PRODUCTS = 1..P ! Set of products
solprod: array(PRODUCTS,TIMES) of real ! Sol. values for var.s produce
solsetup: array(PRODUCTS,TIMES) of real ! Sol. values for var.s setup
DEMAND: array(PRODUCTS,TIMES) of integer ! Demand per period
modELS: array(RM) of Model ! Models
NEWSOL = 2 ! Identifier for "sol. found" event
Msg: Event ! Messages sent by models
params: text ! Submodel runtime parameters
ifsol: boolean ! Whether a solution has been found
tid: integer ! Timer identifier
end-declarations
! Compile, load, and run models M1 and M2
M1:= 1; M2:=3
res:= compile("elsp.mos") ! Compile the submodel
load(modELS(M1), "elsp.bim") ! Load submodels
load(modELS(M2), "elsp.bim")
forall(m in RM) modELS(m).uid:= m
setmodpar(params, "DATAFILE", DATAFILE)
setmodpar(params, "T", T); setmodpar(params, "P", P)
tid:=settimer(0,MAXDUR,true) ! Start the timer (true=heartbeat)
forall(m in [M1,M2]) do
setworkdir(modELS(m), ".")
setmodpar(params, "ALG", m); ! Configure submodels
run(modELS(m), params) ! Start submodel runs
end-do
ifsol:= false
objval:= MAX_REAL
algsol:= -1; algopt:= -1
repeat
wait ! Wait for the next event
Msg:= getnextevent ! Get the event
if getclass(Msg)=NEWSOL then ! Get the event class
if getvalue(Msg) < objval then ! Value of the event (= obj. value)
algsol:= Msg.fromuid ! ID of model sending the event
objval:= getvalue(Msg)
writeln("Improved solution ", objval, " found by model ", algsol)
forall(m in RM | m <> algsol) send(modELS(m), NEWSOL, objval)
ifsol:= true
else
writeln("Solution ", getvalue(Msg), " found by model ", Msg.fromuid)
end-if
elif getclass(Msg)=EVENT_TIMER then ! Handle timer events
if ifsol then
writeln("Time limit reached, stopping all submodels.")
canceltimer(tid) ! Stopping the timer
break ! Exit from 'repeat' loop
else ! Timer enters new iteration to continue
writeln("No solution found - relaunching timer.")
end-if
elif getclass(Msg)=EVENT_END then
algopt:= Msg.fromuid ! Retrieve ID of terminated model
end-if
until getclass(Msg)=EVENT_END ! A model has finished
forall(m in RM) stop(modELS(m)) ! Stop all running models
if algsol=-1 then
writeln("No solution available")
exit(1)
else ! Retrieve the best solution from shared memory
initializations from "bin:shmem:sol"+algsol
solprod
solsetup
end-initializations
initializations from DATAFILE
DEMAND
end-initializations
! Solution printing
writeln("Best solution found by model ", algsol)
writeln("Optimality proven by model ", algopt)
writeln("Objective value: ", objval)
write("Period setup ")
forall(p in PRODUCTS) write(strfmt(p,-7))
forall(t in TIMES) do
write("\n ", strfmt(t,2), strfmt(sum(p in PRODUCTS) solsetup(p,t),8), " ")
forall(p in PRODUCTS) write(strfmt(solprod(p,t),3), " (",DEMAND(p,t),")")
end-do
writeln
end-if
! Cleaning up temporary files
fdelete("elsp.bim")
fdelete("shmem:sol"+M1); fdelete("shmem:sol"+M2)
end-model
|
| elsp.mos |
(!*******************************************************
Mosel Example Problems
======================
file elsp.mos
`````````````
Economic lot sizing, ELS, problem
(Cut generation algorithm adding (l,S)-inequalities
in one or several rounds at the root node or in
tree nodes)
Parallel version.
*** Can be run standalone - intended to be run from runels.mos ***
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. There is a
set-up up cost SETUPCOST[t] associated with
production in period t. The unit production cost
in period t is PRODCOST[p,t]. There is no inventory
or stock-holding cost.
(c) 2008 Fair Isaac Corporation
author: S. Heipcke, Nov. 2004, rev. Sep. 2014
*******************************************************!)
model Els
uses "mmxprs","mmsystem","mmjobs"
parameters
ALG = 0 ! Default algorithm: no user cuts
CUTDEPTH = 3 ! Maximum tree depth for cut generation
DATAFILE = "els4.dat"
T = 60
P = 4
end-parameters
forward public function cb_node: boolean
forward procedure tree_cut_gen
forward public function cb_updatebnd: boolean
forward public procedure cb_intsol
declarations
NEWSOL = 2 ! "New solution" event class
EPS = 1e-6 ! Zero tolerance
TIMES = 1..T ! Time periods
PRODUCTS = 1..P ! 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 real ! 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
Msg: Event ! An event
end-declarations
initializations from DATAFILE
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) )
! Production in period t must not exceed the total demand for the
! remaining periods; 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,T) * setup(p,t)
! Production in periods 0 to t must satisfy the total demand
! during this period of time
forall(p in PRODUCTS,t in TIMES)
sum(s in 1..t) produce(p,s) >= sum (s in 1..t) DEMAND(p,s)
! Capacity limits
forall(t in TIMES) 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
setparam("zerotol", EPS/100) ! Set Mosel comparison tolerance
starttime:=gettime
setparam("XPRS_THREADS", 1) ! No parallel threads for optimization
! Uncomment to get detailed MIP output
! setparam("XPRS_VERBOSE", true)
setparam("XPRS_LPLOG", 0)
setparam("XPRS_MIPLOG", -1000)
writeln("**************ALG=",ALG,"***************")
SEVERALROUNDS:=false; TOPONLY:=false
case ALG of
1: do
setparam("XPRS_CUTSTRATEGY", 0) ! No cuts
setparam("XPRS_HEURSTRATEGY", 0) ! No heuristics
end-do
2: do
setparam("XPRS_CUTSTRATEGY", 0) ! No cuts
setparam("XPRS_HEURSTRATEGY", 0) ! No heuristics
setparam("XPRS_PRESOLVE", 0) ! No presolve
end-do
3: tree_cut_gen ! User branch-and-cut (single round)
4: do ! User branch-and-cut (several rounds)
tree_cut_gen
SEVERALROUNDS:=true
end-do
5: do ! User cut-and-branch (several rounds)
tree_cut_gen
SEVERALROUNDS:=true
TOPONLY:=true
end-do
end-case
! Parallel setup
setcallback(XPRS_CB_PRENODE, "cb_updatebnd")
setcallback(XPRS_CB_INTSOL, "cb_intsol")
! Solve the problem
minimize(MinCost)
!*************************************************************************
! Cut generation loop:
! get the solution values
! identify violated constraints and add them as cuts to the problem
! re-solve the modified problem
!*************************************************************************
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):=getsol(produce(p,t))
solsetup(p,t):=getsol(setup(p,t))
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
! produce(p,t) and the maximum potential production D(p,t,l)*setup(t)
! in periods 1 to l must at least equal the total demand in periods
! 1 to l.
! sum(t=1:l) min(produce(p,t), D(p,t,l)*setup(p,t)) >= D(p,1,l)
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);
writeln("Model ", ALG, ": Cuts added : ", ncut,
" (depth ", depth, ", node ", getparam("XPRS_NODES"),
", obj. ", getparam("XPRS_LPOBJVAL"), ")")
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_HEURSTRATEGY", 0) ! Switch heuristics off
setparam("XPRS_CUTSTRATEGY", 0) ! Switch automatic cuts off
setparam("XPRS_PRESOLVE", 0) ! Switch presolve off
setparam("XPRS_EXTRAROWS", 5000) ! Reserve extra rows in matrix
! setparam("XPRS_EXTRAELEMS", 10*5000) ! Reserve extra matrix elements
setcallback(XPRS_CB_CUTMGR, "cb_node") ! Define the cut manager callback
end-procedure
!*************************************************************************
! Setup for parallel solving:
! check whether cutoff update required at every node
! store and communicate any new solution found
!*************************************************************************
! Update cutoff value
public function cb_updatebnd: boolean
if not isqueueempty then
repeat
Msg:= getnextevent
until isqueueempty
newcutoff:= getvalue(Msg)
cutoff:= getparam("XPRS_MIPABSCUTOFF")
writeln("Model ", ALG, ": New cutoff: ", newcutoff,
" old: ", cutoff)
if newcutoff<cutoff then
setparam("XPRS_MIPABSCUTOFF", newcutoff)
end-if
if (newcutoff < getparam("XPRS_LPOBJVAL")) then
returned:= true ! Node becomes infeasible
end-if
end-if
end-function
! Store and communicate new solution
public procedure cb_intsol
objval:= getparam("XPRS_LPOBJVAL") ! Retrieve current objective value
cutoff:= getparam("XPRS_MIPABSCUTOFF")
writeln("Model ", ALG, ": ", objval, " cutoff: ", cutoff)
if(cutoff > objval) then
setparam("XPRS_MIPABSCUTOFF", objval)
end-if
! Get the solution values
forall(t in TIMES, p in PRODUCTS) do
solprod(p,t):=getsol(produce(p,t))
solsetup(p,t):=getsol(setup(p,t))
end-do
! Store the solution in shared memory
initializations to "bin:shmem:sol"+ALG
solprod
solsetup
end-initializations
! Send "solution found" signal
send(NEWSOL, objval)
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.
|
| elspg.mos |
(!*******************************************************
Mosel Example Problems
======================
file elspg.mos
``````````````
Economic lot sizing, ELS, problem
(Cut generation algorithm adding (l,S)-inequalities
in one or several rounds at the root node or in
tree nodes)
Parallel version with gap notification.
*** Can be run standalone - intended to be run from runels_graph.mos ***
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. There is a
set-up up cost SETUPCOST[t] associated with
production in period t. The unit production cost
in period t is PRODCOST[p,t]. There is no inventory
or stock-holding cost.
(c) 2008 Fair Isaac Corporation
author: S. Heipcke, Nov. 2004, rev. July 2017
*******************************************************!)
model Els
uses "mmxprs","mmsystem","mmjobs"
parameters
ALG = 0 ! Default algorithm: no user cuts
CUTDEPTH = 3 ! Maximum tree depth for cut generation
DATAFILE = "els4.dat"
T = 60
P = 4
end-parameters
forward public function cb_node: boolean
forward procedure tree_cut_gen
forward public function cb_updatebnd: boolean
forward public procedure cb_intsol
forward public procedure cb_gapnotify(rt,at,aot,abt:real)
declarations
NEWSOL = 2 ! "New solution" event class
INITVAL = 3 ! "Initial LP value" event class
NEWGAP = 4 ! "Gap notify" event class
EPS = 1e-6 ! Zero tolerance
TIMES = 1..T ! Time periods
PRODUCTS = 1..P ! 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 real ! 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
Msg: Event ! An event
end-declarations
initializations from DATAFILE
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) )
! Production in period t must not exceed the total demand for the
! remaining periods; 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,T) * setup(p,t)
! Production in periods 0 to t must satisfy the total demand
! during this period of time
forall(p in PRODUCTS,t in TIMES)
sum(s in 1..t) produce(p,s) >= sum (s in 1..t) DEMAND(p,s)
! Capacity limits
forall(t in TIMES) 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
setparam("zerotol", EPS/100) ! Set Mosel comparison tolerance
starttime:=gettime
setparam("XPRS_THREADS", 1) ! No parallel threads for optimization
! Uncomment to get detailed MIP output
! setparam("XPRS_VERBOSE", true)
setparam("XPRS_LPLOG", 0)
setparam("XPRS_MIPLOG", -1000)
writeln("**************ALG=",ALG,"***************")
SEVERALROUNDS:=false; TOPONLY:=false
case ALG of
1: do
setparam("XPRS_CUTSTRATEGY", 0) ! No cuts
setparam("XPRS_HEURSTRATEGY", 0) ! No heuristics
end-do
2: do
setparam("XPRS_CUTSTRATEGY", 0) ! No cuts
setparam("XPRS_HEURSTRATEGY", 0) ! No heuristics
setparam("XPRS_PRESOLVE", 0) ! No presolve
end-do
3: tree_cut_gen ! User branch-and-cut (single round)
4: do ! User branch-and-cut (several rounds)
tree_cut_gen
SEVERALROUNDS:=true
end-do
5: do ! User cut-and-branch (several rounds)
tree_cut_gen
SEVERALROUNDS:=true
TOPONLY:=true
end-do
end-case
! Parallel setup
setcallback(XPRS_CB_PRENODE, "cb_updatebnd")
setcallback(XPRS_CB_INTSOL, "cb_intsol")
mipgap:= 0.30
setparam("XPRS_MIPRELGAPNOTIFY", mipgap)
setcallback(XPRS_CB_GAPNOTIFY, "cb_gapnotify")
! Solve the problem
minimize(XPRS_LPSTOP,MinCost)
send(INITVAL, getobjval) ! Send "initial value" signal
minimize(XPRS_CONT,MinCost)
!*************************************************************************
! Cut generation loop:
! get the solution values
! identify violated constraints and add them as cuts to the problem
! re-solve the modified problem
!*************************************************************************
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):=getsol(produce(p,t))
solsetup(p,t):=getsol(setup(p,t))
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
! produce(p,t) and the maximum potential production D(p,t,l)*setup(t)
! in periods 1 to l must at least equal the total demand in periods
! 1 to l.
! sum(t=1:l) min(produce(p,t), D(p,t,l)*setup(p,t)) >= D(p,1,l)
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);
writeln("Model ", ALG, ": Cuts added : ", ncut,
" (depth ", depth, ", node ", getparam("XPRS_NODES"),
", obj. ", getparam("XPRS_LPOBJVAL"), ")")
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_HEURSTRATEGY", 0) ! Switch heuristics off
setparam("XPRS_CUTSTRATEGY", 0) ! Switch automatic cuts off
setparam("XPRS_PRESOLVE", 0) ! Switch presolve off
setparam("XPRS_EXTRAROWS", 5000) ! Reserve extra rows in matrix
! setparam("XPRS_EXTRAELEMS", 10*5000) ! Reserve extra matrix elements
setcallback(XPRS_CB_CUTMGR, "cb_node") ! Define the cut manager callback
end-procedure
!*************************************************************************
! Setup for parallel solving:
! check whether cutoff update required at every node
! store and communicate any new solution found
!*************************************************************************
! Update cutoff value
public function cb_updatebnd: boolean
if not isqueueempty then
repeat
Msg:= getnextevent
until isqueueempty
newcutoff:= getvalue(Msg)
cutoff:= getparam("XPRS_MIPABSCUTOFF")
writeln("Model ", ALG, ": New cutoff: ", newcutoff,
" old: ", cutoff)
if newcutoff<cutoff then
setparam("XPRS_MIPABSCUTOFF", newcutoff)
end-if
if (newcutoff < getparam("XPRS_LPOBJVAL")) then
returned:= true ! Node becomes infeasible
end-if
end-if
end-function
! Store and communicate new solution
public procedure cb_intsol
objval:= getparam("XPRS_LPOBJVAL") ! Retrieve current objective value
cutoff:= getparam("XPRS_MIPABSCUTOFF")
writeln("Model ", ALG, ": ", objval, " cutoff: ", cutoff)
if(cutoff > objval) then
setparam("XPRS_MIPABSCUTOFF", objval)
end-if
! Get the solution values
forall(t in TIMES, p in PRODUCTS) do
solprod(p,t):=getsol(produce(p,t))
solsetup(p,t):=getsol(setup(p,t))
end-do
! Store the solution in shared memory
initializations to "bin:shmem:sol"+ALG
solprod
solsetup
end-initializations
! Send "solution found" signal
send(NEWSOL, objval)
! Send "gapnotify" signal
send(NEWGAP, getparam("XPRS_BESTBOUND"))
end-procedure
! Notify about gap changes
! With the setting XPRS_MIPRELGAPNOTIFY=0.20 this routine will be called first
! when gap reaches 20%. We then reset the target, so that it gets called
! once more at a 1% smaller gap
public procedure cb_gapnotify(rt,at,aot,abt:real)
writeln("Model ", ALG, ": Reached ", 100*mipgap, "% gap.")
mipobj:= getparam("XPRS_MIPOBJVAL")
bbound:= getparam("XPRS_BESTBOUND")
relgap:= abs( (mipobj-bbound)/mipobj )
if relgap<=mipgap then
! Call "setgndata" to return new target value to the Optimizer
mipgap-=0.0025
setgndata(XPRS_GN_RELTARGET, mipgap)
! Send "gapnotify" signal
send(NEWGAP, bbound)
end-if
if relgap<=0.005 then
setgndata(XPRS_GN_RELTARGET, -1) ! Don't call gapnotify callback any more
end-if
end-procedure
end-model
|
| elsc.mos |
(!*******************************************************
Mosel Example Problems
======================
file elsc.mos
`````````````
Economic lot sizing, ELS, problem
(Cut generation algorithm adding (l,S)-inequalities
in one or several rounds at the root node or in
tree nodes)
Standard version (configurable cutting plane algorithm).
Economic lot sizing (ELS) considers production planning
over a given planning horizon. In every period, there is
a given demand for every product that must be satisfied by
the production in this period and by inventory carried over
from previous periods.
A set-up cost is associated with production in a period,
and the total production capacity per period is limited.
The unit production cost per product and time period is
given. There is no inventory or stock-holding cost.
*** This model cannot be run with a Community Licence
for the provided data instance ***
(c) 2008 Fair Isaac Corporation
author: S. Heipcke, 2001, rev. June 2010
*******************************************************!)
model "ELS"
uses "mmxprs","mmsystem"
parameters
ALG = 6 ! Algorithm choice [default settings: 0]
CUTDEPTH = 10 ! Maximum tree depth for cut generation
EPS = 1e-6 ! Zero tolerance
end-parameters
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 "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
end-case
minimize(MinCost) ! Solve the problem
writeln("Time: ", gettime-starttime, "sec, Nodes: ", getparam("XPRS_NODES"),
", Solution: ", getobjval)
write("Period setup ")
forall(p in PRODUCTS) write(strfmt(p,-7))
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(getsol(produce(p,t)), " (",DEMAND(p,t),") ")
end-do
writeln
!*************************************************************************
! 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):=getsol(produce(p,t))
solsetup(p,t):=getsol(setup(p,t))
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_CUTMGR, "cb_node") ! Set the cut-manager callback
end-procedure
end-model
|
| elscb.mos |
(!*******************************************************
Mosel Example Problems
======================
file elscb.mos
``````````````
Economic lot sizing, ELS, problem
Basic version with large data set and logging callbacks.
Economic lot sizing (ELS) considers production planning
over a given planning horizon. In every period, there is
a given demand for every product that must be satisfied by
the production in this period and by inventory carried over
from previous periods.
A set-up cost is associated with production in a period,
and the total production capacity per period is limited.
The unit production cost per product and time period is
given. There is no inventory or stock-holding cost.
- Using the GAPNOTIFY callback -
(c) 2012 Fair Isaac Corporation
author: S. Heipcke, Sep. 2012, rev. Sep. 2014
*******************************************************!)
model "ELS with logging callbacks"
uses "mmxprs","mmsystem"
parameters
DATAFILE = "els4.dat"
T = 60
P = 4
end-parameters
forward public function cb_node: boolean
forward public procedure cb_intsol
forward public procedure cb_gapnotify(rt,at,aot,abt:real)
declarations
TIMES = 1..T ! Range of time
PRODUCTS = 1..P ! 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,logtime, objval, mipgap: real
end-declarations
initializations from DATAFILE
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)
! Setting callbacks for logging
setcallback(XPRS_CB_INTSOL, "cb_intsol")
setcallback(XPRS_CB_PRENODE, "cb_node")
mipgap:= 0.10
setparam("XPRS_MIPRELGAPNOTIFY", mipgap)
setcallback(XPRS_CB_GAPNOTIFY, "cb_gapnotify")
starttime:=gettime
logtime:=starttime
minimize(MinCost) ! Solve the problem
writeln("Time: ", gettime-starttime, "sec, Nodes: ", getparam("XPRS_NODES"),
", Solution: ", getobjval)
write("Period setup ")
forall(p in PRODUCTS) write(strfmt(p,-7))
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(getsol(produce(p,t)), " (",DEMAND(p,t),") ")
end-do
writeln
!*****************************************************************
! Function called at every B&B node, return value 'true' marks node as infeasible
public function cb_node: boolean
timeNow:=gettime
if timeNow-logtime>=5 then
bbound:= getparam("XPRS_BESTBOUND")
actnodes:= getparam("XPRS_ACTIVENODES")
writeln(timeNow-starttime, "sec. Best bound:", bbound, " best sol.:",
if(getparam("XPRS_MIPSTATUS")=XPRS_MIP_SOLUTION,
text(objval), text(" - ")),
" active nodes: ", actnodes)
logtime:=timeNow
end-if
returned:= false
end-function
! Store and display new solution
public procedure cb_intsol
lastobjval:=objval
objval:= getparam("XPRS_LPOBJVAL") ! Retrieve current objective value
writeln(gettime-starttime, "sec. New solution: ", objval)
! If model runs for more than 60sec and new solution is just slightly
! better, then interrupt search
if gettime-starttime>60 and abs(lastobjval-objval)<=5 then
writeln("Stopping search")
stopoptimize(XPRS_STOP_USER)
end-if
end-procedure
! Notify about gap changes
! With the setting XPRS_MIPRELGAPNOTIFY=0.10 this routine will be called first
! when gap reaches 10%. We then reset the target, so that it gets called
! once more at a 2% smaller gap
public procedure cb_gapnotify(rt,at,aot,abt:real)
writeln(gettime-starttime, "sec. Reached ",
100*mipgap, "% gap.")
mipobj:= getparam("XPRS_MIPOBJVAL")
bbound:= getparam("XPRS_BESTBOUND")
relgap:= abs( (mipobj-bbound)/mipobj )
if relgap<=0.1 then
! Call "setgndata" to return new target value to the Optimizer
mipgap-=0.02
setgndata(XPRS_GN_RELTARGET, mipgap)
end-if
if relgap<=0.02 then
setgndata(XPRS_GN_RELTARGET, -1) ! Don't call gapnotify callback any more
end-if
end-procedure
end-model
|
| elscb_graph.mos |
(!*******************************************************
Mosel Example Problems
======================
file elscb_graph.mos
````````````````````
Economic lot sizing, ELS, problem
Basic version with large data set and logging callbacks.
Economic lot sizing (ELS) considers production planning
over a given planning horizon. In every period, there is
a given demand for every product that must be satisfied by
the production in this period and by inventory carried over
from previous periods.
A set-up cost is associated with production in a period,
and the total production capacity per period is limited.
The unit production cost per product and time period is
given. There is no inventory or stock-holding cost.
- Using the GAPNOTIFY callback -
- Drawing progress graph of MIP gap -
(c) 2012 Fair Isaac Corporation
author: S. Heipcke, Sep. 2012, rev. Sep. 2017
*******************************************************!)
model "ELS with logging callbacks"
uses "mmxprs","mmsystem","mmsvg"
parameters
DATAFILE = "els4.dat"
T = 60
P = 4
end-parameters
forward public function cb_node: boolean
forward public procedure cb_intsol
forward public procedure cb_gapnotify(rt,at,aot,abt:real)
declarations
TIMES = 1..T ! Range of time
PRODUCTS = 1..P ! 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,logtime, objval, mipgap: real
end-declarations
initializations from DATAFILE
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)
! Setting callbacks for logging
setcallback(XPRS_CB_INTSOL, "cb_intsol")
setcallback(XPRS_CB_PRENODE, "cb_node")
mipgap:= 0.10
setparam("XPRS_MIPRELGAPNOTIFY", mipgap)
setcallback(XPRS_CB_GAPNOTIFY, "cb_gapnotify")
! Configuration of graphical output
svgaddgroup("msg", "Model output")
svgsetstyle(SVG_FONTSIZE, "6pt")
svgaddgroup("bbound", "Best bound")
svgaddgroup("sol", "Solutions")
xoffset:=-150.0
lastxb:=0.0; lastxs:=0.0;
svgsetgraphlabels("Time (in sec)", "MIP gap")
svgsetreffreq(5) ! High update frequency
svgrefresh
starttime:=gettime
logtime:=starttime
! Solve the problem
minimize(XPRS_LPSTOP,MinCost)
initval:=getobjval
indl:=initval*1.01; lastys:=initval; lastyb:=initval;
svgsetgraphviewbox(xoffset-10,initval,300,400)
minimize(XPRS_CONT,MinCost)
case getparam("XPRS_MIPSTATUS") of
XPRS_MIP_SOLUTION:
svgaddtext("msg", xoffset, indl, "Search incomplete")
XPRS_MIP_OPTIMAL:
svgaddtext("msg", xoffset, indl, "Optimality proven")
end-case
svgrefresh
writeln("Time: ", gettime-starttime, "sec, Nodes: ", getparam("XPRS_NODES"),
", Solution: ", getobjval)
write("Period setup ")
forall(p in PRODUCTS) write(strfmt(p,-7))
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(getsol(produce(p,t)), " (",DEMAND(p,t),") ")
end-do
writeln
svgwaitclose("Close browser window to terminate model execution.", 1)
!*****************************************************************
! Function called at every B&B node, return value 'true' marks node as infeasible
public function cb_node: boolean
if svgclosing then
writeln("Stopped by closing display window")
stopoptimize(XPRS_STOP_USER)
end-if
timeNow:=gettime
if timeNow-logtime>=5 then
bbound:= getparam("XPRS_BESTBOUND")
grtime:=2*gettime
if lastxs<>0.0 then
svgaddline("sol", lastxs, lastys,grtime, objval)
end-if
lastxs:=grtime; lastys:=objval
svgaddline("bbound", lastxb, lastyb,grtime, bbound)
lastxb:=grtime; lastyb:=bbound
svgrefresh
actnodes:= getparam("XPRS_ACTIVENODES")
writeln(timeNow-starttime, "sec. Best bound:", bbound, " best sol.:",
if(getparam("XPRS_MIPSTATUS")=XPRS_MIP_SOLUTION,
text(objval), text(" - ")),
" active nodes: ", actnodes)
logtime:=timeNow
end-if
returned:= false
end-function
! Store and display new solution
public procedure cb_intsol
lastobjval:=objval
objval:= getparam("XPRS_LPOBJVAL") ! Retrieve current objective value
writeln(gettime-starttime, "sec. New solution: ", objval)
bbound:= getparam("XPRS_BESTBOUND")
svgaddtext("msg", xoffset, indl, "New solution "+ objval)
indl+=10
grtime:=2*gettime
svgaddpoint("sol", grtime, objval)
svgsetstyle(svggetlastobj, SVG_STROKE, SVG_GREEN)
svgsetstyle(svggetlastobj, SVG_FILL, SVG_GREEN)
if lastxs<>0.0 then
svgaddline("sol", [lastxs, lastys, grtime, lastys, grtime, objval])
end-if
lastxs:=grtime; lastys:=objval
svgaddline("bbound", lastxb, lastyb,grtime, bbound)
lastxb:=grtime; lastyb:=bbound
svgrefresh
! If model runs for more than 30sec and new solution is just slightly
! better, then interrupt search
if gettime-starttime>30 and abs(lastobjval-objval)<=5 then
writeln("Stopping search")
svgaddtext("msg", xoffset, indl, "Stopping search at time limit")
indl+=10
svgrefresh
stopoptimize(XPRS_STOP_USER)
end-if
end-procedure
! Notify about gap changes
! With the setting XPRS_MIPRELGAPNOTIFY=0.10 this routine will be called first
! when gap reaches 10%. We then reset the target, so that it gets called
! once more at a 2% smaller gap
public procedure cb_gapnotify(rt,at,aot,abt:real)
writeln(gettime-starttime, "sec. Reached ",
100*mipgap, "% gap.")
svgaddtext("msg", xoffset, indl, "Reached "+100*mipgap+ "% gap")
svgrefresh
indl+=10
mipobj:= getparam("XPRS_MIPOBJVAL")
bbound:= getparam("XPRS_BESTBOUND")
relgap:= abs( (mipobj-bbound)/mipobj )
if relgap<=0.1 then
! Call "setgndata" to return new target value to the Optimizer
mipgap-=0.02
setgndata(XPRS_GN_RELTARGET, mipgap)
end-if
if relgap<=0.02 then
setgndata(XPRS_GN_RELTARGET, -1) ! Don't call gapnotify callback any more
end-if
end-procedure
end-model
|
| elsglobal.mos |
(!*******************************************************
Mosel Example Problems
======================
file elsglobal.mos
``````````````````
Economic lot sizing, ELS, problem
(Cut generation algorithm adding (l,S)-inequalities
in several rounds at tree nodes)
Model version with (optional) global tree cuts
Economic lot sizing (ELS) considers production planning
over a given planning horizon. In every period, there is
a given demand for every product that must be satisfied by
the production in this period and by inventory carried over
from previous periods.
A set-up cost is associated with production in a period,
and the total production capacity per period is limited.
The unit production cost per product and time period is
given. There is no inventory or stock-holding cost.
*** This model cannot be run with a Community Licence
for the provided data instance ***
(c) 2008 Fair Isaac Corporation
author: N. Verma, 2006, rev. Aug. 2018
*******************************************************!)
model "ELS (global cuts)"
uses "mmxprs","mmsystem"
parameters
CUTDEPTH = 10 ! Maximum tree depth for cut generation
EPS = 1e-6 ! Zero tolerance
GLOBAL = true ! Apply global cuts
end-parameters
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
num_cuts_added:integer
end-declarations
initializations from "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
! Without this setting behaviour of B&B is somewhat non-deterministic
setparam("XPRS_THREADS", 1)
! 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)
starttime:=gettime
tree_cut_gen ! User branch-and-cut (several rounds)
setparam("XPRS_CUTSTRATEGY", 0) ! No automatic cuts
minimize(MinCost) ! Solve the problem
writeln("Time: ", gettime-starttime, "sec, Nodes: ", getparam("XPRS_NODES"),
", Solution: ", getobjval)
write("Period setup ")
forall(p in PRODUCTS) write(strfmt(p,-7))
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(getsol(produce(p,t)), " (",DEMAND(p,t),") ")
end-do
writeln
!*************************************************************************
! 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
CutRange: range ! Counter for cuts
cut: array(CutRange) of linctr ! Cuts
cutid: array(CutRange) of integer ! Cut type identification
type: array(CutRange) of integer ! Cut constraint type
ndx_cuts: set of integer ! Index of cuts
active_cuts,violated_cuts: set of integer
objval,ds: real
cut_indices: set of integer
viol: dynamic array(range) of real
end-declarations
node:=getparam("XPRS_NODES")
depth:=getparam("XPRS_NODEDEPTH")
if (depth<=CUTDEPTH) then
ncut:=0
! Get the solution values
forall(t in TIMES, p in PRODUCTS) do
solprod(p,t):=getsol(produce(p,t))
solsetup(p,t):=getsol(setup(p,t))
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
if GLOBAL then
storecuts(0,cutid,type,cut,ndx_cuts)
writeln("Stored cuts at node ", node, " -> ", getsize(ndx_cuts))
getcplist(1,1,0, violated_cuts,viol)
writeln("Current violated cuts in the cutpool after solving node ",
node, " -> ", getsize(violated_cuts))
loadcuts(1,1,violated_cuts)
else
addcuts(cutid, type, cut);
num_cuts_added+=ncut
if(getparam("XPRS_VERBOSE")=true) then
writeln("Cuts added : ", ncut, " (depth ", depth, ", node ",
node, ", obj. ", getparam("XPRS_LPOBJVAL"), ") count=",
num_cuts_added)
end-if
end-if
returned:=true ! Repeat until no new cuts generated
else
getcplist(1,1,0, violated_cuts,viol)
if getsize(violated_cuts)>0 then
writeln("All violated cuts in the cutpool after solving node ", node,
" -> ", getsize(violated_cuts))
writeln("######################### VIOLATED #########################")
end-if
delcuts(true,0,-1,-MAX_REAL) ! Delete all cuts from problem at current node that are not required
! writeln("Cuts with non-basic slack are deleted at node ", node)
end-if
end-if
end-function
! ****Optimizer settings for using the cut manager****
procedure tree_cut_gen
setparam("XPRS_PRESOLVEOPS",2270) ! Non-destructive presolve only
! setparam("XPRS_PRESOLVE", 0) ! Switch presolve off
setparam("XPRS_EXTRAROWS", 5000) ! Reserve extra rows in matrix
setcallback(XPRS_CB_CUTMGR, "cb_node") ! Set the cut-manager callback
end-procedure
end-model
|
© 2001-2019 Fair Isaac Corporation. All rights reserved. This documentation is the property of Fair Isaac Corporation (“FICO”). Receipt or possession of this documentation does not convey rights to disclose, reproduce, make derivative works, use, or allow others to use it except solely for internal evaluation purposes to determine whether to purchase a license to the software described in this documentation, or as otherwise set forth in a written software license agreement between you and FICO (or a FICO affiliate). Use of this documentation and the software described in it must conform strictly to the foregoing permitted uses, and no other use is permitted.