BCH Facility

	[ Home \| Support \| Sales \| Solvers \| Documentation \| Model Libraries \| Search \| Contact Us ]

The GAMS Branch-and-Cut-and-Heuristic Facility

Introduction

It is a known fact that hard mixed-integer programming (MIP) problems can significantly benefit from user supplied routines that generate cutting planes and good integer feasible solutions. GAMS users could supply cutting planes and an integer feasible point as part of the model given to the solver, by adding a set of constraints representing likely to be violated cuts and a feasible solution in combination with the tryint parameter or solver specific options like mipstart in GAMS/CPLEX. A true dynamic interaction between a running branch-and-cut (B&C) solver like CPLEX and XPRESS and user supplied routines was not possible. This new GAMS Branch-and-Cut-and-Heuristic (BCH) facility closes this gap.

The GAMS BCH facility automates all major steps necessary to define, excute and control the use of user defined routines within the framework of general purpose MIP codes. This allows GAMS users to apply complex solution strategies without having to have intimate knowledge about the inner workings of a specific MIP system. BCH strategies can now be implemented rapidly and reliably within a matter of days rather than weeks.

For those of you who learn best by example, we have prepared a multi knapsack problem with cover cut generation in the Annex of this page

A presentation about BCH with performance comparisons can be found on the GAMS presentation page

Design Overview

After a node in the branch-and-bound (B&B) tree is processed by the solver a user supplied routine can be called. All available information from this node, like the LP solution, and the B&B tree, like the current incumbent, is exported in the original namespace to a GAMS database. After that the solver can call two different user supplied GAMS programs, one for finding violated cuts (cut generator) and one to find a better incumbent (heuristic). These GAMS programs must import the information from the solver, find cuts or an integer solution and export their findings to another GAMS database which is imported by the solver after the cut generator or heuristic GAMS program terminates. In case the heuristic GAMS program found a solution, the solver checks this solution for infeasibilties and in case this check is passed and the solution is better than the best known solution, the solver updates it's incumbent. In case violated cuts are found, the solver adds these cuts to it's cut pool, resolves the LP at the node and calls the user supplied cut generating GAMS model again. If no cutting planes are found, the solver will process the next node. Please note that the solver cannot perform any checks on the new cuts. Hence, it is possible to cut off areas of the feasible space and even the optimal solution with wrong cuts.

The Interface

Since data flows between the B&C solver and the user supplied cut generator and heuristic, both ends need to comply to a prespecified data contract. The GAMS databases make use of the GAMS Data eXchange (GDX) facility that allows convenient import and export of data. The B&C solver exports the LP solution, i.e. level values (.l) and the bounds (.lo and .up), to the GDX file bchout.gdx. The bounds are the actual bounds in this node. Hence for discrete variables, they reflect branching decisons made in the B&B tree. In a similar way, the incumbent solution is exported to the GDX file bchout_i.gdx. The bounds for the incumbent solution reflect global bounds, i.e. the original bounds or globally feasible bounds improved by the preprocessor of the B&C solver. The solution can be imported using the compile time $load or run time execute_load routines.

At the end of the heuristic model the newly found solution should be stored in the level values of the variables corresponing to the original variables. For example, if the original model uses binary varible open(i,t), then at the end of the heuristic open.l(i,t) should contain a 0 (zero) or a 1. If the heuristic cannot find an integer feasible solution, the model can terminate with an execution error triggered by abort to prevent the B&C solver from reading the results from the heuristic run.

Exporting violated cuts is a little more complicated because we need the complete constraints of the cuts. The B&C solver needs to know how many cuts have been found by the cut generator, it has to read the constraint matrix, the sense (=l=, =g=, or =e=), and the right hand side of these cuts. The cut generator model has to define and fill the following symbols appropriately.

$set MaxCuts 100
Set        cc          'cuts'  /1*%MaxCuts%/
Parameter  numcuts     'number of cuts to be added' /0/
           rhs_c(cc)   'cut rhs'
           sense_c(cc) 'the sense of the cuts 1 =L=, 2 =E=, 3 =G=';

The only thing missing is the constraint matrix. For each variable which is part of a cut, we need a parameter in the cut generator model. For example, variable open(i,t) is part of a cut. We need to define

Parameter  open_c(cc,i,t) 'coefficients of open in cut cc';

This parameter has an extra cut index in the first position. The values of this parameter should reflect the coeffiecients in the cut constraint matrix. The B&C solver reads all parameters that end in _c, takes the basename and looks for a variable with that name and indicies and builds up the cut matrix. A cut cannot introduce a new variable into the model. All cuts added to the model are global cuts, they apply to the entire problem, not just to a subtree.

Solvers aware of BCH

With the current GAMS system (21.3) only two solvers are aware of the BCH facility: GAMS/CPLEX and GAMS/SBB. We will shortly add GAMS/XPRESS to this list. GAMS/CPLEX implements the full interface described in this document, heuristic and cut generation. GAMS/SBB only uses the heuristic.

Options

There are three different option type that come with the BCH facility. One set option helps to determine when to call the heuristic model and cut generator model. The second set allows to overwrite the filenames for the GDX files (to avoid name clashes) and the options to define the heuristic and cut generation GAMS calls. These options go like other options into the option file of the solver (e.g. cplex.opt or sbb.op9).

Name Description Default

userheurfreq Determines the frequency of the heuristic model calls. 10

userheurmult Determines the multiplier for the frequency of the heuristic model calls. 2

userheurinterval Determines the interval when to apply the multiplier for the frequency of the heuristic model calls. For example, for the first 100 (userheurinterval) nodes, the solver call every 10th (userheurfreq) node the heuristic. After 100 nodes, the frequency gets multiplied by 10 (uerheurmult), so that for the next 100 node the solver calls the heuristic every 20th node. For nodes 200-300, the heuristic get called every 40th node, for nodes 300-400 every 80th node and after node 400 every 100th node. 100

userheurfirst Calls the heuristic for the first n nodes. 10

userheurobjfirst Similar to userheurfirst but only calls the heuristic if the relaxed objective value promises a significant improvement of the current incumbent, i.e. the LP value of the node has to be closer to the best bound than the current incumbent. 0

userheurnewint Calls the heuristic if the solver found a new integer feasible solution (yes/no). no

userheurcall The GAMS command line (minus the gams executable name) to call the heuristic. no

usercutfreq Determines the frequency of the cut generator model calls. 10

usercutmult Determines the multiplier for the frequency of the cut generator model calls. 2

usercutinterval Determines the interval when to apply the multiplier for the frequency of the cut generator model calls. See userheurinterval for details. 100

usercutfirst Calls the cut generator for the first n nodes. 10

usercutnewint Calls the cut generator if the solver found a new integer feasible solution (yes/no). no

usercutcall The GAMS command line (minus the gams executable name) to call the cut generator. no

userincbcall The GAMS command line (minus the gams executable name) to call the incumbent checking routine. The incumbent is rejected if the GAMS program terminates normally. In case of a compilation or execution error, the incumbent is accepted. no

usergdxin The name of the GDX file read back into the solver. bchin.gdx

usergdxname The name of the GDX file exported from the solver with the solution at the node. bchout.gdx

usergdxnameinc The name of the GDX file exported from the solver with the incumbent solution. bchout_i.gdx

usergdxprefix Prefixes usergdxin, usergdxname, and usergdxnameinc empty

userincbcall The GAMS command line to call the incumbent checking program. empty

userincbicall The GAMS command line to call the incumbent reporting program empty

userjobid Postfixes listing and log file and adds --userjobid to the user*calls. Postfixes gdxname, gdxnameinc and gdxin empty

userkeep Calls gamskeep instead of gams 0

Name	Description	Default
userheurfreq	Determines the frequency of the heuristic model calls.	10
userheurmult	Determines the multiplier for the frequency of the heuristic model calls.	2
userheurinterval	Determines the interval when to apply the multiplier for the frequency of the heuristic model calls. For example, for the first 100 (userheurinterval) nodes, the solver call every 10th (userheurfreq) node the heuristic. After 100 nodes, the frequency gets multiplied by 10 (uerheurmult), so that for the next 100 node the solver calls the heuristic every 20th node. For nodes 200-300, the heuristic get called every 40th node, for nodes 300-400 every 80th node and after node 400 every 100th node.	100
userheurfirst	Calls the heuristic for the first n nodes.	10
userheurobjfirst	Similar to userheurfirst but only calls the heuristic if the relaxed objective value promises a significant improvement of the current incumbent, i.e. the LP value of the node has to be closer to the best bound than the current incumbent.	0
userheurnewint	Calls the heuristic if the solver found a new integer feasible solution (yes/no).	no
userheurcall	The GAMS command line (minus the gams executable name) to call the heuristic.	no
usercutfreq	Determines the frequency of the cut generator model calls.	10
usercutmult	Determines the multiplier for the frequency of the cut generator model calls.	2
usercutinterval	Determines the interval when to apply the multiplier for the frequency of the cut generator model calls. See userheurinterval for details.	100
usercutfirst	Calls the cut generator for the first n nodes.	10
usercutnewint	Calls the cut generator if the solver found a new integer feasible solution (yes/no).	no
usercutcall	The GAMS command line (minus the gams executable name) to call the cut generator.	no
userincbcall	The GAMS command line (minus the gams executable name) to call the incumbent checking routine. The incumbent is rejected if the GAMS program terminates normally. In case of a compilation or execution error, the incumbent is accepted.	no
usergdxin	The name of the GDX file read back into the solver.	bchin.gdx
usergdxname	The name of the GDX file exported from the solver with the solution at the node.	bchout.gdx
usergdxnameinc	The name of the GDX file exported from the solver with the incumbent solution.	bchout_i.gdx
usergdxprefix	Prefixes usergdxin, usergdxname, and usergdxnameinc	empty
userincbcall	The GAMS command line to call the incumbent checking program.	empty
userincbicall	The GAMS command line to call the incumbent reporting program	empty
userjobid	Postfixes listing and log file and adds --userjobid to the user*calls. Postfixes gdxname, gdxnameinc and gdxin	empty
userkeep	Calls gamskeep instead of gams	0

For the Oil Pipeline Network Design problem a typical GAMS/CPLEX option file with BCH options could look like this:

userheurcall bchoil_h.inc mip cplex optcr 0 reslim 10
userheurfirst 5
userheurfreq 20
userheurinterval 1000

usercutcall bchoil_c.inc mip cplex
usercutfirst 0
usercutfreq 0
usercutnewint yes

Examples

We have gathered a few examples that show how to use the new BCH facility:

Trimloss Optimization with SBB

This model implements a very simple LP/MIP based heuristic for the trimloss minimization problem.

Single Commodity Uncapacitated Fixed Charge Network Flow Problem

This model implements simple but difficult to separate cuts for a network design problem. The global solver BARON is used to find violated cuts by solving a nonconvex MINLP.

Oil Pipeline Design Problem

This is the most complex example. It implements three different heuristics: An initial heuristic based on a simplified cost structure, a rounding heuristic, and local branching. In addition complex cuts are generated by solving regionalized versions of the original problem.

MIP Decomposition and Parallel Grid Submission - DICE Example

This example uses many of the renameing and jobid options since we have multiple jobs running that need to work on different filenames. This example also uses the incumbent reporting call userincbicall.

Cplex Solution Pool for a Simple Facility Location Problem

This example uses the incumbent checking call userincbcall as an advanced filter for accepting/rejecting solutions found by Cplex.

Annex

The simple multi knapsack problem will be used to illustrate how to implement user cuts in the BCH facility. The model can also be found in the model library under the name bchmknap.

Base Model

$Title Multi knapsack problem using BCH Facility (BCHMKNAP,SEQ=289)
$ontext

This multiknapsack problem illustrates the use of user supplied cutting planes
in the GAMS BCH (branch-and-cut-and-heuristic) facility. Please note, that cover
cuts used in this example, are already implemented in modern MIP solvers.

Example taken from the OR-Library
http://people.brunel.ac.uk/~mastjjb/jeb/orlib/mknapinfo.html
first example in mknap1.


Beasley, J E, OR-Library: Distributing Test Problems by Electronic
Mail. Journal of the Operational Research Society 41 (11) (1990),
1069-1072.

Petersen, C C, Computational experience with variants of the Balas
algorithm applied to the selection of R&D projects. Management Science
13 (9) (1967), 736-750.

Linderoth, J, IE418 Lecture Notes - Lecture #19, 2003. Lehigh University,
http://www.lehigh.edu/~jtl3/teaching/ie418/lecture19.pdf

$offtext

set j columns /c1*c6/
    i rows    /r1*r10/

Parameters obj(j) / c1  100, c2  600, c3 1200, c4 2400, c5  500, c6 2000 /
           rhs(i) / r1   80, r2   96, r3   20, r4   36, r5   44, r6   48,
                    r7   10, r8   18, r9   22, r10  24 /;

Table a(i,j)
      c1   c2   c3   c4   c5   c6
 r1    8   12   13   64   22   41
 r2    8   12   13   75   22   41
 r3    3    6    4   18    6    4
 r4    5   10    8   32    6   12
 r5    5   13    8   42    6   20
 r6    5   13    8   48    6   20
 r7                        8
 r8    3         4         8
 r9    3    2    4         8    4
 r10   3    2    4    8    8    4

Binary   variables x(j)
Positive variables slack(i)
         Variable  z;

Equations mk(i), defobj;

defobj.. z =e= sum(j, obj(j)*x(j));

mk(i).. sum(j, a(i,j)*x(j)) + slack(i) =e= rhs(i);

model m /all/;

* Export base data
execute_unload 'mkdata', j, i, a, rhs;

$ifi %system.mip% == cplex   $goto cont
$abort 'BCH Facility not available for MIP solver %system.mip%.'
$label cont

m.optfile = 1;
m.optcr   = 0;

* We activate the user supplied cut generator and turn all advanced CPLEX options off
$onecho > cplex.opt
usercutcall bchcover.inc mip=cplex
cuts no
preind 0
heurfreq -1
mipinterval 1
$offecho

solve m max z using mip;

Cover cut Generator

$Title Simple cover inequalities for the multi-knapsack problem with integer data

* Declare and get selected data from base model
Set j, i;
Parameter a(i,j), rhs(i);

$gdxin mkdata
$load j i a rhs

* Declare selected variables from base MIP model
Binary   variables x(j)
Positive variable  slack(i);

* Get current values from MIP solver
$gdxin bchout
$load x slack

* Define separation model
Parameter ai(j), rhsi slice of the multi knapsack problem;

Binary variable y(j) membership in the cover
       variable obj;

Equations k, defobj;

defobj.. obj =e= sum(j, (1-x.l(j))*y(j));

* rhsi+1 only works if all ai are integer
k.. sum(j, ai(j)*y(j)) =g= rhsi+1;

model kn /all/;

* In order to communicate with the MIP solver we need certain conventions
* 1. Cut matrix interface requirement with fixed names
Set        cut           potential cuts  / 1*5 /
Parameters rhs_c(cut)    cut rhs
           sense_c(cut)  'the sense of the cuts 1 =L=, 2 =E=, 3 =G='
           numcuts       number of cuts passed back
* 2. Nonzeros in cut matrix using the original variable name plus _c
           x_c(cut,j)    coefficient of the x variables in the cut

* We solve one knapsack type MIP to solve the cover separation problem for each
* row in the original problem that is binding
* The actual cover cut is sum(cover(j), x(j)) =l= sum(cover(j),1)-1;

option optcr=0, limrow=0, limcol=0;
option solprint=off, solvelink=%solvelink.CallModule%;

Set c(cut) current cut /1/;
numcuts = 0;
loop(i$(numcuts < card(cut) and slack.l(i) < 1e-6),
  ai(j) = a(i,j);
  rhsi = rhs(i);
  solve kn min obj using mip;
  if ((kn.modelstat = %modelstat.Optimal% or
       kn.modelstat = %modelstat.IntegerSolution%) and obj.l < 0.95,
     numcuts    = numcuts + 1;
     x_c(c,j)   = round(y.l(j));
     rhs_c(c)   = sum(j,round(y.l(j))) - 1;
     sense_c(c) = 1;
     c(cut)     = c(cut-1);
  )
);

* One can debug the cut generator by solving the RMIP of the base model (change
* mip to rmip) and creating a GDX file with the name bchout.gdx:
*  gams bchmknap rmip=cplex gdx=bchout
* Start the cut generator and analyze the cuts generated
*  gams bchcover.inc mip=cplex
display numcuts, x_c, rhs_c, sense_c;

Solution log

--- Job bchmknap.gms Start 08/06/10 04:36:58 WEX-VS8 23.5.1 x86/MS Windows        
GAMS Rev 235  Copyright (C) 1987-2010 GAMS Development. All rights reserved
Licensee: GAMS Development Corporation, Washington, DC   G871201/0000CA-ANY
          Free Demo,  202-342-0180,  sales@gams.com,  www.gams.com   DC0000
--- Starting compilation
--- bchmknap.gms(77) 3 Mb
--- Starting execution: elapsed 0:00:00.003
--- bchmknap.gms(64) 4 Mb
--- Generating MIP model m
--- bchmknap.gms(75) 4 Mb
---   11 rows  17 columns  68 non-zeroes
---   6 discrete-columns
--- Executing CPLEX: elapsed 0:00:00.008

IBM ILOG CPLEX   Jul  4, 2010 23.5.1 WIN 18414.18495 VS8 x86/MS Windows
Cplex 12.2.0.0, GAMS Link 34 

Reading parameter(s) from "C:\tmp\cplex.opt"
>>  usercutcall bchcover.inc mip=cplex
>>  cuts no
>>  preind 0
>>  heurfreq -1
>>  mipinterval 1
Finished reading from "C:\tmp\cplex.opt"
Cplex MIP uses 1 of 8 parallel threads. Change default with option THREADS.
Reading data...
Starting Cplex...
Warning: Control callbacks may disable some MIP features.
Clique table members: 24.
MIP emphasis: balance optimality and feasibility.
MIP search method: traditional branch-and-cut.
Parallel mode: none, using 1 thread.
Initializing dual steep norms . . .

Iteration      Dual Objective            In Variable           Out Variable
     1            4400.000000                  x(c4)              slack(r6)
     2            4216.666667              slack(r7)              slack(r1)
     3            4134.074074              slack(r6)              slack(r5)
Root relaxation solution time =    0.00 sec.

        Nodes                                         Cuts/
   Node  Left     Objective  IInf  Best Integer     Best Node    ItCnt     Gap

      0     0     4134.0741     2                   4134.0741        3         
*** Calling cut generator. Added     2 cuts
      0     0     4090.6542     2                     User: 1        4         
*** Calling cut generator. Added     2 cuts
      0     0     4001.4270     4                     User: 1        5         
*** Calling cut generator. Added     2 cuts
      0     0     3950.4202     4                     User: 1        6         
*** Calling cut generator. Added     1 cut 
      0     0     3871.4286     2                     User: 1        7         
*** Calling cut generator. Added     1 cut 
      0     0     3841.6244     6                     User: 1        9         
*** Calling cut generator. No cuts found
*** Calling cut generator. No cuts found
      0     2     3841.6244     6                   3800.0000        9         
Elapsed real time =   0.58 sec. (tree size =  0.00 MB, solutions = 0)
*** Calling cut generator. Added     1 cut 
*** Calling cut generator. No cuts found
      1     3     3800.0000     3                   3800.0000       11         
                                                      User: 1                  
*** Calling cut generator. No cuts found
*     2     0      integral     0     3800.0000     3800.0000       12    0.00%

User cuts applied:  6
MIP status(101): integer optimal solution
Fixing integer variables, and solving final LP...

Iteration      Dual Objective            In Variable           Out Variable
     1 sI            0.000000                      z           defobj artif
     2            8000.000000              slack(r8)                  x(c3)
     3            3938.461538             slack(r10)              slack(r5)
     4            3800.000000              slack(r5)                  x(c2)
     5            3800.000000              slack(r9)           mk(r9) artif
Fixed MIP status(1): optimal

Proven optimal solution.

MIP Solution:         3800.000000    (12 iterations, 3 nodes)
Final Solve:          3800.000000    (5 iterations)

Best possible:        3800.000000
Absolute gap:            0.000000
Relative gap:            0.000000

--- Restarting execution
--- bchmknap.gms(75) 0 Mb
--- Reading solution for model m
*** Status: Normal completion
--- Job bchmknap.gms Stop 08/06/10 04:36:59 elapsed 0:00:00.758

Selected Portion of the Listing File

MODEL STATISTICS

BLOCKS OF EQUATIONS           2     SINGLE EQUATIONS           11
BLOCKS OF VARIABLES           3     SINGLE VARIABLES           17
NON ZERO ELEMENTS            68     DISCRETE VARIABLES          6


GENERATION TIME      =        0.000 SECONDS      4 Mb  WIN235-235 Jul  2, 2010


EXECUTION TIME       =        0.000 SECONDS      4 Mb  WIN235-235 Jul  2, 2010
GAMS Rev 235  WEX-VS8 23.5.1 x86/MS Windows             08/06/10 04:36:58 Page 5
Multi knapsack problem using BCH Facility (BCHMKNAP,SEQ=289)
Solution Report     SOLVE m Using MIP From line 75


               S O L V E      S U M M A R Y

     MODEL   m                   OBJECTIVE  z
     TYPE    MIP                 DIRECTION  MAXIMIZE
     SOLVER  CPLEX               FROM LINE  75

**** SOLVER STATUS     1 Normal Completion         
**** MODEL STATUS      1 Optimal                   
**** OBJECTIVE VALUE             3800.0000

 RESOURCE USAGE, LIMIT          0.712      1000.000
 ITERATION COUNT, LIMIT        17    2000000000

IBM ILOG CPLEX   Jul  4, 2010 23.5.1 WIN 18414.18495 VS8 x86/MS Windows
Cplex 12.2.0.0, GAMS Link 34 

Reading parameter(s) from "C:\tmp\cplex.opt"
>>  usercutcall bchcover.inc mip=cplex
>>  cuts no
>>  preind 0
>>  heurfreq -1
>>  mipinterval 1
Finished reading from "C:\tmp\cplex.opt"
Cplex MIP uses 1 of 8 parallel threads. Change default with option THREADS.
*** Calling cut generator. =2
Added     2 cuts
*** Calling cut generator. =2
Added     2 cuts
*** Calling cut generator. =2
Added     2 cuts
*** Calling cut generator. =2
Added     1 cut 
*** Calling cut generator. =2
Added     1 cut 
*** Calling cut generator. =2
No cuts found
*** Calling cut generator. =2
No cuts found
*** Calling cut generator. =2
Added     1 cut 
*** Calling cut generator. =2
No cuts found
*** Calling cut generator. =2
No cuts found
MIP status(101): integer optimal solution
Fixed MIP status(1): optimal
Proven optimal solution.

MIP Solution:         3800.000000    (12 iterations, 3 nodes)
Final Solve:          3800.000000    (5 iterations)

Best possible:        3800.000000
Absolute gap:            0.000000
Relative gap:            0.000000


---- EQU mk  

       LOWER     LEVEL     UPPER    MARGINAL

r1     80.000    80.000    80.000      EPS       
r2     96.000    96.000    96.000      EPS       
r3     20.000    20.000    20.000      EPS       
r4     36.000    36.000    36.000      EPS       
r5     44.000    44.000    44.000      EPS       
r6     48.000    48.000    48.000      EPS       
r7     10.000    10.000    10.000      EPS       
r8     18.000    18.000    18.000      EPS       
r9     22.000    22.000    22.000      EPS       
r10    24.000    24.000    24.000      EPS       

                       LOWER     LEVEL     UPPER    MARGINAL

---- EQU defobj          .         .         .        1.000      

---- VAR x  

      LOWER     LEVEL     UPPER    MARGINAL

c1      .         .        1.000   100.000      
c2      .        1.000     1.000   600.000      
c3      .        1.000     1.000  1200.000      
c4      .         .        1.000  2400.000      
c5      .         .        1.000   500.000      
c6      .        1.000     1.000  2000.000      

---- VAR slack  

       LOWER     LEVEL     UPPER    MARGINAL

r1       .       14.000     +INF       .         
r2       .       30.000     +INF       .         
r3       .        6.000     +INF       .         
r4       .        6.000     +INF       .         
r5       .        3.000     +INF       .         
r6       .        7.000     +INF       .         
r7       .       10.000     +INF       .         
r8       .       14.000     +INF       .         
r9       .       12.000     +INF       .         
r10      .       14.000     +INF       .         

                       LOWER     LEVEL     UPPER    MARGINAL

---- VAR z              -INF   3800.000     +INF       .         


**** REPORT SUMMARY :        0     NONOPT
                             0 INFEASIBLE
                             0  UNBOUNDED