Initializing help system before first use

Column generation for a cutting stock problem: solving different models in sequence


Type: Cutting stock
Rating: 4 (medium-difficult)
Description: The first version (master model: paperp.mos, submodel: knapsack.mos using the 'bin' and 'shmem' I/O drivers, master model: paperpr.mos, submodel: knapsackr.mos, using the 'raw' and 'shmem' I/O drivers) shows how to implement the column generation 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
The second version (papersn.mos) works with multiple problems within a single model, recreating a new problem for each subproblem solving iteration.
File(s): paperp.mos, knapsack.mos (submodel), paperpr.mos, knapsackr.mos (submodel), papersn.mos


paperp.mos
(!*******************************************************
   Mosel User Guide Examples
   =========================

   file paperp.mos
   ``````````````
   Column generation algorithm for finding (near-)optimal
   solutions for a paper cutting example.
   
   Defining several (parallel) models.
   
   (c) 2008 Fair Isaac Corporation
       author: S. Heipcke, 2004, rev. May 2019
*******************************************************!)

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

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

(!
The column generation algorithm requires the solution of a knapsack problem.
This problem is defined in a separate model file.
!)
 
 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]

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

 fdelete("knapsack.bim")               ! Cleaning up
 fdelete("shmem:indata");  fdelete("shmem:resdata")

!************************************************************************
!  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

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

    write("Pass ", npass, ": ")   
    if zbest = 0 then                  ! Same as zbest <= EPS
      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.
                                       ! (useful when PRESOLVE is off)

      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

!************************************************************************
!  Solve the integer knapsack problem                              
!    z = max{cx : ax<=b, x<=d, x in Z^N}                                     
!  with b=MAXWIDTH, N=NWIDTHS, d=DEMAND
! 
! All data is exchanged between the two models via shared memory.
!************************************************************************
 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

!************************************************************************
 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


knapsack.mos
(!*******************************************************
   Mosel User Guide Examples
   =========================

   file knapsack.mos
   `````````````````
   Knapsack (sub)model of the for a paper cutting example,
   reading data from shared memory.

   *** Can be run standalone - or run from paperp.mos ***

   Solve the integer knapsack problem                              
     z = max{Cx : Ax<=B, x<=D, x in Z^N}  

   (c) 2008 Fair Isaac Corporation
       author: S. Heipcke, 2004, rev. Aug. 2011
*******************************************************!)

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)

! Integrality condition and bounds
 forall(j in WIDTHS) x(j) is_integer
 forall(j in WIDTHS) x(j) <= D(j)    ! These bounds can be omitted

 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

paperpr.mos
(!*******************************************************
   Mosel User Guide Examples
   =========================

   file paperpr.mos
   ````````````````
   Column generation algorithm for finding (near-)optimal
   solutions for a paper cutting example.
   
   Defining several (parallel) models.
   Communication of data from/to master model
   using 'shmem' IO driver with 'raw'.
  
   (c) 2008 Fair Isaac Corporation
       author: S. Heipcke, 2004, rev. May 2019
*******************************************************!)

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

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

(!
The column generation algorithm requires the solution of a knapsack problem.
This problem is defined in a separate model file.
!)
 
 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]

                                       ! 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)))
 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)  
setparam("XPRS_verbose",true)
 res:= compile("knapsackr.mos")        ! Compile the knapsack model
 load(Knapsack, "knapsackr.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   

 fdelete("knapsackr.bim")              ! Cleaning up

!************************************************************************
!  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

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

    write("Pass ", npass, ": ")   
    if zbest = 0 then                  ! Same as zbest <= EPS
      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.
                                       ! (useful when PRESOLVE is off)

      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

!************************************************************************
!  Solve the integer knapsack problem                              
!    z = max{cx : ax<=b, x<=d, x in Z^N}                                     
!  with b=MAXWIDTH, N=NWIDTHS, d=DEMAND
! 
! All data is exchanged between the two models via shared memory.
!************************************************************************
 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

  initializations to "raw:noindex"
   A as "shmem:A" B as "shmem:B" C as "shmem:C" D as "shmem:D"
  end-initializations
 
  run(Knapsack, "NWIDTHS="+NWIDTHS)    ! Start solving knapsack subproblem
  wait                                 ! Wait until subproblem finishes
  dropnextevent                        ! Ignore termination message
 
  initializations from "raw:"
   xbest as "shmem:xbest" returned as "shmem:zbest"
  end-initializations
    
 end-function

!************************************************************************
 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


knapsackr.mos
(!*******************************************************
   Mosel User Guide Examples
   =========================

   file knapsackr.mos
   ``````````````````
   Knapsack (sub)model of the for a paper cutting example,
   reading data from shared memory.
   Communication of data from/to master model
   using 'shmem' IO driver with 'raw'.

   *** Not intended to be run standalone - run from paperpr.mos ***

   Solve the integer knapsack problem                              
     z = max{Cx : Ax<=B, x<=D, x in Z^N}  

   (c) 2008 Fair Isaac Corporation
       author: S. Heipcke, 2004, rev. July 2010
*******************************************************!)

model "Knapsack"
 uses "mmxprs"

 parameters
  NWIDTHS=5                            ! Number of different widths
 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 "raw:noindex"
  A as "shmem:A" B as "shmem:B" C as "shmem:C"  D as "shmem: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)

! Integrality condition and bounds
 forall(j in WIDTHS) x(j) is_integer
 forall(j in WIDTHS) x(j) <= D(j)    ! These bounds can be omitted

 maximize(KnapObj)
 z:=getobjval
 forall(j in WIDTHS) xbest(j):=round(getsol(x(j)))

 initializations to "raw:"
  xbest as "shmem:xbest" z as "shmem:zbest"
 end-initializations

end-model

papersn.mos
(!*******************************************************
   Mosel User Guide Examples
   =========================

   file papersn.mos
   ````````````````
   Column generation algorithm for finding (near-)optimal
   solutions for a paper cutting example.
   
   - Working with multiple problems within a single model,
     recreating a new problem -
   
   (c) 2008 Fair Isaac Corporation
       author: S. Heipcke, Nov. 2008, rev. May 2019
*******************************************************!)

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): 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).
!)
 
 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
 end-declarations

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

 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   

!************************************************************************
!  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

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

    write("Pass ", npass, ": ")   
    if zbest = 0 then                  ! Same as zbest <= EPS
      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.
                                       ! (useful when PRESOLVE is off)

      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

!************************************************************************
!  Solve the integer knapsack problem                              
!    z = max{cx : ax<=b, x<=d, x in Z^N}                                     
!  with b=MAXWIDTH, N=NWIDTHS, d=DEMAND
!************************************************************************
 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)     ! These bounds can be omitted
   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

!************************************************************************

 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