# An application of the Xpress callbacks to a Benders decomposition algorithm.
# Uses preintsol callback for cut injection.
#
# Solves the problem:
# min c1*x + c2*y
# st A1*x + A2*y <=b (m constraints)
# x binary n1-dimensional vector
# y >=0 continuous n2-dimensional vector
#
# NOTE: This example is for illustration/tutorial purposes only. There is no
# benefit in solving the generic problem defined below using Benders
# decomposition.
# Courtesy of Georgios Patsakis (UC Berkeley, Amazon) and Richard L.-Y. Chen (Amazon).
#
# (C) 2019-2026 Fair Isaac Corporation
# Import necessary packages.
import xpress as xp
import sys
# Initialize problem parameters.
c1 = [1,6,5,7] # n1 x 1
c2 = [9,3,0,2,3] # n2 x 1
b = [-3,-4,1,4,5] # m x 1
A1 = [[0, -2, 3, 2],
[-5, 0, -3, 1],
[1, 0, 4, -2],
[0, -3, 4, -1],
[-5, -4, 3, 0]]
A2 = [[3, 4, 2, 0, -5],
[0, 2, 3, -2, 1],
[2, 0, 1, -3, -5],
[-5, 3, -2, -3, 0],
[-2, 3, -1, 2, -4]]
m = len(b)
n1 = len(c1)
n2 = len(c2)
# Absolute accuracy for optimal objective.
ObjAbsAccuracy=0.00001
# SOLVE ORIGINAL PROBLEM WITHOUT DECOMPOSING
print("**** Solving original problem without decomposing ****")
p = xp.problem()
# Define Variables.
x = [p.addVariable(vartype=xp.binary) for _ in range(n1)]
y = [p.addVariable() for _ in range(n2)] # positive real variable
# Define Constraints.
constr = [xp.Sum(A1[ii][jj]*x[jj] for jj in range(n1)) + \
xp.Sum(A2[ii][jj]*y[jj] for jj in range(n2)) \
<= b[ii] for ii in range(m)]
p.addConstraint(constr)
# Define Objective.
p.setObjective(xp.Sum(c1[jj]*x[jj] for jj in range(n1)) + \
xp.Sum(c2[jj]*y[jj] for jj in range(n2)) , \
sense=xp.minimize)
p.controls.outputlog = 0 # Suppress Xpress output
# Solve and retrieve solution.
p.optimize()
if p.attributes.solvestatus != xp.SolveStatus.COMPLETED or \
p.attributes.solstatus != xp.SolStatus.OPTIMAL:
raise RuntimeError('Problem could not be solved to MIP optimality')
xopt = p.getSolution(x)
yopt = p.getSolution(y)
xopt_rounded = [0.0 if round(val, 2) == 0 else round(val, 2) for val in xopt]
yopt_rounded = [0.0 if round(val, 2) == 0 else round(val, 2) for val in yopt]
print(f"Solution without decomposing: obj: {p.attributes.objval:.2f}")
print(f" x: {xopt_rounded}")
print(f" y: {yopt_rounded}")
print()
# SOLVE DECOMPOSED PROBLEM
# We define the following functions:
# - subproblem: function that solves the subproblem for a given first stage
# variable xhat. If the subproblem has an optimal solution, it returns the
# the optimal dual multiplier associated with the non anticipativity
# constraint, the optimal objective, and the label "Optimal". If the
# subproblem is infeasible, it returns the dual ray of unboundedness and
# the label "Infeasible".
#
# - preintsol_callback: A pre-integer solution callback is triggered every time
# an integer solution is found, BEFORE it is accepted by the solver. This
# callback validates the solution against the subproblem and can inject cuts
# directly when soltype=0 (integer solution from node relaxation).
# Subproblem Solver Function.
def subproblem(xhat):
# Input:
# - xhat: n1x1 array is the first stage variable,
# passed to the subproblem.
#
# Output:
# return (lamb: n1x1 array, beta: scalar, flag:string)
# - If the subproblem is Infeasible, the flag is 'Infeasible', and
# lamb (n1x1) and beta (1x1) are the components of the dual ray of
# unboundedness necessary to write a feasibility cut of the form:
# sum(x[i]*lamb[i] for all i) +beta >= 0
# - If the subproblem is Optimal, the flag is 'Optimal', beta (1x1) is the
# optimal objective, and lamb (n1x1) is the optimal dual multiplier
# associated with the non anticipativity constraint, so that the
# optimality cut will have the form:
# theta >= sum(x[i]*lamb[i] for all i) + beta
# IMPORTANT NOTE: The subproblem needs to ensure that upper bounds or non
# zero lower bounds should ONLY BE DEFINED THROUGH CONSTRAINTS, not as part
# of the variable definition. This is to ensure the feasibility cuts are
# correct.
# Side Note: If there are violated feasibility cuts, then the subproblem will be infeasible.
# In such a case, it is a good practice to ensure that the subproblem is always
# feasible through modeling, for instance by adding slacks to all the constraints.
# This often improves performance because it tends to separate stronger feasibility cuts.
# However, this is an advanced technique that is not implemented in this illustrative example.
# Define and Solve Subproblem.
# Initialize Problem.
r=xp.problem()
# Define Variables.
y=[r.addVariable() for _ in range(n2)] # Positive real variable - second stage.
z=[r.addVariable(lb=-xp.infinity) for _ in range(n1)] # Dummy copy variable,
# must have the exact shape as xhat.
epsilon=r.addVariable(lb=-xp.infinity) # Dummy variable to extract part of the
# infeasibility ray (if infeasible).
# Define Constraints.
# Dummy constraint for optimality/feasibility cuts.
# Note: make sure the index of the variable and of the position in the
# constraint array is the same (or has a known mapping).
dummy1= [z[i]==xhat[i] for i in range(n1)]
# Dummy constraint for feasibility cut.
dummy2= epsilon==1
# Second stage constraints, where RHS b has been multiplied by epsilon.
constr=[xp.Sum(A1[ii][jj]*z[jj] for jj in range(n1)) + \
xp.Sum(A2[ii][jj]*y[jj] for jj in range(n2)) \
- epsilon*b[ii]<=0 for ii in range(m)]
r.addConstraint(constr,dummy1,dummy2)
# Define Objective.
r.setObjective(xp.Sum(c2[jj]*y[jj] for jj in range(n2)),
sense=xp.minimize)
# Disable presolve. The only reason to do this is because in the case of
# an infeasible problem, the presolve might identify infeasibility of the
# problem without running simplex. Because of that, no dual infeasibility
# certificate will be found, hence a dual ray of unboundedness will not be
# returned, which means that a feasibility cut can not be formed.
r.controls.presolve = 0
# Silence output.
r.controls.outputlog = 0
# Solve optimization.
r.optimize()
# Find the indices of constraint dummy1, which will be used to access the
# dual multipliers corresponding to the entries of xhat.
xind1=[dummy1[ii].index for ii in range(n1)]
if r.attributes.solvestatus != xp.SolveStatus.COMPLETED:
print("ERROR: Subproblem was not solved. Terminating.")
sys.exit()
# Take cases depending on subproblem status.
if r.attributes.solstatus == xp.SolStatus.OPTIMAL:
# Optimal subproblem.
# Retrieve optimal dual multiplier and objective and return.
dualmult=r.getDuals()
lamb=[dualmult[ii] for ii in xind1]
beta=r.attributes.objval
return(lamb,beta,'Optimal')
elif r.attributes.solstatus == xp.SolStatus.INFEAS:
# Infeasible subproblem.
dray = r.getDualRay()
if dray is None:
print ("Could not retrieve a dual ray, return no good cut instead:")
# This is just a cheat if the subproblem is found infeasible and
# the optimizer fails to find a dual ray. Since all the first stage
# variables are binary, we add a No-Good-Cut to exclude the point
# xhat instead of an infeasibility cut.
lamb=[2*xhat[ii]-1 for ii in range(n1)]
beta=-sum(xhat)+1
else:
# Extract the dual ray.
print ("Dual Ray:", dray)
# Extract the part of the dual ray corresponding to the first stage
# variables xhat.
lamb=[dray[ii] for ii in xind1]
# Extract the constant from the dual ray entry of constraint dummy 2.
beta=dray[dummy2.index]
return(lamb,beta,'Infeasible')
else:
print("ERROR: Subproblem not optimal or infeasible. Terminating.")
sys.exit()
# Pre-Integer Solution Callback Function
def preintsol_callback(p, data, soltype, cutoff):
# Input:
# NOTE: the input is populated by the solver when an integer solution is
# found, BEFORE it is accepted.
# - p: an xp.problem() (a thread-local copy of the problem from which the callback was triggered)
# - data: structure that is passed to the callback from the problem,
# specified in the addPreIntsolCallback function. No data is needed here,
# so the addPreIntsolCallback() call uses None.
# - soltype: 0 if the integer solution was found by a node relaxation,
# 1 if found by a heuristic, 2 if provided by user
# - cutoff: the cutoff value that the Optimizer will use if the solution
# is accepted.
#
# Output:
# return (reject, newcutoff)
# - reject: 1 to reject the solution, 0 to accept it
# - newcutoff: new cutoff value (or None to keep the current one)
# Retrieve xhat and thetahat from the candidate integer solution.
xhat = p.getCallbackSolution([x[ii] for ii in range(n1)])
thetahat = p.getCallbackSolution(theta)
xhat_rounded = [0.0 if round(val, 2) == 0 else round(val, 2) for val in xhat]
print(f"Integer solution found! x={xhat_rounded}, theta={thetahat:.2f}")
# Solve the subproblem to check feasibility.
(dual_mult, opt, status) = subproblem(xhat)
print(f" Solving subproblem... subproblem obj = {opt:.2f}")
if status == 'Infeasible':
# The subproblem is infeasible. Add a feasibility cut.
# For heuristic solutions (soltype != 0), we cannot add cuts,
# and we can only reject the solution if it is not feasible.
# Depending on the specific application, the heuristics could be
# very unlikely to find valid feasible solutions (they don't
# respect the Benders cuts we're generating), and the Benders
# subproblem could be expensive to solve.
# In such cases, it could be more convenient to reject all
# heuristic solutions without solving the subproblem, or to
# simply disable heuristics completely.
if soltype != 0:
print(" -> Subproblem infeasible - Heuristic solution, rejecting")
return (1, None)
print(" -> Subproblem infeasible - Adding feasibility cut")
xind = [x[ii].index for ii in range(n1)]
coefficients = dual_mult
rhs = -opt
p.addCuts([1], ['L'], [rhs], [0, len(xhat)], xind, coefficients)
# IMPORTANT: Return 0 (accept) to keep the cuts. Returning reject=1 would drop the cuts.
# When soltype=0, addCuts() causes the node LP to be resolved with the new cuts.
return (0, None)
elif status == 'Optimal':
# Check if theta is sufficient for the optimal subproblem value.
if thetahat >= opt - ObjAbsAccuracy:
# The solution is feasible to the full problem.
print(f" -> Subproblem obj <= theta -> Valid solution!")
return (0, None)
else:
# The solution is infeasible; theta is too small. Add an optimality cut.
# For heuristic solutions (soltype != 0), we cannot add cuts.
if soltype != 0:
print(" -> Subproblem obj > theta - Heuristic solution, rejecting")
return (1, None)
print(f" -> Subproblem obj > theta -> Adding optimality cut")
xind = [x[ii].index for ii in range(n1)]
thetaind = [theta.index]
coefficients = [1] + [-dual_mult[ii] for ii in range(len(dual_mult))]
rhs = opt - sum(dual_mult[ii]*xhat[ii] for ii in range(len(dual_mult)))
p.addCuts([1], ['G'], [rhs], [0, len(xhat)+1],
thetaind+xind, coefficients)
# IMPORTANT: Return 0 (accept) to keep the cuts. Returning reject=1 would drop the cuts.
# When soltype=0, addCuts() causes the node LP to be resolved with the new cuts.
return (0, None)
else:
print("ERROR: Unexpected subproblem status. Rejecting solution.")
return (1, None)
# Main Problem.
print("**** Benders decomposition algorithm ****")
# Define main problem.
p=xp.problem()
# Define first stage variables.
x=[p.addVariable(vartype=xp.binary) for _ in range(n1)]
theta=p.addVariable() # (positive) second stage cost - change the default lower
# bound if second stage cost may not be positive.
# Define objective.
p.setObjective(xp.Sum(c1[jj]*x[jj] for jj in range(n1)) + theta,
sense=xp.minimize)
# Add Pre-Integer Solution callback.
# This callback is called when an integer solution is found, BEFORE it is
# accepted by the solver. When soltype=0 (node relaxation solution), we can
# add cuts directly from this callback, which will cause the node LP to be
# resolved.
p.addPreIntsolCallback(preintsol_callback, None, 0)
# Deactivate presolve and symmetry detection. This is to ensure that the solver will not eliminate
# variables, so the indices we retrieve for the variables and the cuts we add
# correspond to/contain existing optimization variables. Presolve can be
# activated, but additional steps are required for that (presolving the cuts
# the same way as the rows of the matrix).
p.setControl({"presolve":0,"mippresolve":0,"symmetry":0})
p.controls.outputlog = 0 # Suppress Xpress output
# Solve the main problem.
p.optimize()
print()
if p.attributes.solvestatus == xp.SolveStatus.COMPLETED and \
p.attributes.solstatus == xp.SolStatus.OPTIMAL:
# Print results.
print("**** Final solution (Benders) ****")
xsol = p.getSolution(x)
# Solve subproblem one final time to get y values
final_sub = xp.problem()
final_sub.controls.outputlog = 0
y_final = [final_sub.addVariable() for _ in range(n2)]
constr_final = [xp.Sum(A1[ii][jj]*xsol[jj] for jj in range(n1)) + \
xp.Sum(A2[ii][jj]*y_final[jj] for jj in range(n2)) <= b[ii]
for ii in range(m)]
final_sub.addConstraint(constr_final)
final_sub.setObjective(xp.Sum(c2[jj]*y_final[jj] for jj in range(n2)), sense=xp.minimize)
final_sub.optimize()
ysol = final_sub.getSolution(y_final)
xsol_rounded = [0.0 if round(val, 2) == 0 else round(val, 2) for val in xsol]
ysol_rounded = [0.0 if round(val, 2) == 0 else round(val, 2) for val in ysol]
print(f"Solution with Benders decomposition: obj: {p.attributes.objval:.2f}")
print(f" x: {xsol_rounded}")
print(f" y: {ysol_rounded}")
else:
print("Could not solve the problem")
|
