| |||||||||||||||
Branch-and-Price for the Generalized Assignment Problem Description The model implements a branch-and-price algorithm that
solves a disaggregated formulation of
the Generalized Assignment Problem (GAP) where columns
represent feasible assignments of batches
to machines. Column generation is applied at every node of the
branch-and-bound tree. The branching algorithm is completely
implemented in Mosel, and the optimizer is used only to solve
the LP relaxation at each node. The model implementation shows the following features of Mosel:
Source Files By clicking on a file name, a preview is opened at the bottom of this page.
Data Files GAPbp3.mos (!************************************************************************** Mosel Example Problems ====================== file GAPbp3.mos ``````````````` TYPE: Branch-and-Price for Generalized Assignment Problem (GAP) DIFFICULTY: 5 FEATURES: mmjobs, Mosel 2.0 records. DESCRIPTION: A set of production batches is assigned to machines with the objective to maximize profit. The model implements a branch- and-price algorithm that solves a disaggregated formulation of the GAP where columns represent feasible assignments of batches to machines. Column generation is applied at every node of the branch-and-bound tree. Branching is based on fixing compact formulation variables that represent assignments of production batches to machines. The branch-and-bound logic is implemented in Mosel, and the optimizer is used only to solve the LP relaxation at each node. FURTHER INFO: Similar Problems: Mosel Parallel Whitepaper 3.1 Example problem: cutting stock 5.1 Example problem: multi-item, multi-period production planning References for Branch-and-Price: Savelsbergh, M.W.P.: "A Branch-and-Price Algorithm for the Generalized Assignment Problem" Oper. Res., v45, pp. 831- 841, 1997. Savelsbergh, M.W.P.: "Branch-and-price: Integer programming with column generation." In: Encyclopedia of Optimization, C. Floudas, pp. Pardalos, eds, 2001 *** ATTENTION: This model will return an error if no *** *** sufficient number of Xpress licences is available. *** The input data required by this model is generated by executing model "genGapDataD.mos" prior to running this model. (c) 2008 Fair Isaac Corporation author: Hernan Wurgaft, 2007, rev. S. Heipcke Sep. 2022 ***************************************************************************!) model GAPmas uses "mmxprs", "mmjobs", "mmsystem" parameters DATAFILE = "Data/Dtestx5_30_1.dat" TOL = 0.00001 DO_BRANCHING=true DO_HEURISTIC=false CUTOFF=0 BRANCH_STRATEGY=0 ! 0 is for variable branching (this is a placeholder for ! other branching strategies) BIG=1000000 SUBMOD="Data/GAPsubDP.mos" MESSLEV=4 ! 0: silent, 1-5: progressivly more detailed output end-parameters setparam("XPRS_PRESOLVE",0) declarations EVENT_GEN_INI_COLS=2 ! Event codes sent to submodels EVENT_GEN_COLS=3 EVENT_STOP_SUBMOD=5 EVENT_NEW_NODE=9 EVENT_SOLVED=6 ! Event codes sent by submodels EVENT_FAILED=7 EVENT_READY=8 EVENT_INFEAS=11 end-declarations ! Declare and read input data forward procedure read_data_size forward procedure read_prob_data declarations NM, NP : integer end-declarations read_data_size declarations MACH = 1..NM ! Set of machines PRODS = 1..NP ! Set of production batches CAP: array(MACH) of integer ! Machine capacities DUR: array(MACH,PRODS) of integer ! Production durations PROFIT: array(MACH,PRODS) of integer ! Production profit end-declarations read_prob_data ! Make data available for submodels to read initializations to "bin:shmem:data" DUR CAP PROFIT end-initializations !**** Branch-and-Price declarations**** declarations excessA: mpvar ! Violation of Assignment weight: array(MACH,range) of mpvar ! Weight for each column nPROP: array(MACH) of integer ! Number of columns from each machine Assign: array(PRODS) of linctr ! Assignment constraints Convex: array(MACH) of linctr ! Convexity constraints Mobjective: linctr ! Objective function of main problem USE: array(MACH,range,PRODS) of integer ! Stores columns as they are generated Price_convex: array(MACH) of real ! Dual prices for convexity constraints Price_assign: array(PRODS) of real ! Dual prices for assignment constraints branchvartype= record ! Data structure representing compact formulation variables: m: integer ! machine p: integer ! product end-record branchvar: branchvartype Node= record ! Data structure representing a node in the tree: ! Pointers refer to indexes in dynamic array nodes value: real ! optimal value of relaxation associated with node left,right,parent: integer ! pointers to parent and children in node tree feasible_sol: boolean integer_sol: boolean solbasis: basis ! basis associated with optimal value of relaxation prev_val, next_val: integer ! pointers to nodes in linked queue of ! active nodes (queue is maintained sorted by ! value) branchvar: branchvartype ! compact formulation variable to ! branch at (variable branching) end-record R: range nodes:dynamic array(R) of Node ! Array containing all nodes node: integer ! Current node being solved by column generation head: integer ! Head of queue of active nodes (0 if queue is empty) tail: integer ! Tail of queue of active nodes (0 if queue is empty) Select_node:integer ! Active node selected for branching Fixed= record ! Data structure representing fixed variables: var: branchvartype ! compact formulation variable that is fixed to_zero: boolean ! true if fixed to zero, false if fixed to one end-record fixed_vars: array(1..NP) of Fixed ! Array with all fixed variables fixedcolumns: set of integer Nfixed: integer ! Number of compact formulation variables fixed NWfix: integer ! Count of constraints to fix variables consfix: array(range) of linctr ! Constraints to fix variables sol_BEST: array(PRODS) of integer ! Incumbent solution in terms of compact ! formulation variables sol_GAP: array(MACH,PRODS) of real ! Solution in terms of compact ! formulation variables ubound: real ! Upper bound from Lagrangian relaxation ! calculated at every iteration Bestubound: real ! Best upper bound from Lagrangian ! relaxation Best_Value: real ! Incumbent value masbasis: basis ! Used to save basis at every column ! generation iteration Totiter: integer ! Count of column generation iterations countCols: integer ! Count of columns in disaggregated ! formulation integer_solution: boolean ! true if latest node solution is ! integer (set by translate_solution) end-declarations ! Column generation procedures forward procedure solve_node forward procedure process_sub_result forward procedure solve_main forward procedure process_main_result ! Branch-and-Price supporting procedures forward procedure translate_solution forward procedure branching_variable forward procedure variable_branching_reset forward procedure fixvar(branchvar:branchvartype,to_zero:boolean) ! (functions used to manage set fixedcolumns): forward function encode(m:integer,k:integer):integer forward function getm(e:integer):integer forward function getk(e:integer):integer ! Placeholder for heuristic forward procedure heuristic ! Reporting procedures forward procedure save_best_solution forward procedure write_sol_GAP ! Writes node solution in terms of ! compact representation variables forward procedure write_sol_BEST ! Writes incumbent solution in terms of ! compact representation variables forward procedure write_new_col(Machine:integer) forward procedure write_msg(msglevel:integer, msg:string) !************************************************* !************************************************* !*************** Get the Sub-models running !************************************************* !************************************************* declarations submod: array(MACH) of Model ! One submodel per machine modid: array(set of integer) of integer ! Model indices end-declarations res:= compile("g",SUBMOD,"tmp:gapsub.bim") ! Compile the submodel file forall(m in MACH) do ! Load & run submodels load(submod(m), "tmp:gapsub.bim") modid(getid(submod(m))):= m run(submod(m), "Machine=" + m + ",NM=" + NM + ",NP=" + NP ) !!All submodels are started in sequence !!They will wait until they receive signal to continue wait ! Wait for submodels to be ready ev:=getnextevent if ev.class=EVENT_END then writeln("*** Impossible to start all necessary models - aborting ***") exit(1) end-if end-do !************************************************* !************************************************* !*************** Solve Root Node !************************************************* !************************************************* declarations Machine: integer sol_use: array(PRODS) of integer sol_Profit: real end-declarations StartTime:= gettime ! Generate initial columns !************************* forall(m in MACH) send(submod(m),EVENT_GEN_INI_COLS,0) ! All submodels are sent signal to run in parallel forall(m in MACH) do ! This piece uses the memory pipe to wait and then processes ! the solution of each sub-model to generate initial columns initializations from "mempipe:noindex,sol" Machine sol_use sol_Profit end-initializations ! Add the new column to the main problem nPROP(Machine)+=1 countCols+=1 create(weight(Machine,nPROP(Machine))) weight(Machine,nPROP(Machine)) is_binary forall(p in PRODS) USE(Machine,nPROP(Machine),p):=sol_use(p) write_new_col(Machine) end-do ! Catches event that submodels completed generation of initial columns forall(m in MACH) do wait; dropnextevent; end-do ! Define and solve main problem with initial columns !*************************************************** forall(m in MACH) Convex(m):= sum (k in 1..nPROP(m)) weight(m,k) <= 1 ! Convexity constraints forall(p in PRODS) Assign(p):= sum(m in MACH, k in 1..nPROP(m)) USE(m,k,p) * weight(m,k)+ excessA = 1 ! Assignment constraints Mobjective:= sum(m in MACH, k in 1..nPROP(m)) (sum (p in PRODS) USE(m,k,p)* PROFIT(m,p))* weight(m,k)- 100 * excessA ! Objective function write_msg(2, "Solving initial main problem: " + (gettime-StartTime) + "sec") maximize(XPRS_LIN, Mobjective) savebasis(masbasis) write_msg(1, "Iter: 0 sol: " + getobjval) process_main_result Best_Value:=CUTOFF solve_node !************************************************* !************************************************* !************************************************* !************************************************* !************************************************* !************************************************* !*************** Begin Branch-and-Price !************************************************* !************************************************* !************************************************* !************************************************* !************************************************* !************************************************* if DO_BRANCHING then !************************************************* !************************************************* !*** Initialize root node in search tree active nodes !************************************************* !************************************************* write_msg(2, "Start branching: " + (gettime-StartTime) + "sec") node := 0 NWfix := 0 Nfixed := 0 translate_solution if integer_solution then save_best_solution else if BRANCH_STRATEGY=0: branching_variable if DO_HEURISTIC: heuristic end-if if getsol(excessA)>0 then nodes(node).feasible_sol:=false write_msg(1, "Infeasible") elif integer_solution then write_msg(1, "Integer Solution") nodes(node).integer_sol:=true else nodes(node).feasible_sol:=true nodes(node).integer_sol:=false nodes(node).value:=getobjval nodes(node).parent:=-1 nodes(node).prev_val:=-1 nodes(node).next_val:=-1 head:=0 tail:=0 Select_node:=0 nodes(node).left:= -1 nodes(node).right:=-1 end-if if nodes(0).feasible_sol and not nodes(0).integer_sol then !***************************************************** !*************** Begin looping over search tree nodes !***************************************************** search_done:=false repeat ! until search_done node+=1 nodes(node).parent:=Select_node nodes(node).left:= -1 nodes(node).right:=-1 !***************************************************** !** Reset the main problem for the current node. !***************************************************** (! Main problem columns that are fixed to zero for current node are first put in set fixedcolumns. Constraints consfix are used to fix the columns to zero. Also, information about fixed variables is passed to knapsack submodels. Procedures used: variable_branching_reset fixvar(branchvar:branchvartype,to_zero:boolean) Functions used (to manage set fixedcolumns): function encode(m:integer,k:integer):integer function getm(e:integer):integer function getk(e:integer):integer !) forall(c in 1..NWfix) consfix(c):=0 NWfix:=0 Nfixed:=0 fixedcolumns:={} if BRANCH_STRATEGY=0: variable_branching_reset initializations to "bin:shmem:fixed" Nfixed fixed_vars end-initializations forall(col in fixedcolumns) do NWfix+=1 consfix(NWfix):=weight(getm(col),getk(col))=0 end-do forall(m in MACH) send(submod(m), EVENT_NEW_NODE, 0) !***************************************************** !** Optimize main problem !***************************************************** loadprob(Mobjective) ! Reload the problem loadbasis(nodes(Select_node).solbasis) ! Load saved basis maximize(XPRS_LIN,Mobjective) savebasis(masbasis) process_main_result solve_node translate_solution if not integer_solution then !Determine how to branch off the current node active nodes if BRANCH_STRATEGY=0: branching_variable if DO_HEURISTIC: heuristic end-if !***************************************************** !** Update search tree based on solution of current node !** (Update status of all nodes using new information) !***************************************************** if getsol(excessA)>0 then nodes(node).feasible_sol:=false write_msg(3, "Infeasible") else nodes(node).feasible_sol:=true nodes(node).value:=Bestubound savebasis(nodes(node).solbasis) if integer_solution then write_msg(2, "Integer Solution Found") nodes(node).value:=getobjval nodes(node).integer_sol:=true if nodes(node).value>Best_Value then save_best_solution Best_Value:=nodes(node).value write_msg(2, "Improved integer value: " + Best_Value) if head=-1 or nodes(head).value-Best_Value<1.0 then search_done:=true end-if end-if else !Current node becomes active in the tree if nodes(node).value>=Best_Value+1 then position:=head !Tail repeat if position<>-1 and nodes(node).value < nodes(position).value then position:=nodes(position).next_val else !p1 will precede node and position will follow !p1 --- node --- position if position<>-1 then p1:=nodes(position).prev_val nodes(position).prev_val:=node else p1:=tail tail:=node end-if if p1<>-1 then nodes(p1).next_val:=node else head:=node end-if nodes(node).prev_val:=p1 nodes(node).next_val:=position break end-if until false else write_msg(3,"Cutoff node: " + node) end-if end-if end-if if MESSLEV>3 then writeln("************NODE REPORT*****************") sn:=node writeln("node ",sn," Value ",nodes(sn).value) sn:=Select_node writeln("parent ", sn, " Value ", nodes(sn).value, " P ", nodes(sn).parent," L ", nodes(sn).left," R ",nodes(sn).right) writeln("branch m ", nodes(sn).branchvar.m, " branch p ", nodes(sn).branchvar.p) writeln("time elapsed ", gettime, " total iter: ", Totiter, " total cols: ", countCols) writeln("******************************************") end-if !***************************************************** !** Select next node using Best Bound Rule !***************************************************** if head=-1 or nodes(head).value-Best_Value<1.0 then search_done:=true else Select_node:=head end-if until search_done end-if !if root node is not infeasible or integer end-if !if DO_BRANCHING !************************************************* !************************************************* forall(m in MACH) send(submod(m), EVENT_STOP_SUBMOD,0) ! Stop all models !Catch event indicating that all submodels stopped. forall(m in MACH) do wait; dropnextevent; end-do write_msg(1, "Total time: " + (gettime-StartTime) + "sec") write_msg(1, "Optimal objective value: " + Best_Value) write_sol_BEST !**************************************************** !**************PROCEDURES************************ !**************************************************** procedure solve_node (! Solve node relaxation applying column generation *************************************************** At each iteration, a knapsack submodel is solved for each machine. Knapsack solutions produce columns for the main problem and they are also used to calculate a Lagrangean upper bound. The loop is broken when no knapsack solution produces a column that improves the objective function (solved to optimality), when the main problem objective value is close enough to the upper bound (solved to almost optimality), or when the upper bound is smaller than the cutoff value (Best bound implies cutoff). Uses procedures: process_sub_result solve_main process_main_result !) Bestubound:=BIG Niter:=0 repeat Niter+=1 forall(m in MACH) ! Start solving all submodels send(submod(m), EVENT_GEN_COLS, 0) improve:=0 forall(m in MACH) do wait ev:=getnextevent if getclass(ev)<>EVENT_INFEAS then ubound+=getvalue(ev) else write_msg(3, "Infeasible knapsack") end-if !!(useful message) if getclass(ev)=EVENT_SOLVED then improve+=1 process_sub_result end-if end-do if ubound<Bestubound: Bestubound:=ubound if MESSLEV>2: writeln("Iter: ", Niter, " sol: ", getobjval, " ubound: ", ubound, " Bestubound ", Bestubound) if improve = 0 then write_msg(3, "solved to optimality") break end-if if Bestubound<Best_Value+1 then write_msg(3, "Bestbound implies cutoff") break end-if !If bound and value are close enough if integer(Bestubound)=integer(getobjval) then write_msg(3, "solved to almost optimality") break end-if solve_main process_main_result until false Totiter+=Niter end-procedure !**************************************************** procedure process_sub_result ! Add the new column to the main problem initializations from "mempipe:noindex,sol" Machine sol_use sol_Profit end-initializations nPROP(Machine)+=1 countCols+=1 create(weight(Machine,nPROP(Machine))) weight(Machine,nPROP(Machine)) is_binary forall(p in PRODS) USE(Machine,nPROP(Machine),p):=sol_use(p) Convex(Machine)+=weight(Machine,nPROP(Machine)) forall(p in PRODS) Assign(p)+= USE(Machine,nPROP(Machine),p) * weight(Machine,nPROP(Machine)) Mobjective+=(sum (p in PRODS)USE(Machine,nPROP(Machine),p)* PROFIT(Machine,p))* weight(Machine,nPROP(Machine)) write_new_col(Machine) end-procedure !**************************************************** procedure solve_main ! Optimize main problem loadprob(Mobjective) ! Reload the problem loadbasis(masbasis) ! Load the saved basis maximize(XPRS_LIN,Mobjective) end-procedure !**************************************************** procedure process_main_result (! Update dual pricing data for subproblems Add dual prices to upper bound calculation !) ubound:=0 savebasis(masbasis) forall(p in PRODS) do Price_assign(p):=getdual(Assign(p)) ubound+=Price_assign(p) end-do forall(m in MACH) do Price_convex(m):=getdual(Convex(m)) ubound+=Price_convex(m) end-do initializations to "bin:shmem:Price" Price_assign Price_convex end-initializations end-procedure !**************************************************** procedure translate_solution (! Translate solution from disaggregated representation variables (weight) to compact representation variables (sol_GAP). Check whether current solution is integer (integer_solution) !) integer_solution:=true if MESSLEV>1: writeln("Node ", node, " Objective: ",getobjval) forall(m in MACH) do ! Loop over machines forall(p in PRODS) sol_GAP(m,p):=0. forall(k in 1..nPROP(m)) do ! Loop over columns for machine if (getsol(weight(m,k))> TOL) then ! Finds active columns if abs(getsol(weight(m,k))-1.)>TOL then integer_solution:=false end-if forall(p in PRODS) do if (USE(m,k,p)=1) then sol_GAP(m,p)+=getsol(weight(m,k)) ! Translates variables end-if end-do end-if end-do end-do end-procedure !**************************************************** procedure branching_variable (! Determine the compact formulation variable to branch on for the current node !) branch_use :=1. branch_price :=0. forall(m in MACH) do ! Loop over machines forall(p in PRODS) do if (abs(sol_GAP(m,p)-0.5)<abs(branch_use-0.5)) then nodes(node).branchvar.m:=m nodes(node).branchvar.p:=p branch_use:=sol_GAP(m,p) branch_price:=PROFIT(m,p) elif (abs(sol_GAP(m,p)-0.5)=abs(branch_use-0.5)) and (PROFIT(m,p)>branch_price) then nodes(node).branchvar.m:=m nodes(node).branchvar.p:=p branch_use:=sol_GAP(m,p) branch_price:=PROFIT(m,p) end-if end-do end-do end-procedure !**************************************************** procedure variable_branching_reset (! Determine all compact formulation variables that need to be fixed at the current node. Right children are fixed to zero and left children are fixed to one. uses procedure: fixvar !) branchvar:= nodes(Select_node).branchvar !Fix variables if nodes(Select_node).left=-1 then ! Left node fixvar(branchvar,false) nodes(Select_node).left:=node else ! Right node fixvar(branchvar,true) nodes(Select_node).right:=node !!and take selected node out of queue head:=nodes(Select_node).next_val nodes(head).prev_val:=-1 end-if !Fix variable for all the parent nodes in the tree. visit:=Select_node while (nodes(visit).parent<>-1) do parent:=nodes(visit).parent branchvar:=nodes(parent).branchvar if visit=nodes(parent).left then fixvar(branchvar,false) else fixvar(branchvar,true) end-if visit:=parent end-do end-procedure !**************************************************** procedure fixvar(branchvar:branchvartype, to_zero: boolean) (! Build array fixed_vars and set fixedcolumns with compact and disaggregated formulation variables to be fixed. Array fixed_vars is to be later passed to knapsack subproblems. The set fixedcolumns is to be used to create consfix constraints. Implementation avoids duplicating consfix constraints by using function encode to place unique pairs m-k in set fixedcolumns. !) Nfixed+=1 fixed_vars(Nfixed).var:=branchvar case to_zero of true: do ! Fix to zero forall(k in 1..nPROP(branchvar.m)) do if USE(branchvar.m,k,branchvar.p)=1 then fixedcolumns+={encode(branchvar.m,k)} end-if end-do fixed_vars(Nfixed).to_zero:=true end-do false: do ! Fix to one fixed_vars(Nfixed).to_zero:=false forall(k in 1..nPROP(branchvar.m)) do if USE(branchvar.m,k,branchvar.p)<>1 then fixedcolumns+={encode(branchvar.m,k)} end-if end-do !Fix the rest to zero forall(m in MACH|m<>branchvar.m) do forall(k in 1..nPROP(m)) do if USE(m,k,branchvar.p)=1 then fixedcolumns+={encode(m,k)} end-if end-do end-do end-do end-case end-procedure !**************************************************** procedure heuristic writeln("Placeholder for heuristic") end-procedure !**************************************************** procedure write_sol_GAP forall(m in MACH) do write("Machine ", m, ": ") forall(p in PRODS) do if sol_GAP(m,p)>0 then write(p, "(",sol_GAP(m,p),")") end-if end-do writeln end-do end-procedure !**************************************************** procedure write_sol_BEST forall(m in MACH) do write("Machine ", m, ": ") forall(p in PRODS | sol_BEST(p)=m) write(p, " ") writeln end-do end-procedure !**************************************************** procedure save_best_solution forall(m in MACH,p in PRODS | sol_GAP(m,p)>0) sol_BEST(p):=m end-procedure !**************************************************** procedure write_new_col(Machine:integer) if MESSLEV>4 then write("Machine ", Machine, ": ") forall(p in PRODS | sol_use(p)=1) write(p, " ") writeln(" (total dur profit: ", sum(p in PRODS) DUR(Machine,p)*sol_use(p), " ",sol_Profit," ", CAP(Machine), ")") end-if end-procedure !**************************************************** procedure write_msg(msglevel:integer, msg:string) if msglevel<= MESSLEV: writeln(msg) end-procedure !**************************************************** function encode(m:integer,k:integer):integer returned:=m+k*1024 end-function function getm(e:integer):integer returned:=bittest(e,1023) end-function function getk(e:integer):integer returned:=e div 1024 end-function !**************************************************** procedure read_data_size initializations from DATAFILE NM NP end-initializations end-procedure procedure read_prob_data initializations from DATAFILE DUR CAP PROFIT end-initializations end-procedure end-model | |||||||||||||||
© Copyright 2024 Fair Isaac Corporation. |