(!******************************************************* Mosel User Guide Examples ========================= file paperm.mos ``````````````` Column generation algorithm for finding (near-)optimal solutions for a paper cutting example. Data in- and output via a static module. Model can only be executed from the file paper.c, and not in standalone mode. Heuristic generation of patterns at start. (c) 2008 Fair Isaac Corporation author: S. Heipcke, 2002, rev. Dec. 2017 *******************************************************!) model Papermill uses "mmxprs" uses "cutstkdata" ! User module for in-memory data input parameters NWIDTHS = 5 ! Number of different widths WDATA='' ! Location of data in memory WSIZE=0 ! Size of the data block (nb of reals) DDATA='' ! Location of data in memory DSIZE=0 ! Size of the data block (nb of integers) end-parameters 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, pass: integer): real forward function pat_heur:integer (! The column generation algorithm requires the solution of a knapsack problem. We define the constraints of the knapsack problem globally because Mosel does not delete them at the end of a function (due to its principle of incremental model formulation). !) declarations 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, range) 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 KnapCtr, KnapObj: linctr ! Knapsack constraint+objective x: array(WIDTHS) of mpvar ! Knapsack variables end-declarations ! WIDTH:: [ 17, 21, 22.5, 24, 29,5] ! DEMAND:: [150, 96, 48, 108, 227] getcutstkdata(WIDTH, WDATA, WSIZE) getcutstkdata(DEMAND, DDATA, DSIZE) ! 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 setparam("zerotol", EPS) ! Set Mosel comparison tolerance bnd:=pat_heur ! Create heuristic solution + ! new non-dominated patterns 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) 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 !************************************************************************ ! 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 ubnd, zbest, objval: real bas: basis ! LP basis end-declarations setparam("XPRS_CUTSTRATEGY", 0) ! Disable automatic cuts setparam("XPRS_PRESOLVE", 0) ! Switch presolve off 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, npass) - 1 write("Pass ", npass, ": ") if zbest = 0 then writeln("no profitable column found.\n") break else printpat(zbest, xbest, WIDTH) ! 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 ubnd:=0 forall(i in WIDTHS) if xbest(i) > 0 then Dem(i)+= xbest(i)*use(npatt) ! Add new var. to demand constr.s ubnd:= maxlist(ubnd, ceil(DEMAND(i)/xbest(i) )) end-if use(npatt) <= ubnd ! Set upper bound on the new var. 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 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 ! ! We reset the knapsack constraints to 0 at the end of this function ! so that they do not unnecessarily increase the size of the main ! problem. The same is true in the other sense: hiding the demand ! constraints while solving the knapsack problem makes life easier ! for the optimizer, but is not essential for getting the correct ! solution. !************************************************************************ function knapsack(C:array(range) of real, A:array(range) of real, B:real, D:array(range) of integer, xbest:array(range) of integer, pass: integer):real ! Hide the demand constraints forall(j in WIDTHS) sethidden(Dem(j), true) ! Define the knapsack problem KnapCtr := sum(j in WIDTHS) A(j)*x(j) <= B KnapObj := sum(j in WIDTHS) C(j)*x(j) ! Integrality condition if(pass=1) then forall(j in WIDTHS) x(j) is_integer forall(j in WIDTHS) x(j) <= D(j) end-if maximize(KnapObj) returned:=getobjval forall(j in WIDTHS) xbest(j):=round(getsol(x(j))) ! Reset knapsack constraint and objective, unhide demand constraints KnapCtr := 0 KnapObj := 0 forall(j in WIDTHS) sethidden(Dem(j), false) end-function !********************************************************************** ! This heuristic assumes that widths are given in ascending order ! * Calculate the theoretical minimum number `minpat' of rolls ! * Assign all items to rolls, starting with the largest items and ! `minpat' rolls. If an item does not fit, start a new roll. ! * If the resulting number of rolls is equal to `minpat': stop, ! optimal solution found. ! * Otherwise, complete new patterns to non-dominated patterns by ! adding the largest-possible item to every new pattern ! * Create the variables `use' for the new patterns !********************************************************************** function pat_heur:integer declarations NI: integer end-declarations NI:= sum(j in WIDTHS) DEMAND(j) declarations CurCont: array(range) of real ALLITEM: array(1..NI) of real INDALLITEM: array(1..NI) of integer mincont:real end-declarations ! Calculate the theoretical minimum number of rolls largest:=WIDTH(NWIDTHS) smallest:=WIDTH(1) newpat:=ceil((sum(j in WIDTHS) WIDTH(j)*DEMAND(j))/MAXWIDTH) minpat:=newpat returned:=newpat writeln("Lower bound: ", minpat) ! Copy all widths into the array ALLITEM, according to their occurrences ct:=1 forall(j in WIDTHS) forall(k in 1..DEMAND(NWIDTHS-j+1)) do ALLITEM(ct):= WIDTH(NWIDTHS-j+1) INDALLITEM(ct):= NWIDTHS-j+1 ct+=1 end-do ! Assign the largest items, one to every roll forall(i in 1..newpat) do PATTERNS(INDALLITEM(i),NWIDTHS+i):=1 CurCont(i)+=ALLITEM(i) end-do ! Add the remaining items to the rolls ct:=newpat+1 while(ct<=NI) do mincont:=MAXWIDTH indi:=0 forall(i in 1..newpat) ! Find the least filled roll `indi' if(CurCont(i)<mincont) then mincont:=CurCont(i) indi:=i end-if>=ALLITEM(ct)) then ! Add the item to roll `indi' PATTERNS(INDALLITEM(ct),NWIDTHS+indi)+=1 CurCont(indi)+=ALLITEM(ct) else ! Start a new roll newpat+=1 PATTERNS(INDALLITEM(ct),NWIDTHS+newpat)+=1 CurCont(newpat)+=ALLITEM(ct) end-if ct+=1 end-do ! Stop if optimal solution found if(newpat=minpat) then writeln("Heuristic solution is optimal: ", newpat, " patterns.") forall(k in 1..newpat) do forall(j in WIDTHS) write( if(PATTERNS(j,NWIDTHS+k)>0, " "+WIDTH(j)+":"+PATTERNS(j,NWIDTHS+k), "")) writeln(" ",sum(j in WIDTHS) PATTERNS(j,NWIDTHS+k)*WIDTH(j)) end-do exit(3) end-if ! Complete to non-dominated patterns forall(k in 1..newpat) do totwidth:=MAXWIDTH-sum(j in WIDTHS) PATTERNS(j,NWIDTHS+k)*WIDTH(j) if(totwidth >=smallest) then ! If roll is not filled repeat if(totwidth>largest) then ! As long as the largest item fits: r:=integer(ceil(NWIDTHS*random)) ! add a randomly chosen item PATTERNS(r,NWIDTHS+k)+=1 totwidth-=WIDTH(r) else ! Add the largest possible item ct:=0 while(ct<NWIDTHS and WIDTH(ct+1)<=totwidth) ct+=1>0) then
       PATTERNS(ct,NWIDTHS+k)+=1
       totwidth-=WIDTH(ct)
      end-if
     end-if 
    until (totwidth<smallest) end-if end-do (! Create the variables `use' for the new patterns ubnd:=max(k in WIDTHS) DEMAND(k) forall(k in 1..newpat) do create(use(NWIDTHS+k)) ! Create additional variables `use' use(NWIDTHS+k) is_integer ! Variables are integer and bounded use(NWIDTHS+k) <= ubnd end-do !) writeln("Heuristic solution: ", newpat," patterns") end-function end-model