Initializing help system before first use

# Column generation: solving different models in sequence

Topics covered in this section:

- Example problem: cutting stock
- Implementation
  - Main model
  - Knapsack model
- Alternative implementation with multiple problems
  - Main problem
  - Knapsack problem
- Results

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 (DEMANDi). 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,..., PATTERNSNWIDTH), 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 usej 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

usej

∑
j ∈ WIDTHS

PATTERNSij · usej ≥ DEMANDi
∀ j ∈ WIDTHS: usej ≤ ⌈DEMANDj / PATTERNSjj⌉ , usej ∈ ℕ

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

Ci · xj

∑
j ∈ WIDTHS

Aj · xj ≤ B
∀ j ∈ WIDTHS: xj 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

Parent Topic Mosel Parallel Whitepaper

© 2001-2026 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.