Description: This notebook describes how to implement Benders decomposition, which is a large scale optimization scheme, in JuMP. Both the classical approach (using loop) and the modern approach (using lazy constraints) are described.
Author: Shuvomoy Das Gupta
License:
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
To illustrate implementation of solver callback in JuMP, we consider applying Benders decomposition to the following general mixed integer problem:
\begin{align} & \text{maximize} \quad &&c_1^T x+c_2^T v \\ & \text{subject to} \quad &&A_1 x+ A_2 v \preceq b \\ & &&x \succeq 0, x \in \mathbb{Z}^n \\ & &&v \succeq 0, v \in \mathbb{R}^p \\ \end{align}where $b \in \mathbb{R}^m$, $A_1 \in \mathbb{R}^{m \times n}$ and $A_2 \in \mathbb{R}^{m \times p}$. Here the symbol $\succeq$ ($\preceq$) stands for element-wise greater (less) than or equal to. Any mixed integer programming problem can be written in the form above.
We want to write the Benders decomposition algorithm for the problem above. Consider the polyhedron $\{u \in \mathbb{R}^m| A_2^T u \succeq 0, u \succeq 0\}$. Assume the set of vertices and extreme rays of the polyhedron is denoted by $P$ and $Q$ respectively. Assume on the $k$th iteration the subset of vertices of the polyhedron mentioned is denoted by $T(k)$ and the subset of extreme rays are denoted by $Q(k)$, which will be generated by the Benders decomposition problem below.
Step 1 (Initialization)
We start with $T(1)=Q(1)=\emptyset$. Let $f_m^{(1)}$ be arbitrarily large and $x^{(1)}$ be any non-negative integer vector and go to Step 3.
** Step 2 (Solving the master problem)**
Solve the master problem:
Let the maximizer corresponding to the objective value $f_\text{m}^{(k)}$ be denoted by $x^{(k)}$. Now there are three possibilities:
If $f_\text{m}^{(k)}=-\infty$, i.e., the master problem is infeasible, then the original proble is infeasible and sadly, we are done.
If $f_\text{m}^{(k)}=\infty$, i.e. the master problem is unbounded above, then we take $f_\text{m}^{(k)}$ to be arbitrarily large and $x^{(k)}$ to be a corresponding feasible solution. Go to Step 3
If $f_\text{m}^{(k)}$ is finite, then we collect $t^{(k)}$ and $x^{(k)}$ and go to Step 3.
** Step 3 (Solving the subproblem and add Benders cut when needed) **
Solve the subproblem
Let the minimizer corresponding to the objective value $f_s(x^{(k)})$ be denoted by $u^{(k)}$. There are three possibilities:
If $f_s(x^{(k)}) = \infty$, the original problem is either infeasible or unbounded. We quit from Benders algorithm and use special purpose algorithm to find a feasible solution if there exists one.
If $f_s(x^{(k)}) = - \infty$, we arrive at an extreme ray $y^{(k)}$. We add the Benders cut corresponding to this extreme ray $(A_1 ^T y^{(k)})^T x \leq b^T y^{(k)}$ to the master problem i.e., $Q(k+1):= Q(k) \cup \{y^{(k)}\}$. Take $k:=k+1$ and go to Step 3.
If $f_s(x^{(k)})$ is finite, then
If $f_s(x^{(k)})=f_m^{(k)}$ we arrive at the optimal solution. The optimum objective value of the original problem is $f_s(x^{(k)})=f_m^{(k)}$, an optimal $x$ is $x^{(k)}$ and an optimal $v$ is the dual values for the second constraints of the subproblem. We are happily done!
If $f_s(x^{(k)}) < f_m^{(k)}$ we get an suboptimal vertex $u^{(k)}$. We add the corresponding Benders cut $u_0 + (A_1^T u^{(k)} - c_1)^T x \leq b^T u^{(k)}$ to the master problem, i.e., $T(k+1) := T(k) \cup \{u^{(k)}\}$. Take $k:=k+1$ and go to Step 3.
The input data is from page 139, Integer programming by Garfinkel and Nemhauser, 1972.
# Data for the problem
#---------------------
c1=[-1;-4]
c2=[-2; -3]
dimX=length(c1)
dimU=length(c2)
b=[-2;-3]
A1=[
1 -3;
-1 -3
]
A2=[
1 -2;
-1 -1
]
M=1000
1000
There are two ways we can implement Benders decomposition in JuMP:
The classical approach might be inferior to the modern one, as the solver
For more details on the comparison between the two approaches, see Paul Rubin's blog on Benders Decomposition.
Let's describe the master problem first. Note that there are no constraints, which we will added later using Benders decomposition.
# Loading the necessary packages
#-------------------------------
using JuMP
using GLPK
using MathOptInterface
const MOI = MathOptInterface
# Master Problem Description
# --------------------------
masterProblemModel = Model(optimizer=GLPK.GLPKOptimizerMIP())
# Variable Definition
# ----------------------------------------------------------------
@variable(masterProblemModel, 0<= x[1:dimX]<= 1e6 , Int)
@variable(masterProblemModel, t<=1e6)
# Objective Setting
# -----------------
@objective(masterProblemModel, Max, t)
iC=1 # iC is the iteration counter
print(masterProblemModel)
A JuMP Model
Here is the loop that checks the status of the master problem and the subproblem and then adds necessary Benders cuts accordingly.
# Trying the entire benders decomposition in one step
while(true)
println("\n-----------------------")
println("Iteration number = ", iC)
println("-----------------------\n")
println("The current master problem is")
print(masterProblemModel)
JuMP.optimize(masterProblemModel)
t_status = JuMP.terminationstatus(masterProblemModel)# == MOI.Success
p_status = JuMP.primalstatus(masterProblemModel)# == MOI.FeasiblePoint
if p_status == MOI.InfeasiblePoint
println("The problem is infeasible :-(")
break
end
if t_status == MOI.InfeasibleOrUnbounded
fmCurrent = M
xCurrent=M*ones(dimX)
end
if p_status == MOI.FeasiblePoint
fmCurrent = JuMP.resultvalue(t)
xCurrent=Float64[]
for i in 1:dimX
push!(xCurrent, JuMP.resultvalue(x[i]))
end
end
println("Status of the master problem is", t_status,
"\nwith fmCurrent = ", fmCurrent,
"\nxCurrent = ", xCurrent)
subProblemModel = Model(optimizer=GLPK.GLPKOptimizerLP())
cSub=b-A1*xCurrent
@variable(subProblemModel, u[1:dimU]>=0)
@constraint(subProblemModel, constrRefSubProblem[j=1:size(A2,2)],sum(A2[i,j]*u[i] for i in 1:size(A2,1))>=c2[j])
# The second argument of @constraint macro, constrRefSubProblem[j=1:size(A2,2)] means that the j-th constraint is
# referenced by constrRefSubProblem[j].
@objective(subProblemModel, Min, dot(c1, xCurrent) + sum(cSub[i]*u[i] for i in 1:dimU))
print("\nThe current subproblem model is \n", subProblemModel)
JuMP.optimize(subProblemModel)
t_status_sub = JuMP.terminationstatus(subProblemModel)# == MOI.Success
p_status_sub = JuMP.primalstatus(subProblemModel)# == MOI.FeasiblePoint
fsxCurrent = JuMP.objectivevalue(subProblemModel)
uCurrent = Float64[]
for i in 1:dimU
push!(uCurrent, JuMP.resultvalue(u[i]))
end
γ=dot(b,uCurrent)
println("Status of the subproblem is ", t_status_sub,
"\nwith fsxCurrent= ", fsxCurrent,
"\nand fmCurrent= ", fmCurrent)
if p_status_sub == MOI.FeasiblePoint && fsxCurrent == fmCurrent # we are done
println("\n################################################")
println("Optimal solution of the original problem found")
println("The optimal objective value t is ", fmCurrent)
println("The optimal x is ", xCurrent)
println("The optimal v is ", JuMP.resultdual.(constrRefSubProblem))
println("################################################\n")
break
end
if p_status_sub == MOI.FeasiblePoint && fsxCurrent < fmCurrent
println("\nThere is a suboptimal vertex, add the corresponding constraint")
cv= A1'*uCurrent - c1
@constraint(masterProblemModel, t+sum(cv[i]*x[i] for i in 1:dimX) <= γ )
println("t + ", cv, "ᵀ x <= ", γ)
end
if t_status_sub == MOI.InfeasibleOrUnbounded
println("\nThere is an extreme ray, adding the corresponding constraint")
ce = A1'*uCurrent
@constraint(masterProblemModel, sum(ce[i]*x[i] for i in 1:dimX) <= γ)
println(ce, "ᵀ x <= ", γ)
end
iC=iC+1
sleep(0.5)
end
----------------------- Iteration number = 1 ----------------------- The current master problem is A JuMP ModelStatus of the master problem isSuccess with fmCurrent = 1.0e6 xCurrent = [0.0, 0.0] The current subproblem model is A JuMP ModelStatus of the subproblem is Success with fsxCurrent= -7.666666666666667 and fmCurrent= 1.0e6 There is a suboptimal vertex, add the corresponding constraint t + [-1.0, -4.0]ᵀ x <= -7.666666666666667 ----------------------- Iteration number = 2 ----------------------- The current master problem is A JuMP ModelStatus of the master problem isSuccess with fmCurrent = 1.0e6 xCurrent = [0.0, 250002.0] The current subproblem model is A JuMP ModelStatus of the subproblem is Success with fsxCurrent= -1.000008e6 and fmCurrent= 1.0e6 There is a suboptimal vertex, add the corresponding constraint t + [1.0, 4.0]ᵀ x <= 0.0 ----------------------- Iteration number = 3 ----------------------- The current master problem is A JuMP ModelStatus of the master problem isSuccess with fmCurrent = -4.0 xCurrent = [4.0, 0.0] The current subproblem model is A JuMP ModelStatus of the subproblem is Success with fsxCurrent= -13.0 and fmCurrent= -4.0 There is a suboptimal vertex, add the corresponding constraint t + [2.5, -0.5]ᵀ x <= -3.0 ----------------------- Iteration number = 4 ----------------------- The current master problem is A JuMP ModelStatus of the master problem isSuccess with fmCurrent = -4.0 xCurrent = [0.0, 1.0] The current subproblem model is A JuMP ModelStatus of the subproblem is Success with fsxCurrent= -4.0 and fmCurrent= -4.0 ################################################ Optimal solution of the original problem found The optimal objective value t is -4.0 The optimal x is [0.0, 1.0] The optimal v is [0.0, 0.0] ################################################
Mixed integer programming solver works on branch-and-bound strategy. Often in a mixed integer programming problem, including all possible constraints might be too space consuming or computationally too expensive. Instead we can start with a manageable set of constraints in a comparatively smaller, hence relaxed version of the original MIP. Once a feasible integer solution is found, we can check the status of the problem by solving some subproblem. The subproblem can be derived from duality, as in our current problem and/or from logic. In the case of suboptimality, we can add lazy constraints to cut off the current feasible solution at that node of the branch-and-bound tree.
Lazy constraints have the following advantages:
Note that, in Step 3 we are solving the subproblem, checking the state of our problem and adding Benders cut if necessary. We write a function addBendersCut(cb)
which will perform the steps. The argument of the function cb
refers to the callback handle, which is a reference to the callback management code inside JuMP.
When we add the lazy constraints we have to use the @addLazyConstraint
macro as follows:
``@addLazyConstraint(cb, description of the constraint)
``.
Description of the constraint will be of the same manner as in adding a normal constraint in JuMP using @constraint
macro.
Note that we have not mentioned which model is associated with the added lazy constraints. So outside the addBendersCut(cb)
function we invoke the command:
``addLazyCallback(name of the master problem model, addBendersCut)
``
, which tells JuMP to add the lazy constraints generated by the function addBendersCut
to the master problem.
# Loading the necessary packages
#-------------------------------
using JuMP
#using CPLEX
using Gurobi
# Master Problem Description
# --------------------------
# Model name
#masterProblemModel = Model(solver = CplexSolver())
masterProblemModel = Model(optimizer = GurobiOptimizer(Heuristics=0, Cuts = 0, OutputFlag = 0)) # If we want to add Benders lazy constraints
# in Gurobi, then we have to turn off Gurobi's own Cuts and Heuristics in the master problem
# Variable Definition (Only CplexSolver() works properly for these)
# ----------------------------------------------------------------
#@variable(masterProblemModel, x[1:dimX] >=0 , Int)
#@variable(masterProblemModel, t)
# ***************ALTERNATIVE VARIABLE DEFINITION FOR GUROBI************
#If we replace the two lines above with the follwoing:
@variable(masterProblemModel, 0<= x[1:dimX] <= 1e6 , Int)
@variable(masterProblemModel, t <= 1e6)
# then all the solvers give the expected solution
#**********************************************************************
# Objective Setting
# -----------------
@objective(masterProblemModel, Max, t)
print(masterProblemModel)
stringOfBenderCuts=AbstractString[] # this is an array of strings which will contain all the
# Benders cuts added to be displayed later
# There are no constraints when we start, so we will add all the constraints in the
# form of Benders cuts lazily
function addBendersCut(cb)
#***************************************************************************
# First we store the master problem solution in conventional data structures
println("----------------------------")
println("ITERATION NUMBER = ", length(stringOfBenderCuts)+1)
println("---------------------------\n")
fmCurrent = getvalue(t)
xCurrent=Float64[]
for i in 1:dimX
push!(xCurrent, JuMP.resultvalue(x[i]))
end
# Display the current solution of the master problem
println("MASTERPROBLEM INFORMATION")
println("-------------------------")
println("The master problem that was solved was:")
print(masterProblemModel)
println("with ", length(stringOfBenderCuts), " added lazy constraints")
println(stringOfBenderCuts)
println("Current Value of x is: ", xCurrent)
println("Current objective value of master problem, fmCurrent is: ", fmCurrent)
println("\n")
#************************************************************************
# ========================================================================
# Now we solve the subproblem
# subProblemModel=Model(solver=CplexSolver())
subProblemModel = Model(optimizer=GurobiOptimizer(Presolve=0, OutputFlag = 0))
cSub=b-A1*xCurrent
@variable(subProblemModel, u[1:dimU]>=0)
@constraint(subProblemModel, constrRefSubProblem[j=1:size(A2,2)], sum(A2[i,j]*u[i] for i in 1:size(A2,1))>=c2[j])
@objective(subProblemModel, Min, dot(c1, xCurrent) + sum(cSub[i]*u[i] for i in 1:dimU))
println("The subproblem is being solved")
statusSubProblem = JuMP.optimize(subProblemModel)
# We store the results achieved from the subproblem in conventional data structures
fsxCurrent = JuMP.objectivevalue(subProblemModel)
uCurrent = Float64[]
for i in 1:dimU
push!(uCurrent, JuMP.resultvalue(u[i]))
end
# Display the solution corresponding to the subproblem
println("SUBPROBLEM INFORMATION")
println("----------------------")
println("The subproblem that was solved was: ")
print(subProblemModel)
println("Current status of the subproblem is ", statusSubProblem)
println("Current Value of u is: ", uCurrent) # JuMP will return an extreme ray
# automatically (if the solver supports it), so we do not need to change the syntax
println("Current Value of fs(xCurrent) is: ", fsxCurrent)
println("\n")
# ==========================================================================
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Now we check the status of the algorithm and add Benders cut when necessary
γ=dot(b,uCurrent)
if statusSubProblem == :Optimal && fsxCurrent==fmCurrent # we are done
println("OPTIMAL SOLUTION OF THE ORIGINAL PROBLEM FOUND :-)")
println("The optimal objective value t is ", fmCurrent)
println("The optimal x is ", xCurrent)
println("The optimal v is ", JuMP.resultdual(constrRefSubProblem))
println("\n")
return
end
println("-------------------ADDING LAZY CONSTRAINT----------------")
if statusSubProblem == :Optimal && fsxCurrent < fmCurrent
println("\nThere is a suboptimal vertex, add the corresponding constraint")
cv= A1'*uCurrent - c1
@lazyconstraint(cb, t+sum(cv[i]*x[i] for i in 1:dimX) <= γ)
println("t + ", cv, "ᵀ x <= ", γ)
push!(stringOfBenderCuts, string("t+", cv, "'x <=", γ))
end
if statusSubProblem == :Unbounded
println("\nThere is an extreme ray, adding the corresponding constraint")
ce = A1'*uCurrent
@lazyconstraint(cb, sum(ce[i]*x[i] for i in 1:dimX) <= γ)
println(ce, "x <= ", γ)
push!(stringOfBenderCuts, string(ce, "ᵀ x <= ", γ))
end
println("\n")
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
end
Now we tell the solver to use the callback function and solve the problem.
addlazycallback(masterProblemModel, addBendersCut) # Telling the solver to use the
# callback function
JuMP.optimize(masterProblemModel)