Initializing help system before first use

# More about Integer Programming

Topics covered in this chapter:

- Cut generation
  - Example problem
  - Model formulation
  - Implementation
  - Cut-and-Branch
  - Comparison tolerance
  - Branch-and-Cut
- Column generation
  - Example problem
  - Model formulation
  - Implementation
  - Alternative implementation: Working with multiple problems

This chapter presents two applications to (Mixed) Integer Programming of the programming facilities in Mosel that have been introduced in the previous chapters.

## Cut generation

Cutting plane methods add constraints (cuts) to the problem that cut off parts of the convex hull of the integer solutions, thus drawing the solution of the LP relaxation closer to the integer feasible solutions and improving the bound provided by the solution of the relaxed problem.

The Xpress Optimizer provides automated cut generation (see the optimizer documentation for details). To show the effects of the cuts that are generated by our example we switch off the automated cut generation.

### Example problem

The problem we want to solve is the following: a large company is planning to outsource the cleaning of its offices at the least cost. The NSITES office sites of the company are grouped into areas (set AREAS={1,...,NAREAS}). Several professional cleaning companies (set CONTR={1,...,NCONTRACTORS}) have submitted bids for the different sites, a cost of 0 in the data meaning that a contractor is not bidding for a site.

To avoid being dependent on a single contractor, adjacent areas have to be allocated to different contractors. Every site s (s in SITES={1,...,NSITES}) is to be allocated to a single contractor, but there may be between LOWCONa and UPPCONa contractors per area a.

### Model formulation

For the mathematical formulation of the problem we introduce two sets of variables:

cleancs indicates whether contractor c is cleaning site s
allocca indicates whether contractor c is allocated any site in area a

The objective to minimize the total cost of all contracts is as follows (where PRICEsc is the price per site and contractor):

minimize

∑
c ∈ CONTR

∑
s ∈ SITES

PRICEsc · cleancs

We need the following three sets of constraints to formulate the problem:

1. Each site must be cleaned by exactly one contractor.

∀ s ∈ SITES:

∑
c ∈ CONTR

cleancs = 1

2. Adjacent areas must not be allocated to the same contractor.

∀ c ∈ CONTR, a,b ∈ AREAS, a > b and ADJACENTab = 1: allocca + alloccb ≤ 1

3. The lower and upper limits on the number of contractors per area must be respected.

∀ a ∈ AREAS:

∑
c ∈ CONTR

allocca ≥ LOWCONa
∀ a ∈ AREAS:

∑
c ∈ CONTR

allocca ≤ UPPCONa

To express the relation between the two sets of variables we need more constraints: a contractor c is allocated to an area a if and only if he is allocated a site s in this area, that is, yca is 1 if and only if some xcs (for a site s in area a) is 1. This equivalence is expressed by the following two sets of constraints, one for each sense of the implication (AREA_s is the area a site s belongs to and NUMSITEa the number of sites in area a):

∀ c ∈ CONTR, a ∈ AREAS: allocca ≤

∑
s ∈ SITES, AREAs=a

cleancs
∀ c ∈ CONTR, a ∈ AREAS: allocca ≥

1
NUMSITEa

·

∑
s ∈ SITES, AREAs=a

cleancs

### Implementation

The resulting Mosel program is the following model clean.mos. The variables cleancs are defined as a dynamic array and are only created if contractor c bids for site s (that is, PRICEsc>0 or, taking into account inaccuracies in the data, PRICEsc>0.01).

Another implementation detail that the reader may notice is the separate initialization of the array sizes: we are thus able to create all arrays with fixed sizes (with the exception of the previously mentioned array of variables that is explicitly declared dynamic). This allows Mosel to handle them in a more efficient way.

    model "Office cleaning"
     uses "mmxprs","mmsystem"

     declarations
       PARAM: array(1..3) of integer
     end-declarations

     initializations from 'clparam.dat'
       PARAM
     end-initializations

     declarations
       NSITES = PARAM(1)                    ! Number of sites
       NAREAS = PARAM(2)                    ! Number of areas (subsets of sites)
       NCONTRACTORS = PARAM(3)              ! Number of contractors
       AREAS = 1..NAREAS
       CONTR = 1..NCONTRACTORS
       SITES = 1..NSITES
       AREA: array(SITES) of integer        ! Area site is in
       NUMSITE: array(AREAS) of integer     ! Number of sites in an area
       LOWCON: array(AREAS) of integer      ! Lower limit on the number of
                                            ! contractors per area
       UPPCON: array(AREAS) of integer      ! Upper limit on the number of
                                            ! contractors per area
       ADJACENT: array(AREAS,AREAS) of integer    ! 1 if areas adjacent, 0 otherwise
       PRICE: array(SITES,CONTR) of real    ! Price per contractor per site

       clean: dynamic array(CONTR,SITES) of mpvar ! 1 iff contractor c cleans site s
       alloc: array(CONTR,AREAS) of mpvar   ! 1 iff contractor allocated to a site
                                            !  in area a
     end-declarations

     initializations from 'cldata.dat'
       [NUMSITE,LOWCON,UPPCON] as 'AREA'
       ADJACENT
       PRICE
     end-initializations

     ct:=1
     forall(a in AREAS) do
       forall(s in ct..ct+NUMSITE(a)-1)
         AREA(s):=a
       ct+= NUMSITE(a)
     end-do

     forall(c in CONTR, s in SITES | PRICE(s,c) > 0.01) create(clean(c,s))

    ! Objective: Minimize total cost of all cleaning contracts
     Cost:= sum(c in CONTR, s in SITES) PRICE(s,c)*clean(c,s)

    ! Each site must be cleaned by exactly one contractor
     forall(s in SITES) sum(c in CONTR) clean(c,s) = 1

    ! Ban same contractor from serving adjacent areas
     forall(c in CONTR, a,b in AREAS | a > b and ADJACENT(a,b) = 1)
       alloc(c,a) + alloc(c,b) <= 1

    ! Specify lower & upper limits on contracts per area
     forall(a in AREAS | LOWCON(a)>0)
       sum(c in CONTR) alloc(c,a) >= LOWCON(a)
     forall(a in AREAS | UPPCON(a)<NCONTRACTORS)
       sum(c in CONTR) alloc(c,a) <= UPPCON(a)

    ! Define alloc(c,a) to be 1 iff some clean(c,s)=1 for sites s in area a
     forall(c in CONTR, a in AREAS) do
       alloc(c,a) <= sum(s in SITES| AREA(s)=a) clean(c,s)
       alloc(c,a) >= 1.0/NUMSITE(a) * sum(s in SITES| AREA(s)=a) clean(c,s)
     end-do

     forall(c in CONTR) do
       forall(s in SITES) clean(c,s) is_binary
       forall(a in AREAS) alloc(c,a) is_binary
     end-do

     minimize(Cost)                  ! Solve the MIP problem
    end-model

In the preceding model, we have chosen to implement the constraints that force the variables allocca to become 1 whenever a variable cleancs is 1 for some site s in area a in an aggregated way (this type of constraint is usually referred to as Multiple Variable Lower Bound, MVLB, constraints). Instead of

     forall(c in CONTR, a in AREAS)
       alloc(c,a) >= 1.0/NUMSITE(a) * sum(s in SITES| AREA(s)=a) clean(c,s)

we could also have used the stronger formulation

     forall(c in CONTR, s in SITES)
       alloc(c,AREA(s)) >= clean(c,s)

but this considerably increases the total number of constraints. The aggregated constraints are sufficient to express this problem, but this formulation is very loose, with the consequence that the solution of the LP relaxation only provides a very bad approximation of the integer solution that we want to obtain. For large data sets the Branch-and-Bound search may therefore take a long time.

### Cut-and-Branch

To improve this situation without blindly adding many unnecessary constraints, we implement a cut generation loop at the top node of the search that only adds those constraints that are violated be the current LP solution.

The cut generation loop (procedure top_cut_gen) performs the following steps:

- solve the LP and save the basis
- get the solution values
- identify violated constraints and add them to the problem
- load the modified problem and load the previous basis

     procedure top_cut_gen
      declarations
        MAXCUTS = 2500              ! Max no. of constraints added in total
        MAXPCUTS = 1000             ! Max no. of constraints added per pass
        MAXPASS  = 50               ! Max no. passes
        ncut, npass, npcut: integer       ! Counters for cuts and passes
        feastol: real                     ! Zero tolerance
        solc: array(CONTR,SITES) of real  ! Sol. values for variables `clean'
        sola: array(CONTR,AREAS) of real  ! Sol. values for variables `alloc'
        objval,starttime: real
        cut: array(range) of linctr
        bas: basis                        ! LP basis
      end-declarations

      starttime:=gettime
      setparam("XPRS_CUTSTRATEGY", 0)     ! Disable automatic cuts
      setparam("XPRS_PRESOLVE", 0)        ! Switch presolve off
      feastol:= getparam("XPRS_FEASTOL")  ! Get the solver zero tolerance
      setparam("ZEROTOL", feastol * 10)   ! Set the comparison tolerance of Mosel
      ncut:=0
      npass:=0

      while (ncut<MAXCUTS and npass<MAXPASS) do
        npass+=1
        npcut:= 0
        minimize(XPRS_LIN, Cost)          ! Solve the LP
        if npass>1 and objval=getobjval: break
        savebasis(bas)                    ! Save the current basis
        objval:= getobjval                ! Get the objective value

        forall(c in CONTR) do             ! Get the solution values
          forall(a in AREAS) sola(c,a):=getsol(alloc(c,a))
          forall(s in SITES) solc(c,s):=getsol(clean(c,s))
        end-do

    ! Search for violated constraints and add them to the problem:
        forall(c in CONTR, s in SITES)
          if solc(c,s) > sola(c,AREA(s)) then
            cut(ncut):= alloc(c,AREA(s)) >= clean(c,s)
            ncut+=1
            npcut+=1
            if npcut>MAXPCUTS or ncut>MAXCUTS: break 2
          end-if

        writeln("Pass ", npass, " (", gettime-starttime, " sec), objective value ",
                objval, ", cuts added: ", npcut, " (total ", ncut,")")

        if npcut=0 then
          break
        else
          loadprob(Cost)                  ! Reload the problem
          loadbasis(bas)                  ! Load the saved basis
        end-if
      end-do
                                          ! Display cut generation status
      write("Cut phase completed: ")
      if (ncut >= MAXCUTS) then writeln("space for cuts exhausted")
      elif (npass >= MAXPASS) then writeln("maximum number of passes reached")
      else writeln("no more violations or no improvement to objective")
      end-if
     end-procedure

Assuming we add the definition of procedure top_cut_gen to the end of our model, we need to add its declaration at the beginning of the model:

    forward procedure topcutgen

and the call to this function immediately before the optimization:

     top_cut_gen                         ! Constraint generation at top node
     minimize(Cost)                      ! Solve the MIP problem

Since we wish to use our own cut strategy, we switch off the default cut generation in Xpress Optimizer:

     setparam("XPRS_CUTSTRATEGY", 0)

We also turn the presolve off since we wish to access the solution to the original problem after solving the LP-relaxations:

     setparam("XPRS_PRESOLVE", 0)

### Comparison tolerance

In addition to the parameter settings we also retrieve the feasibility tolerance used by Xpress Optimizer: the solver works with tolerance values for integer feasibility and solution feasibility that are typically of the order of 10-6 by default. When evaluating a solution, for instance by performing comparisons, it is important to take into account these tolerances.

After retrieving the feasibility tolerance of the solver we set the comparison tolerance of Mosel (ZEROTOL) to this value. This means, for example, the test x = 0 evaluates to true if x lies between -ZEROTOL and ZEROTOL, x ≤ 0 is true if the value of x is at most ZEROTOL, and x > 0 is fulfilled if x is greater than ZEROTOL.

Comparisons in Mosel always use a tolerance, with a very small default value. By resetting this parameter to the solver feasibility tolerance Mosel evaluates solution values just like Xpress Optimizer.

### Branch-and-Cut

The cut generation loop presented in the previous subsection only generates violated inqualities at the top node before entering the Branch-and-Bound search and adds them to the problem in the form of additional constraints. We may do the same using the cut manager of Xpress Optimizer. In this case, the violated constraints are added to the problem via the cut pool. We may even generate and add cuts during the Branch-and-Bound search. A cut added at a node using addcuts only applies to this node and its descendants, so one may use this functionality to define local cuts (however, in our example, all generated cuts are valid globally).

The cut manager is set up with a call to procedure tree_cut_gen before starting the optimization (preceded by the declaration of the procedure using forward earlier in the program). To avoid initializing the solution arrays and the feasibility tolerance repeatedly, we now turn these into globally defined objects:

     declarations
       feastol: real                        ! Zero tolerance
       solc: array(CONTR,SITES) of real     ! Sol. values for variables `clean'
       sola: array(CONTR,AREAS) of real     ! Sol. values for variables `alloc'
     end-declarations

     tree_cut_gen                           ! Set up cut generation in B&B tree
     minimize(Cost)                         ! Solve the MIP problem

As we have seen before, procedure tree_cut_gen disables the default cut generation and turns presolve off. It also indicates the number of extra rows to be reserved in the matrix for the cuts we are generating:

     procedure tree_cut_gen
       setparam("XPRS_CUTSTRATEGY", 0)      ! Disable automatic cuts
       setparam("XPRS_PRESOLVE", 0)         ! Switch presolve off
       setparam("XPRS_EXTRAROWS", 5000)     ! Reserve extra rows in matrix

       feastol:= getparam("XPRS_FEASTOL")   ! Get the zero tolerance
       setparam("zerotol", feastol * 10)    ! Set the comparison tolerance of Mosel

       setcallback(XPRS_CB_OPTNODE, ->cb_node)
     end-procedure

The last line of this procedure defines the optimal node (optnode) callback function that will be called by the optimizer from every node of the Branch-and-Bound search tree if the node LP relaxation has been solved to optimality. This cut generation routine (function cb_node) performs the following steps:

- get the solution values
- identify violated inequalities and add them to the problem

It is implemented as follows (we restrict the generation of cuts to the first three levels, i.e. depth < 4, of the search tree):

     function cb_node:boolean

      declarations
        ncut: integer                           ! Counters for cuts
        cut: dynamic array(range) of linctr     ! Cuts
        cutid: dynamic array(range) of integer  ! Cut type identification
        type: dynamic array(range) of integer   ! Cut constraint type
      end-declarations

      depth:=getparam("XPRS_NODEDEPTH")
      node:=getparam("XPRS_NODES")

      if depth<4 then
        ncut:=0

      ! Get the solution values
        forall(c in CONTR) do
          forall(a in AREAS) sola(c,a):=getsol(alloc(c,a))
          forall(s in SITES) solc(c,s):=getsol(clean(c,s))
        end-do

      ! Search for violated constraints
        forall(c in CONTR, s in SITES | solc(c,s) > sola(c,AREA(s))) do
          cut(ncut):= alloc(c,AREA(s)) - clean(c,s)
          cutid(ncut):= 1
          type(ncut):= CT_GEQ
          ncut+=1
        end-do

      ! Add cuts to the problem
        if ncut>0 then
          addcuts(cutid, type, cut);
          writeln("Cuts added : ", ncut, " (depth ", depth, ", node ", node,
                  ", obj. ", getparam("XPRS_LPOBJVAL"), ")")
        end-if
      end-if

     end-function

The prototype of this function is prescribed by the type of the callback (see the Xpress Optimizer Reference Manual and the documentation of setcallback in the chapter on mmxprs in the Mosel Language Reference Manual). At every node this function is called repeatedly, followed by a re-solution of the current LP.

An alternative syntax for defining callbacks relies on the subroutine name: with this syntax the function implementing the callback must be declared as public in order to make its name accessible to the external program that invokes the callback:

      forward public function cb_node:boolean
     ...
      setcallback(XPRS_CB_OPTNODE, "cb_node")

## Column generation

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

### 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,..., 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. This type of problem is called a cutting stock problem.

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 ≤ ceil(DEMANDj / PATTERNSjj), usej ∈ ℕ

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 columngen
     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 shownewpat(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)

     columngen                               ! 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 columngen 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 usej 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 columngen
       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
           shownewpat(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 columngen calls the following auxiliary function knapsack to solve an integer knapsack problem of the form

maximize z =

∑
j ∈ WIDTHS

Ci · xj

∑
j ∈ WIDTHS

Aj · xj ≤ B
∀ j ∈ WIDTHS: xj integer
∀ j ∈ WIDTHS: xj ≤ Dj

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 shownewpat to print every new pattern we find.

     procedure shownewpat(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:

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

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

Parent Topic Advanced language features

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