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 stragies 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
$ontext

This multiknapsack problem illustrates the use of user supplied cutting planes
in the GAMS BCH (branch-and-vut-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://www.ms.ic.ac.uk/info.html
First example in mknap1 which comes from

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

J.Linderoth "IE418 Lecture Notes - Lecture #19", Lehigh University, 2003.
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;

option mip=cplex;
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, solprint=off;

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 = 1 or kn.modelstat = 8) 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

GAMS Rev 137  Copyright (C) 1987-2003 GAMS Development. All rights reserved
Licensee: GAMS Development Corporation, Washington, DC   G871201:0000XX-XXX
          Free Demo,  202-342-0180,  sales@gams.com,  www.gams.com   DC9999
--- Starting compilation
--- mk.gms(68) 3 Mb
--- Starting execution
--- mk.gms(57) 4 Mb
--- Generating model m
--- mk.gms(68) 4 Mb
---    11 rows, 17 columns, and 68 non-zeroes.
--- Executing CPLEX

GAMS/Cplex    BETA 12Dec03 WIN.CP.CP 21.3 025.027.041.VIS For Cplex 9.0
Cplex 9.0.0, GAMS Link 25a6


User supplied options:
usercutcall mkc.inc mip=cplex
cuts no
preind 0
heurfreq -1
mipinterval 1

Reading data...
Starting Cplex...
Clique table members: 5
MIP emphasis: balance optimality and feasibility
Initializing dual steep norms . . .

Iteration      Dual Objective            In Variable           Out Variable
     1            5063.636364              slack(r7)              slack(r1)
     2            4325.000000                  x(c4)                  x(c5)
     3            4134.074074                  x(c5)              slack(r5)
Root relaxation solution time =    0.05 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
                  4090.6542     2                    User:  2        4
*** Calling cut generator. Added     2 cuts
                  4001.4270     4                    User:  2        5
*** Calling cut generator. Added     2 cuts
                  3950.4202     4                    User:  2        6
*** Calling cut generator. Added     1 cut
                  3871.4286     2                    User:  1        7
*** Calling cut generator. Added     1 cut
                  3841.6244     6                    User:  1        9
*** Calling cut generator. No cuts found
      1     1     3800.0000     3                   3800.0000       10
*** Calling cut generator. Added     1 cut
                                                     User:  1
*** Calling cut generator. No cuts found
*     2     1                   0     3800.0000     3800.0000       12    0.00%
*** Calling cut generator. No cuts found

User cuts applied: 6
Fixing integer variables, and solving final LP...
Using devex.

Iteration           Objective            In Variable           Out Variable
     1            3800.000000                      z                  x(c1)
     2            3800.000000              slack(r8)                  x(c6)
     3            3800.000000             slack(r10)                  x(c3)
     4            3800.000000              slack(r9)           mk(r9) artif

Proven optimal solution.

MIP Solution:         3800.000000    (12 iterations, 2 nodes)
Final Solve:          3800.000000    (4 iterations)

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

--- Restarting execution
--- mk.gms(68) 0 Mb
--- Reading solution for model m
*** Status: Normal completion

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.047 SECONDS    3.9 Mb  WIN213-137 Dec 17, 2003


EXECUTION TIME       =        0.078 SECONDS    3.9 Mb  WIN213-137 Dec 17, 2003

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

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

**** SOLVER STATUS     1 NORMAL COMPLETION         
**** MODEL STATUS      1 OPTIMAL                   
**** OBJECTIVE VALUE             3800.0000

 RESOURCE USAGE, LIMIT          1.281      1000.000
 ITERATION COUNT, LIMIT        16         10000

GAMS/Cplex    BETA 12Dec03 WIN.CP.CP 21.3 025.027.041.VIS For Cplex 9.0
Cplex 9.0.0, GAMS Link 25a6 


User supplied options: 
usercutcall bchcover.inc mip=cplex
cuts no
preind 0
heurfreq -1
mipinterval 1

*** 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
*** Calling cut generator. =2
No cuts found
Proven optimal solution.

MIP Solution:         3800.000000    (12 iterations, 2 nodes)
Final Solve:          3800.000000    (4 iterations)

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


---- EQU mk  

           LOWER          LEVEL          UPPER         MARGINAL

r1         80.0000        80.0000        80.0000         EPS         
r2         96.0000        96.0000        96.0000         EPS         
r3         20.0000        20.0000        20.0000         EPS         
r4         36.0000        36.0000        36.0000         EPS         
r5         44.0000        44.0000        44.0000         EPS         
r6         48.0000        48.0000        48.0000         EPS         
r7         10.0000        10.0000        10.0000         EPS         
r8         18.0000        18.0000        18.0000         EPS         
r9         22.0000        22.0000        22.0000         EPS         
r10        24.0000        24.0000        24.0000         EPS         

                           LOWER          LEVEL          UPPER         MARGINAL

---- EQU defobj              .              .              .             EPS         

---- VAR x  

          LOWER          LEVEL          UPPER         MARGINAL

c1          .              .             1.0000       100.0000      
c2          .             1.0000         1.0000       600.0000      
c3          .             1.0000         1.0000      1200.0000      
c4          .              .             1.0000      2400.0000      
c5          .              .             1.0000       500.0000      
c6          .             1.0000         1.0000      2000.0000      

---- VAR slack  

           LOWER          LEVEL          UPPER         MARGINAL

r1           .            14.0000        +INF             .          
r2           .            30.0000        +INF             .          
r3           .             6.0000        +INF             .          
r4           .             6.0000        +INF             .          
r5           .             3.0000        +INF             .          
r6           .             7.0000        +INF             .          
r7           .            10.0000        +INF             .          
r8           .            14.0000        +INF             .          
r9           .            12.0000        +INF             .          
r10          .            14.0000        +INF             .          

                           LOWER          LEVEL          UPPER         MARGINAL

---- VAR z                 -INF         3800.0000        +INF             .          

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