(!******************************************************* Mosel User Guide Examples ========================= file paperp.mos `````````````` Column generation algorithm for finding (near-)optimal solutions for a paper cutting example. Defining several (parallel) models. (c) 2008 Fair Isaac Corporation author: S. Heipcke, 2004, rev. May 2019 *******************************************************!) model "Papermill (multi-prob)" uses "mmxprs", "mmjobs", "mmsystem" forward procedure column_gen forward function knapsack(C:array(range) of real, A:array(range) of real, B:real, D:array(range) of integer, xbest:array(range) of integer): real forward procedure show_new_pat(dj:real, vx: array(range) of integer) (! The column generation algorithm requires the solution of a knapsack problem. This problem is defined in a separate model file. !) declarations NWIDTHS = 5 ! Number of different widths WIDTHS = 1..NWIDTHS ! Range of widths RP: range ! Range of cutting patterns MAXWIDTH = 94 ! Maximum roll width EPS = 1e-6 ! Zero tolerance WIDTH: array(WIDTHS) of real ! Possible widths DEMAND: array(WIDTHS) of integer ! Demand per width PATTERNS: array(WIDTHS, WIDTHS) of integer ! (Basic) cutting patterns use: dynamic array(RP) of mpvar ! Rolls per pattern soluse: dynamic array(RP) of real ! Solution values for variables `use' Dem: array(WIDTHS) of linctr ! Demand constraints MinRolls: linctr ! Objective function Knapsack: Model ! Reference to the knapsack model end-declarations WIDTH:: [ 17, 21, 22.5, 24, 29.5] DEMAND:: [150, 96, 48, 108, 227] ! Make basic patterns forall(j in WIDTHS) PATTERNS(j,j) := floor(MAXWIDTH/WIDTH(j)) (! Alternative starting patterns, particularly for small demand values: forall(j in WIDTHS) PATTERNS(j,j) := minlist(floor(MAXWIDTH/WIDTH(j)),DEMAND(j)) forall(j in WIDTHS) PATTERNS(j,j) := 1 !) forall(j in WIDTHS) do create(use(j)) ! Create NWIDTHS variables `use' use(j) is_integer ! Variables are integer and bounded use(j) <= integer(ceil(DEMAND(j)/PATTERNS(j,j))) end-do MinRolls:= sum(j in WIDTHS) use(j) ! Objective: minimize no. of rolls ! Satisfy all demands forall(i in WIDTHS) Dem(i):= sum(j in WIDTHS) PATTERNS(i,j) * use(j) >= DEMAND(i) res:= compile("knapsack.mos") ! Compile the knapsack model load(Knapsack, "knapsack.bim") ! Load the knapsack model column_gen ! Column generation at top node minimize(MinRolls) ! Compute the best integer solution ! for the current problem (including ! the new columns) writeln("Best integer solution: ", getobjval, " rolls") write(" Rolls per pattern: ") forall(i in RP) write(getsol(use(i)),", ") writeln fdelete("knapsack.bim") ! Cleaning up fdelete("shmem:indata"); fdelete("shmem:resdata") !************************************************************************ ! Column generation loop at the top node: ! solve the LP and save the basis ! get the solution values ! generate new column(s) (=cutting pattern) ! load the modified problem and load the saved basis !************************************************************************ procedure column_gen declarations dualdem: array(WIDTHS) of real xbest: array(WIDTHS) of integer dw, zbest, objval: real bas: basis end-declarations setparam("XPRS_PRESOLVE", 0) ! Switch presolve off setparam("zerotol", EPS) ! Set comparison tolerance of Mosel npatt:=NWIDTHS npass:=1 while(true) do minimize(XPRS_LIN, MinRolls) ! Solve the LP savebasis(bas) ! Save the current basis objval:= getobjval ! Get the objective value ! Get the solution values forall(j in 1..npatt) soluse(j):=getsol(use(j)) forall(i in WIDTHS) dualdem(i):=getdual(Dem(i)) ! Solve a knapsack problem zbest:= knapsack(dualdem, WIDTH, MAXWIDTH, DEMAND, xbest) - 1.0 write("Pass ", npass, ": ") if zbest = 0 then ! Same as zbest <= EPS writeln("no profitable column found.\n") break else show_new_pat(zbest, xbest) ! Print the new pattern npatt+=1 create(use(npatt)) ! Create a new var. for this pattern use(npatt) is_integer MinRolls+= use(npatt) ! Add new var. to the objective dw:=0 forall(i in WIDTHS) if xbest(i) > 0 then Dem(i)+= xbest(i)*use(npatt) ! Add new var. to demand constr.s dw:= maxlist(dw, ceil(DEMAND(i)/xbest(i) )) end-if use(npatt) <= dw ! Set upper bound on the new var. ! (useful when PRESOLVE is off) loadprob(MinRolls) ! Reload the problem loadbasis(bas) ! Load the saved basis end-if npass+=1 end-do writeln("Solution after column generation: ", objval, " rolls, ", getsize(RP), " patterns") write(" Rolls per pattern: ") forall(i in RP) write(soluse(i),", ") writeln setparam("XPRS_PRESOLVE", 1) ! Switch presolve on end-procedure !************************************************************************ ! Solve the integer knapsack problem ! z = max{cx : ax<=b, x<=d, x in Z^N} ! with b=MAXWIDTH, N=NWIDTHS, d=DEMAND ! ! All data is exchanged between the two models via shared memory. !************************************************************************ function knapsack(C:array(range) of real, A:array(range) of real, B:real, D:array(range) of integer, xbest:array(range) of integer):real INDATA := "bin:shmem:indata" RESDATA := "bin:shmem:resdata" initializations to INDATA A B C D end-initializations run(Knapsack, "NWIDTHS=" + NWIDTHS + ",INDATA=" + INDATA + ",RESDATA=" + RESDATA) ! Start knapsack (sub)model wait ! Wait until subproblem finishes dropnextevent ! Ignore termination message initializations from RESDATA xbest returned as "zbest" end-initializations end-function !************************************************************************ procedure show_new_pat(dj:real, vx: array(range) of integer) declarations dw: real end-declarations writeln("new pattern found with marginal cost ", dj) write(" Widths distribution: ") dw:=0 forall(i in WIDTHS) do write(WIDTH(i), ":", vx(i), " ") dw += WIDTH(i)*vx(i) end-do writeln("Total width: ", dw) end-procedure end-model