Obtaining extreme ray from unbounded problem for use with benders decomposition

[Eric B.]

I am trying to obtain an extreme rays from a dual subproblem for the purpose of adding feasibility cuts to the master problem of a Benders decomposition approach for an MIP, but the method of obtaining these rays is unclear. The type of problem (simplified) I have is:

subproblem primal: max  2x
                   s.t. 3x <= -5
                         x >=  0

subproblem dual:   min  -5y
                   s.t.  3y >= 2
                          y >= 0

It is easy to see the primal problem is infeasible and the dual problem is unbounded.

If I run (with 5 instead of -5 so both primal and dual are bounded and feasible):

import pyomo.environ
import pyomo.opt

c = 5

pm = pyomo.environ.AbstractModel()

pm.c = c
pm.x = pyomo.environ.Var(domain = pyomo.environ.NonNegativeReals)
pm.dual = pyomo.environ.Suffix(direction = pyomo.environ.Suffix.IMPORT)

def pObjFunc(pm):
    return 2*pm.x
pm.Obj = pyomo.environ.Objective(rule= pObjFunc,sense=pyomo.environ.maximize)

def pCons(pm):
    return 3*pm.x <= pm.c
pm.Const = pyomo.environ.Constraint(rule= pCons)

opt=pyomo.opt.SolverFactory('cbc')

inst=pm.create_instance()
results=opt.solve(inst)
inst.solutions.load_from(results)

print 'primal:',results.solver.termination_condition
print 'primal constraint dual:',inst.dual[inst.Const]

dm = pyomo.environ.AbstractModel()

dm.c = c
dm.y = pyomo.environ.Var(domain = pyomo.environ.NonNegativeReals)

def dObjFunc(dm):
    return dm.c*dm.y
dm.Obj = pyomo.environ.Objective(rule= dObjFunc,sense=pyomo.environ.minimize)

def dCons(pm):
    return 3*pm.y >= 2
dm.Const = pyomo.environ.Constraint(rule= dCons)

inst=dm.create_instance()
results=opt.solve(inst)
inst.solutions.load_from(results)

print 'dual:',results.solver.termination_condition
print 'dual solution:',inst.y.value

The output is:

primal: optimal
primal constraint dual: -0.66666667
dual: optimal
dual solution: 0.66666667

However when I make c = -5 (so I have the primal infeasible and dual unbounded), the output is:

primal: infeasible
(second print commented out: it causes error, since nothing is imported to suffix from solver)
dual: unbounded
dual solution: None

How would I go about getting an extreme ray associated with the dual problem? I am trying to solve an MIP with Benders decomposition and generate feasibility cuts, so there must be a way to get an unbounded ray either through solving either the infeasible primal or unbounded dual. It seems straightforward to obtain optimality cuts (since solvers return solutions), but not the feasibility cuts.

I see some mention in this thread with respect to Gurobi: https://groups.google.com/forum/#!searchin/pyomo-forum/unbounded$20ray%7Csort:date/pyomo-forum/GLIwAZtMeTQ/680CUZ3PBQAJ also in: https://groups.google.com/forum/#!msg/pyomo-forum/5b-UizHe0co/xYO4pK4iauMJ;context-place=searchin/pyomo-forum/extreme$20ray%7Csort:date but there did not seem to be a resolution here. I did not find much on stackexchange.

I also looked at the Benders decomposition example in the book from the ampl examples, but it seems this was only adding optimality cuts. But the suffixes of rc, urc, and lrc were unexplained. Is there any documentation or could there also be explanation of what these represent, or is there a more thorough list of all the suffix option for various solvers? I suspect the method to do what I wantinvolves suffixes but it is unclear.

[Bryan Arguello]

Eric, did you ever get an answer to this question?