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 LOWCON_a and UPPCON_a contractors per area a.

Model formulation

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

clean_cs indicates whether contractor c is cleaning site s
alloc_ca 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 PRICE_sc is the price per site and contractor):

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

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

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

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

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, y_ca is 1 if and only if some x_cs (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 NUMSITE_a the number of sites in area a):

Implementation

The resulting Mosel program is the following model clean.mos. The variables clean_cs are defined as a dynamic array and are only created if contractor c bids for site s (that is, PRICE_sc>0 or, taking into account inaccuracies in the data, PRICE_sc>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 alloc_ca to become 1 whenever a variable clean_cs 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) then break; end-if
    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) then break 2; end-if
     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_CUTMGR, "cb_node")
 end-procedure

The last line of this procedure defines the cut manager entry callback function that will be called by the optimizer from every node of the Branch-and-Bound search tree. 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):

 public 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

  returned:=false                     ! Call this function once per node

  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)
    if solc(c,s) > sola(c,AREA(s)) then
     cut(ncut):= alloc(c,AREA(s)) - clean(c,s)
     cutid(ncut):= 1
     type(ncut):= CT_GEQ
     ncut+=1
    end-if

! Add cuts to the problem
   if ncut>0 then
    returned:=true                     ! Call this function again
    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 chapter on mmxprs in the Mosel Language Reference Manual). We declare the function as public to make sure that our model continues to work if it is compiled with the -s (strip) option. At every node this function is called repeatedly, followed by a re-solution of the current LP, as long as it returns true.