Column generation: solving different models in sequence

Topics covered in this section:

The cutting stock example we are working with in this section is taken from the `Mosel User Guide'. The reader is refered to this manual for further detail on the column generation algorithm and its implementation with Mosel.

Column generation algorithms are typically 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.

Our column generation algorithm for the cutting stock problem requires us to solve a knapsack problem based on the dual value of the current solution to determine a new column (= cutting pattern). The difference between the User Guide implementation and the one shown below consists in the handling of this knapsack (sub)problem. In the User Guide implementation Mosel's constraint hiding functionality is used to blend out subsets of constraints; in the version shown below the subproblem is implemented in a model on its own. Both versions implement exactly the same algorithm and their performance is comparable. On larger instances, however, the two-model version is likely to be slightly more efficient, since every model defines exactly the problem to be solved, without any selection of (un)hidden constraints.

In this example, the changes to the problems are such that they cause complete re-loading of the problems for every optimization run. A clearer advantage of the multi-model version would show up if there were only slight changes (bound updates) to the main (cutting stock) problem so that this problem did not have to be reloaded into the solver for every new run.

Example problem: cutting stock

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.

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.

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 ≤ ⌈DEMAND_j / PATTERNS_jj⌉ , use_j ∈ ℕ

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.

Our heuristic performs a column generation loop at the top node, before starting the MIP search. Every iteration of the column generation loop executes the following steps:

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

Step 3 of this loop requires us 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

This second optimization problem is independent of the main, cutting stock problem since the two have no variables in common.

Implementation

The implementation is divided into two parts: the main model (file paperp.mos) with the definition of the cutting stock problem and the column generation algorithm, and the knapsack model (file knapsack.mos) that is run as a submodel.

Main model

The main part of the cutting stock model looks as follows:

model "Papermill (multi-model)"
 uses "mmxprs", "mmjobs"

 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

 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]

 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

                                       ! 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
end-model

Before starting the column generation heuristic (the definition of procedure column_gen is left out here since it remains unchanged from the User Guide example) the knapsack model is compiled and loaded so that at every column generation loop we merely need to run it with new data. The knapsack model is run from the function knapsack that takes as its parameters the data for the knapsack problem and its solution values. The function saves all data to shared memory, then runs the knapsack model and retrieves the solution from shared memory. Its return value is the objective value (zbest) of the knapsack problem.

 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

To enforce a sequential execution of the two models (we need to retrieve the results from the knapsack problem before we may continue with the main problem) we must add a call to the procedure wait immediately after the run statement. Otherwise the execution of the parent model continues concurrently to the child model. On termination, the child model sends a `termination' event (an event of class EVENT_END). Since our algorithm does not require this event we simply remove it from the model's event queue with a call to dropnextevent.

Knapsack model

The implementation of the knapsack model is straightforward. All problem data is obtained from shared memory and after solving the problem its solution is saved into shared memory.

model "Knapsack"
 uses "mmxprs"

 parameters
  NWIDTHS=5                            ! Number of different widths
  INDATA = "knapsack.dat"              ! Input data file
  RESDATA = "knresult.dat"             ! Result data file
 end-parameters

 declarations
  WIDTHS = 1..NWIDTHS                  ! Range of widths
  A,C: array(WIDTHS) of real           ! Constraint + obj. coefficients
  B: real                              ! RHS value of knapsack constraint
  D: array(WIDTHS) of integer          ! Variables bounds (demand quantities)
  KnapCtr, KnapObj: linctr             ! Knapsack constraint+objective
  x: array(WIDTHS) of mpvar            ! Knapsack variables
  xbest: array(WIDTHS) of integer      ! Solution values
 end-declarations

 initializations from INDATA
  A  B  C  D
 end-initializations

! Define the knapsack problem
 KnapCtr:= sum(j in WIDTHS) A(j)*x(j) <= B
 KnapObj:= sum(j in WIDTHS) C(j)*x(j)

 forall(j in WIDTHS) x(j) is_integer
 forall(j in WIDTHS) x(j) <= D(j)

! Solve the problem and retrieve the solution
 maximize(KnapObj)
 z:=getobjval
 forall(j in WIDTHS) xbest(j):=round(getsol(x(j)))

 initializations to RESDATA
  xbest  z as "zbest"
 end-initializations

end-model

Alternative implementation with multiple problems

The solving of the main problem and the knapsack subproblem is sequential. An alternative implementation of the column generation algorithm therefore uses several problems (in the place of several models). Since the different problems are contained in a single model file, all Mosel code referring to data transfer between models is removed in this version; the mathematical part of the model remains the same. The complete Mosel program is somewhat shorter

Main problem

The main model body now simply defines the cutting stock problem, starts the column generation loop followed by the MIP optimization, and prints out the results, just like the model version documented in the Mosel User Guide.

The definition of the column generation algorithm in subroutine column_gen remains entirely unchanged and we therefore omit its listing here .

Knapsack problem

The complete formulation of the knapsack problem is now contained in the function knapsack, where the submodel is re-defined and re-solved at every call to the function. The line with mpproblem do indicates that we switch to a new problem, created locally at this point. After the do block we return automatically to the main (cutting stock) problem.

 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

  declarations
   x: array(WIDTHS) of mpvar            ! Knapsack variables
  end-declarations

  with mpproblem do                     ! Create a local subproblem

 ! Define the knapsack problem
   forall(j in WIDTHS) x(j) is_integer
   forall(j in WIDTHS) x(j) <= D(j)
   sum(j in WIDTHS) A(j)*x(j) <= B

   maximize(sum(j in WIDTHS) C(j)*x(j))
   returned:=getobjval
   forall(j in WIDTHS) xbest(j):=round(getsol(x(j)))

  end-do

 end-function

Slightly more efficient is the following version where the declaration of the knapsack problem has been moved into the main model body, that is, the problem is declared globally instead of being re-created and deleted at every call to the subroutine knapsack.

 declarations
  Knapsack: mpproblem                  ! Knapsack subproblem
  KnapCtr, KnapObj: linctr             ! Knapsack constraint+objective
  x: array(WIDTHS) of mpvar            ! Knapsack variables
 end-declarations


 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

Results

With the data in the model above, the column generation algorithm generates 6 new patterns, taking the value of the LP-relaxation of the cutting stock problem from originally 177.67 down to 160.95. The MIP finds a solution with 161 rolls using the following patterns:

	Widths
Pattern	17	21	22.5	24	29.5	Usage
3	0	0	4	0	0	1
5	0	0	0	0	3	15
6	0	1	0	3	0	32
8	2	0	0	0	2	75
10	0	2	1	0	1	32
11	0	0	2	2	0	6