FICO
FICO Xpress Optimization Examples Repository
FICO Optimization Community FICO Xpress Optimization Home
Back to examples browserPrevious exampleNext example

Cut generation for an economic lot-sizing (ELS) problem

Description
This model implements various forms of cut-and-branch and branch-and-cut algorithms. In its simplest form (looping over LP solving) it illustrates the following features:
  • adding new constraints and resolving the LP-problem (cut-and-branch)
  • basis in- and output
  • if statement
  • repeat-until statement
  • procedure
The model els.mos also implements a configurable cutting plane algorithm:
  • defining the cut manager node callback function,
  • defining and adding cuts during the MIP search (branch-and-cut), and
  • using run-time parameters to configure the solution algorithm.
The version elsglobal.mos shows how to implement global cuts. And the model version elscb.mos defines additional callbacks for extended logging and user stopping criteria based on the MIP gap.

Another implementation (main model: runels.mos, submodel: elsp.mos) parallelizes the execution of several model instances, showing the following features:
  • parallel execution of submodels
  • communication between different models (for bound updates on the objective function)
  • sending and receiving events
  • stopping submodels
The fourth implementation (main model: runelsd.mos, submodel: elsd.mos) is an extension of the parallel version in which the solve of each submodels are distributed to various computing nodes.

Further explanation of this example: elscb.mos, elsglobal.mos, runels.mos: Xpress Whitepaper 'Multiple models and parallel solving with Mosel', Section 'Solving several model instances in parallel'.


Source Files

Data Files





els.mos

(!*******************************************************
  * Mosel Example Problems                              *
  * ======================                              *
  *                                                     *
  * file els.mos                                        *
  * ````````````                                        *
  * Example for the use of the Mosel language           *
  * (Economic lot sizing, ELS, problem, solved by       *
  *  adding(l,S)-inequalities in one or several rounds  *
  *  looping over the root node or at tree nodes)       *
  *                                                     *
  * ELS considers production planning over a horizon    *
  * of T periods. In period t, t=1,...,T, there is a    *
  * given demand DEMAND[p,t] that must be satisfied by  *
  * production produce[p,t] in period t and by          *
  * inventory carried over from previous periods.       *
  * A set-up cost SETUPCOST[t] is associated with       *
  * production in period t, and the total production    *
  * capacity per period is limited. The unit production *
  * cost in period t is PRODCOST[p,t]. There is no      *
  * inventory or stock-holding cost.                    *
  * The model implements a configurable cutting plane   *
  * algorithm for this problem.                         *
  *                                                     *
  * (c) 2008 Fair Isaac Corporation                     *
  *     author: S. Heipcke, 2001, rev. July 2023        *
  *******************************************************!)

model Els                             ! Start a new model
 uses "mmxprs","mmsystem"             ! Load the optimizer library

 parameters
  ALG = 6                             ! Algorithm choice [default settings: 0]
  CUTDEPTH = 10                       ! Maximum tree depth for cut generation
  EPS = 1e-6                          ! Zero tolerance
 end-parameters 

 forward procedure lp_cut_gen         ! Declare a procedure that is
                                      ! defined later
 forward function cb_node:boolean
 forward procedure tree_cut_gen

 declarations
  TIMES = 1..15                              ! Range of time
  PRODUCTS = 1..4                            ! Set of products

  DEMAND: array(PRODUCTS,TIMES) of integer   ! Demand per period
  SETUPCOST: array(TIMES) of integer         ! Setup cost per period
  PRODCOST: array(PRODUCTS,TIMES) of integer ! Production cost per period
  CAP: array(TIMES) of integer               ! Production capacity per period
  D: array(PRODUCTS,TIMES,TIMES) of integer  ! Total demand in periods t1 - t2

  produce: array(PRODUCTS,TIMES) of mpvar    ! Production in period t
  setup: array(PRODUCTS,TIMES) of mpvar      ! Setup in period t

  solprod: array(PRODUCTS,TIMES) of real     ! Sol. values for var.s produce
  solsetup: array(PRODUCTS,TIMES) of real    ! Sol. values for var.s setup
  starttime: real
 end-declarations

 initializations from "Data/els.dat"
  DEMAND SETUPCOST PRODCOST CAP
 end-initializations

 forall(p in PRODUCTS,s,t in TIMES) D(p,s,t):= sum(k in s..t) DEMAND(p,k)

! Objective: minimize total cost
 MinCost:= sum(t in TIMES) (SETUPCOST(t) * sum(p in PRODUCTS) setup(p,t) + 
                            sum(p in PRODUCTS) PRODCOST(p,t) * produce(p,t) )

! Satisfy the total demand
 forall(p in PRODUCTS,t in TIMES) 
   Dem(p,t):= sum(s in 1..t) produce(p,s) >= sum (s in 1..t) DEMAND(p,s)

! If there is production during t then there is a setup in t
 forall(p in PRODUCTS, t in TIMES) 
  ProdSetup(p,t):= produce(p,t) <= D(p,t,getlast(TIMES)) * setup(p,t)

! Capacity limits
 forall(t in TIMES) Capacity(t):= sum(p in PRODUCTS) produce(p,t) <= CAP(t)

! Variables setup are 0/1
 forall(p in PRODUCTS, t in TIMES) setup(p,t) is_binary 

! Uncomment to get detailed MIP output
 setparam("XPRS_VERBOSE", true)
 
! All cost data are integer, we therefore only need to search for integer
! solutions
 setparam("XPRS_MIPADDCUTOFF", -0.999)

! Set Mosel comparison tolerance to a sufficiently small value
 setparam("ZEROTOL", EPS/100)
 
 writeln("**************ALG=",ALG,"***************")

 starttime:=gettime
 SEVERALROUNDS:=false; TOPONLY:=false

 case ALG of
  1: setparam("XPRS_CUTSTRATEGY", 0)  ! No cuts
  2: setparam("XPRS_PRESOLVE", 0)     ! No presolve
  3: tree_cut_gen                     ! User branch-and-cut + automatic cuts
  4: do                               ! User branch-and-cut (several rounds),
      tree_cut_gen                    ! no automatic cuts
      setparam("XPRS_CUTSTRATEGY", 0)
      SEVERALROUNDS:=true
     end-do
  5: do                               ! User cut-and-branch (several rounds)
      tree_cut_gen                    ! + automatic cuts
      SEVERALROUNDS:=true
      TOPONLY:=true
     end-do
  6: do                               ! User branch-and-cut (several rounds)
      tree_cut_gen                    ! + automatic cuts
      SEVERALROUNDS:=true
     end-do
  7: lp_cut_gen                       ! User cut-and-branch using loop
                                      ! around LP solving and save/load basis
 end-case

 minimize(MinCost)                    ! Solve the problem
                                       
 writeln("Time: ", gettime-starttime, "sec,  Nodes: ", getparam("XPRS_NODES"),
         ",  Solution: ", getobjval) 
 writeln(" "*15, "|    Prod. Quantity (Demand)")
 write("Period  Setup  |   ")
 forall(p in PRODUCTS) write(strfmt(p,-8))
 write("\n", "-"*50)
 forall(t in TIMES) do
  write("\n ", strfmt(t,2), strfmt(getsol(sum(p in PRODUCTS) setup(p,t)),8), "    |")
  forall(p in PRODUCTS) 
   write(strfmt(produce(p,t).sol,4,0), " (", DEMAND(p,t), ")")
 end-do
 writeln 

!*************************************************************************
!  Cut generation loop at the top node:
!    solve the LP and save the basis
!    get the solution values
!    identify and set up violated constraints
!    load the modified problem and load the saved basis
!*************************************************************************

 procedure lp_cut_gen
  declarations
   ncut,npass,npcut:integer            ! Counters for cuts and passes
   objval,ds: real
   cut: array(range) of linctr
   bas: basis
  end-declarations

  setparam("XPRS_CUTSTRATEGY", 0)      ! Disable automatic cuts (optional)
  setparam("XPRS_PRESOLVE", 0)         ! Switch presolve off (required)
  ncut:=0 
  npass:=0

  repeat
    npass+=1
    npcut:= 0
    minimize(XPRS_LIN+XPRS_PRI, MinCost) ! Solve the LP using primal simplex
    savebasis(bas)                       ! Save the current basis
    objval:= getobjval                   ! Get the objective value

    forall(p in PRODUCTS,t in TIMES) do  ! Get the solution values
      solprod(p,t):=produce(p,t).sol
      solsetup(p,t):=setup(p,t).sol
    end-do
  
      ! Search for violated constraints:
    forall(p in PRODUCTS,l in TIMES) do
      ds:=0 
      forall(t in 1..l)
        if(solprod(p,t) < D(p,t,l)*solsetup(p,t) + EPS) then ds += solprod(p,t)
        else  ds += D(p,t,l)*solsetup(p,t)
        end-if
  
      ! Add the violated inequality: the minimum of the actual production
      ! prod(p,t) and the maximum potential production D(p,t,l)*setup(p,t)
      ! in periods 1 to l must at least equal the total demand in periods
      ! 1 to l.
    
      if(ds < D(p,1,l) - EPS) then
        ncut+=1
        npcut+=1
        cut(ncut):= sum(t in 1..l) 
         if(solprod(p,t)<(D(p,t,l)*solsetup(p,t))+EPS, produce(p,t),
	    D(p,t,l)*setup(p,t)) >= D(p,1,l)
      end-if   
    end-do
   
    writeln("Pass ", npass, " (", gettime-starttime, " sec), objective value ", 
            objval, ", cuts added: ", npcut, " (total ", ncut,")")

    if(npcut=0) then
      writeln("Optimal integer solution found:")
    else
      loadprob(MinCost)                ! Reload the problem
      loadbasis(bas)                   ! Load the saved basis
    end-if
  until (npcut<=0 or npass>15)

 end-procedure

!*************************************************************************
!  Cut generation loop:
!    get the solution values
!    identify and set up violated constraints
!    load the modified problem and load the saved basis
!*************************************************************************

 function cb_node:boolean
  declarations
   ncut:integer                        ! Counter for cuts
   cut: array(range) of linctr         ! Cuts
   cutid: array(range) of integer      ! Cut type identification
   type: array(range) of integer       ! Cut constraint type
   ds: real
  end-declarations

  returned:=false                      ! OPTNODE: This node is not infeasible

  depth:=getparam("XPRS_NODEDEPTH")
  cnt:=getparam("XPRS_CALLBACKCOUNT_OPTNODE")

  if ((TOPONLY and depth<1) or (not TOPONLY and depth<=CUTDEPTH)) and 
     (SEVERALROUNDS or cnt<=1) then
   ncut:=0 

 ! Get the solution values
   forall(t in TIMES, p in PRODUCTS) do
     solprod(p,t):=produce(p,t).sol
     solsetup(p,t):=setup(p,t).sol
   end-do
  
 ! Search for violated constraints
   forall(p in PRODUCTS,l in TIMES) do
    ds:=0 
    forall(t in 1..l)
      if(solprod(p,t) < D(p,t,l)*solsetup(p,t) + EPS) then ds += solprod(p,t)
      else  ds += D(p,t,l)*solsetup(p,t)
      end-if
  
   ! Generate the violated inequality
    if(ds < D(p,1,l) - EPS) then
      cut(ncut):= sum(t in 1..l) 
       if(solprod(p,t)<(D(p,t,l)*solsetup(p,t))+EPS, produce(p,t), 
          D(p,t,l)*setup(p,t)) - D(p,1,l)
      cutid(ncut):= 1
      type(ncut):= CT_GEQ
      ncut+=1
    end-if   
  end-do

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

! ****Optimizer settings for using the cut manager****

 procedure tree_cut_gen
  setparam("XPRS_PRESOLVE", 0)        ! Switch presolve off
  setparam("XPRS_EXTRAROWS", 5000)    ! Reserve extra rows in matrix
  setcallback(XPRS_CB_OPTNODE, ->cb_node)  ! Set the optnode callback function
 end-procedure

end-model

Economic Lot-Sizing (ELS)
=========================

A well-known class of valid inequalities for ELS are the
(l,S)-inequalities.  Let D(pq) denote the total demand in periods p 
to q and y(t) be a binary variable indicating whether there is any 
production in period t.  For each period l and each subset of periods S 
of 1 to l, the (l,S)-inequality is

    sum (t=1:l | t in S) x(t) + sum (t=1:l | t not in S) D(tl) * y(t)
        >= D(1l)

It says that actual production x(t) in periods included S plus maximum 
potential production D(tl)*y(t) in the remaining periods (those not in 
S) must at least equal total demand in periods 1 to l.  Note that in 
period t at most D(tl) production is required to meet demand up to 
period l.

Based on the observation that

    sum (t=1:l | t in S) x(t) + sum (t=1:l | t not in S) D(tl) * y(t)
        >= sum (t=1:l) min(x(t), D(tl) * y(t))
        >= D(1l)

it is easy to develop a separation algorithm and thus a cutting plane
algorithm based on these (l,S)-inequalities.

Back to examples browserPrevious exampleNext example