(!*******************************************************
Mosel User Guide Examples
=========================
file paper.mos
``````````````
Column generation algorithm for finding (near-)optimal
solutions for a paper cutting example.
(c) 2008 Fair Isaac Corporation
authors: Bob Daniel, S. Heipcke, J. Tebboth, 2001
rev. June 2022
*******************************************************!)
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)
(!
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
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]
! 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)
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)),", ")
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 columngen
declarations
dualdem: array(WIDTHS) of real
xbest: array(WIDTHS) of integer
dw, zbest, objval: real
bas: basis
end-declarations
setparam("XPRS_verbose", true)
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 ! Same as zbest <= EPS
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.
! (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
!
! 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
setparam("XPRS_PRESOLVE", 1) ! Switch presolve on
! 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 and bounds
if pass=1 then
forall(j in WIDTHS) x(j) is_integer
forall(j in WIDTHS) x(j) <= D(j) ! These bounds can be omitted
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)
setparam("XPRS_PRESOLVE", 0) ! Switch presolve off
end-function
!************************************************************************
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
|
(!*******************************************************
Mosel User Guide Examples
=========================
file papers.mos
```````````````
Column generation algorithm for finding (near-)optimal
solutions for a paper cutting example.
- Working with multiple problems within a single model -
(c) 2008 Fair Isaac Corporation
author: S. Heipcke, Nov. 2008, rev. Jun. 2022
*******************************************************!)
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)
(!
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
Knapsack: mpproblem ! Knapsack subproblem
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]
! 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)
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)),", ")
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 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 ! Same as zbest <= EPS
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.
! (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,
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) ! These bounds can be omitted
end-if
maximize(KnapObj)
returned:=getobjval
forall(j in WIDTHS) xbest(j):=round(getsol(x(j)))
end-do
end-function
!************************************************************************
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
|
(!*******************************************************
Mosel User Guide Examples
=========================
file papers_log.mos
```````````````````
Column generation algorithm for finding (near-)optimal
solutions for a paper cutting example.
- Working with multiple problems within a single model -
- Detailed solution printout -
(c) 2010 Fair Isaac Corporation
author: S. Heipcke, July 2010, rev. Sep. 2022
*******************************************************!)
model Papermill
uses "mmxprs"
uses "mmsystem"
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)
forward procedure printallpat(ifsol:boolean)
(!
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, RP) 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: mpproblem ! Knapsack subproblem
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]
! 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
!)
writeln("Starting patterns: ")
printallpat(false)
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)
columngen ! Column generation at top node
writeln("All patterns:")
printallpat(false)
minimize(MinRolls) ! Compute the best integer solution
! for the current problem (including
! the new columns)
writeln("Best integer solution: ", getobjval, " rolls")
TotalCut:= sum(i in WIDTHS) (sum(j in RP) WIDTH(i)*PATTERNS(i,j)*use(j))
writeln("Total width of production: ", getsol(TotalCut),
" (excess width: ",
getsol(TotalCut)-sum(i in WIDTHS) WIDTH(i)*DEMAND(i),
" num.pat.: ", round(sum(i in WIDTHS) -getslack(Dem(i))), ")")
write(" Rolls per pattern: ")
forall(i in RP) write(strfmt(getsol(use(i)),5))
write("\n Cutting waste: ")
forall(j in RP)
write( strfmt((MAXWIDTH-sum(i in WIDTHS) PATTERNS(i,j)*WIDTH(i)) * use(j).sol, 5))
writeln
writeln("Detailed solution:")
printallpat(true)
!************************************************************************
! 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 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 ! Same as zbest <= EPS
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
forall(i in WIDTHS) PATTERNS(i,npatt):=xbest(i)
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,
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) ! These bounds can be omitted
end-if
maximize(KnapObj)
returned:=getobjval
forall(j in WIDTHS) xbest(j):=round(getsol(x(j)))
end-do
end-function
!************************************************************************
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
!************************************************************************
procedure printallpat(ifsol: boolean)
write("Pattern no.: ")
forall(j in RP) write(strfmt(j,4))
write(" | Demand")
if ifsol: write(" | Produced | Excess")
writeln("\n", "-"*(if(ifsol,42,22)+4*getsize(RP)))
forall(i in WIDTHS) do
write(formattext("Width %4g: ", WIDTH(i)))
forall(j in RP) write(strfmt(PATTERNS(i,j),4))
write(" |", strfmt(DEMAND(i),5))
if ifsol: write(formattext(" | %6g | %3d",
sum(j in RP) PATTERNS(i,j)*use(j).sol, round(-getslack(Dem(i)))))
writeln
end-do
writeln("-"*(if(ifsol,42,22)+4*getsize(RP)))
localsetparam("REALFMT","%4g")
write("Waste: ")
forall(j in RP) write(MAXWIDTH-sum(i in WIDTHS) PATTERNS(i,j)*WIDTH(i))
writeln("\n", "-"*(14+4*getsize(RP)))
if ifsol then
write("Times used: ")
forall(j in RP) write(use(j).sol)
writeln
end-if
end-procedure
end-model
|
