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

• adding new constraints and resolving the LP-problem (cut-and-branch)
• basis in- and output
• if statement
• repeat-until statement
• procedure

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

• parallel execution of submodels
• communication between different models (for bound updates on the objective function)
• sending and receiving events
• stopping submodels

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