| |||||||||
Cut generation for an economic lot-sizing (ELS) problem Description This model implements various forms of cut-and-branch and branch-and-cut algorithms. In its simplest form (looping over LP solving) it illustrates the following
features:
Another implementation (main model: runels.mos, submodel: elsp.mos) parallelizes the execution of several model instances, showing the following features:
Further explanation of this example: elscb.mos, elsglobal.mos, runels.mos: Xpress Whitepaper 'Multiple models and parallel solving with Mosel', Section 'Solving several model instances in parallel'.
Source Files By clicking on a file name, a preview is opened at the bottom of this page. Data Files 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 2023 *******************************************************!) 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_HEUREMPHASIS", 0) ! No heuristics end-do 2: do setparam("XPRS_CUTSTRATEGY", 0) ! No cuts setparam("XPRS_HEUREMPHASIS", 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 ds: real end-declarations returned:=false ! OPTNODE: This node is not infeasible depth:=getparam("XPRS_NODEDEPTH") cnt:=getparam("XPRS_CALLBACKCOUNT_OPTNODE") if ((TOPONLY and depth<1) or (not TOPONLY and depth<=CUTDEPTH)) and (SEVERALROUNDS or cnt<=1) 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 ! 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"), ")") end-if end-if end-function ! ****Optimizer settings for using the cut manager**** procedure tree_cut_gen setparam("XPRS_HEUREMPHASIS", 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 setcallback(XPRS_CB_OPTNODE, "cb_node") ! Set the optnode callback function 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 | |||||||||
© Copyright 2024 Fair Isaac Corporation. |