FICO Xpress Optimization Examples Repository
FICO Optimization Community FICO Xpress Optimization Home
Back to examples browserPrevious exampleNext example

Column generation for a cutting stock problem

  • creating new variables and adding them to constraints and objective
  • basis in- and output
  • sorting and single knapsack algorithms
  • repeat-until statement
  • use of functions
  • hiding / unhiding constraints
  • resetting (deleting) constraints
The first version of this model (cutstk.mos) solves the knapsack subproblem heuristically, the second version (paper.mos) solves it as a second optimization problem, hiding (blending out) the constraints of the master problem. The third version defines the subproblem as a separate problem with the same model, repeatedly creating new problems (papersn.mos) or re-using/redefining the same problem (papers.mos), that is:
  • defining multiple problems within a model
  • switching between problems
A fourth version (master model: paperp[r].mos, submodel: knapsack[r].mos) shows how to implement the algorithm with two separate models, illustrating the following features of Mosel:
  • working with multiple models
  • executing a sequence of models
  • passing data via shared memory

Further explanation of this example: paper.mos, papers.mos, papersn.mos: 'Mosel User Guide', Section 11.2 Column generation, and the Xpress Whitepaper 'Embedding Optimization Algorithms', Section 'Column generation for cutting-stock' (also discusses a generalization to bin-packing problems); for the multiple model versions paperp.mos, paperms.mos see the Xpress Whitepaper 'Multiple models and parallel solving with Mosel', Section 'Column generation: solving different models in sequence'.

Source Files


   Mosel User Guide Examples

   file paper.mos
   Column generation algorithm for finding (near-)optimal
   solutions for a paper cutting example.
   (c) 2008 Fair Isaac Corporation
       authors: Bob Daniel, S. Heipcke, J. Tebboth, 2001
                rev. May 2019

model Papermill
 uses "mmxprs"

 forward procedure column_gen
 forward function knapsack(C:array(range) of real, 
                           A:array(range) of real, 
                           D:array(range) of integer, 
                           xbest:array(range) of integer,
                           pass: 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.
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).
  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
  KnapCtr, KnapObj: linctr             ! Knapsack constraint+objective
  x: array(WIDTHS) of mpvar            ! Knapsack variables

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

 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)),", ")

!  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

   dualdem: array(WIDTHS) of real 
   xbest: array(WIDTHS) of integer
   dw, zbest, objval: real 
   bas: basis

  setparam("XPRS_verbose", true)

  setparam("XPRS_PRESOLVE", 0)         ! Switch presolve off
  setparam("zerotol", EPS)             ! Set comparison tolerance of Mosel

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

    write("Pass ", npass, ": ")   
    if zbest = 0 then                  ! Same as zbest <= EPS
      writeln("no profitable column found.\n")
      show_new_pat(zbest, xbest)       ! Print the new pattern
      create(use(npatt))               ! Create a new var. for this pattern
      use(npatt) is_integer

      MinRolls+= use(npatt)            ! Add new var. to the objective
      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) )) 
      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

  writeln("Solution after column generation: ", objval, " rolls, ",
          getsize(RP), " patterns")
  write("   Rolls per pattern: ")
  forall(i in RP) write(soluse(i),", ")

  setparam("XPRS_PRESOLVE", 1)         ! Switch presolve on

!  Solve the integer knapsack problem
!    z = max{cx : ax<=b, x<=d, x in Z^N}
! 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, 
                   D:array(range) of integer, 
                   xbest:array(range) of integer,
                   pass: integer):real

  setparam("XPRS_PRESOLVE", 1)         ! Switch presolve on

! 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 and bounds
  if pass=1 then
   forall(j in WIDTHS) x(j) is_integer
   forall(j in WIDTHS) x(j) <= D(j)    ! These bounds can be omitted

  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)

  setparam("XPRS_PRESOLVE", 0)         ! Switch presolve off


 procedure show_new_pat(dj:real, vx: array(range) of integer)
   dw: real
  writeln("new pattern found with marginal cost ", dj)
  write("   Widths distribution: ")
  forall(i in WIDTHS) do 
    write(WIDTH(i), ":", vx(i), "  ")
    dw += WIDTH(i)*vx(i)
  writeln("Total width: ", dw)


Back to examples browserPrevious exampleNext example