Column generation

Topics covered in this section:

The technique of column generation is used for solving linear problems with a huge number of variables for which it is not possible to generate explicitly all columns of the problem matrix. Starting with a very restricted set of columns, after each solution of the problem a column generation algorithm adds one or several columns that improve the current solution. These columns must have a negative reduced cost (in a minimization problem) and are calculated based on the dual value of the current solution.

For solving large MIP problems, column generation typically has to be combined with a Branch-and-Bound search, leading to a so-called Branch-and-Price algorithm. The example problem described below is solved by solving a sequence of LPs without starting a tree search.

Example problem

A paper mill produces rolls of paper of a fixed width MAXWIDTH that are subsequently cut into smaller rolls according to the customer orders. The rolls can be cut into NWIDTHS different sizes. The orders are given as demands for each width i (DEMAND_i). The objective of the paper mill is to satisfy the demand with the smallest possible number of paper rolls in order to minimize the losses.

Model formulation

The objective of minimizing the total number of rolls can be expressed as choosing the best set of cutting patterns for the current set of demands. Since it may not be obvious how to calculate all possible cutting patterns by hand, we start off with a basic set of patterns (PATTERNS_1,..., PATTERNS_NWIDTH), that consists of cutting small rolls all of the same width as many times as possible (and at most the demanded quantity) out of the large roll. This type of problem is called a cutting stock problem.

If we define variables use_j to denote the number of times a cutting pattern j (j ∈ WIDTHS = {1,...,NWIDTH}) is used, then the objective becomes to minimize the sum of these variables, subject to the constraints that the demand for every size has to be met.

minimize

∑

_{j ∈ WIDTHS}

use_j

∑

_{j ∈ WIDTHS}

PATTERNS_ij · use_j ≥ DEMAND_i
∀ j ∈ WIDTHS: use_j ≤ ceil(DEMAND_j / PATTERNS_jj), use_j ∈ ℕ

Function ceil means rounding to the next larger integer value.

Implementation

The first part of the Mosel model paper.mos implementing this problem looks as follows:

model Papermill
 uses "mmxprs"

 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 procedure show_new_pat(dj:real, vx: array(range) of integer)

 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

  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]

 forall(j in WIDTHS)                    ! Make basic patterns
   PATTERNS(j,j) := minlist(floor(MAXWIDTH/WIDTH(j)),DEMAND(j))

 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

 forall(i in WIDTHS)                    ! Satisfy all demands
  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)),", ")

The paper mill can satisfy the demand with just the basic set of cutting patterns, but it is likely to incur significant losses through wasting more than necessary of every large roll and by cutting more small rolls than its customers have ordered. We therefore employ a column generation heuristic to find more suitable cutting patterns.

The following procedure column_gen defines a column generation loop that is executed at the top node (this heuristic was suggested by M. Savelsbergh for solving a similar cutting stock problem). The column generation loop performs the following steps:

solve the LP and save the basis
get the solution values
compute a more profitable cutting pattern based on the current solution
generate a new column (= cutting pattern): add a term to the objective function and to the corresponding demand constraints
load the modified problem and load the saved basis

To be able to increase the number of variables use_j in this function, these variables have been declared at the beginning of the program as a dynamic array without specifying any index range.

By setting Mosel's comparison tolerance to EPS, the test zbest = 0 checks whether zbest lies within EPS of 0 (see explanation in Section Cut generation).

Switching off presolve for the column generation problem generally helps to improve performance when iteratively resolving the problem after adding a new column and warm-starting it with the previous 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, npass) - 1.0

    write("Pass ", npass, ": ")
    if zbest = 0 then
      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.

      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

The preceding procedure column_gen calls the following auxiliary function knapsack to solve an integer knapsack problem of the form

maximize z =

∑

_{j ∈ WIDTHS}

C_i · x_j

∑

_{j ∈ WIDTHS}

A_j · x_j ≤ B
∀ j ∈ WIDTHS: x_j integer
∀ j ∈ WIDTHS: x_j ≤ D_j

The function knapsack solves a second optimization problem that is independent of the main, cutting stock problem since the two have no variables in common. We thus effectively work with two problems in a single Mosel model.

For efficiency reasons we have defined the knapsack variables and constraints globally. The integrality condition on the knapsack variables remains unchanged between several calls to this function, so we establish it when solving the first knapsack problem. On the other hand, the knapsack constraint and the objective function have different coefficients at every execution, so we need to replace them every time the function is called.

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

To complete the model, we add the following procedure show_new_pat to print every new pattern we find.

 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

Alternative implementation: Working with multiple problems

The implementation of the function knapsack in the previous section uses the sethidden functionality to blend out parts of the problem definition. The two parts of the problem (the main cutting stock problem and the problem solved in the knapsack routine) do not have any elements in common, that is, we really are solving two different problems within a single model.

An alternative implementation (available since Mosel 3.0) consists in formulating this model as two separate problems within the same model file.

The implementation as two separate problems in the model file papers.mos requires only few changes to the previous model formulation:

The declaration of a subproblem 'Knapsack' is added to the global declarations at the start of the model definition.

 declarations
  Knapsack: mpproblem                  ! Knapsack subproblem
 end-declarations

The implementation of function knapsack now works within the subproblem 'Knapsack' instead of hiding and unhiding subsets of the constraints. The scope of the subproblem is marked by the keywords with mpproblem [do ... end-do].

 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

  with Knapsack do

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

  end-do

 end-function

© 2001-2021 Fair Isaac Corporation. All rights reserved. This documentation is the property of Fair Isaac Corporation (“FICO”). Receipt or possession of this documentation does not convey rights to disclose, reproduce, make derivative works, use, or allow others to use it except solely for internal evaluation purposes to determine whether to purchase a license to the software described in this documentation, or as otherwise set forth in a written software license agreement between you and FICO (or a FICO affiliate). Use of this documentation and the software described in it must conform strictly to the foregoing permitted uses, and no other use is permitted.

Contents

Index

Glossary

Search Results

Column generation

Example problem

Model formulation

Implementation

Alternative implementation: Working with multiple problems