Jia Guo
Jia Guo

Reputation: 1

When calling Gurobi by pyomo, if the primal is infeasible, how to find a random dual extreme points and extreme ray without coding the dual?

I am using pyomo to call gurobi to solve a linear programming model P. Let D be the dual problem of P. When P is infeasible, how can we get a random (not fixed) extreme point and extreme ray of D without coding D in pyomo?

For simplicity, I used a tiny example with only two non-negative variables (x1 and x2). The objective is to minimize 30x1 - 20x2, subject to -6x1 + 4x2 >= 12 and 5x1 - 5x2 >= 5.

Primal P is infeasible. By geometry, we know that the dual extreme points are (0, 6) and (0, 4). Dual extreme rays are (1, 0.8) and (1, 1.2). However, from the code below, I always get dual extreme point (0, 4) and dual extreme ray (-1, -0.8). I cannot get dual extreme point (0, 6) and dual extreme ray (1, 1.2). It seems that we can only get one fixed dual extreme point and one fixed dual extreme ray.

The code is the following:

from pyomo.environ import *
import pyomo.environ as pe
import pyomo.opt
import gurobipy as gp

# Define model
model = ConcreteModel()
model.x = Var([1, 2], within=NonNegativeReals)

# Objective function
model.obj = Objective(expr=30*model.x[1] - 20 * model.x[2], sense=minimize)

# Infeasible constraints
model.a_rule = Constraint(expr=-6 * model.x[1] + 4 * model.x[2] >= 12)
model.b_rule = Constraint(expr= 5 * model.x[1] - 5 * model.x[2] >= 5)

    # Solve with Gurobi
solver = SolverFactory('gurobi')
solver.options['Method'] = 1  # Use dual simplex
solver.options['InfUnbdInfo'] = 1  # Ensure infeasibility certificate is available
results = solver.solve(model, tee=True, load_solutions=False)  # Do NOT load solutions

# Find dual extreme points and extreme rays when the model is infeasible
if results.solver.termination_condition == TerminationCondition.infeasible:
     lp_filename = "pyomo_model.lp"
     gurobi_model = gp.read(lp_filename)
     gurobi_model.setParam("Method", 1)  # Ensure Dual Simplex is used
     gurobi_model.setParam("InfUnbdInfo", 1)  # Enable infeasibility certificate
     gurobi_model.setParam("DualReductions", 0)  # Ensure extreme ray is available
     gurobi_model.setParam("PreSolve", 0)   # Disable preprocessing
     gurobi_model.optimize()

 # Find dual extreme points
 print ('Dual extreme points:')
 dual_values = gurobi_model.getAttr("Pi", gurobi_model.getConstrs())
 for constr, dual_value in zip(gurobi_model.getConstrs(), dual_values):
    print (constr, dual_value)
     
 # Find dual extreme rays
 print ('Dual extreme rays:')
 for constr in gurobi_model.getConstrs():
      farkas_dual = constr.getAttr("FarkasDual")
      print (constr, farkas_dual)

Upvotes: 0

Views: 12

Answers (0)

Related Questions