| 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
|
|
