Disclaimer

This notebook is only working under the versions:

  • JuMP 0.19 (unreleased, but currently in master)

  • MathOptInterface 0.4.1

  • GLPK 0.6.0

Disclaimer 2

The second part od this notebook (using Lazy Constraints) is not working

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: Creative Commons License
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.

Using Julia+JuMP for optimization - Benders decomposition



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.

Benders decomposition algorithm


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:

\begin{align} & f_\text{m}^{(k)}= \\ &\text{maximize} &&\quad t \\ &\text{subject to} \quad &&\forall \bar{u} \in T(k) \qquad t + (A_1^T \bar{u} - c_1)^T x \leq b^T \bar{u} \\ & && \forall \bar{y} \in Q(k) \qquad (A_1 ^T \bar{y})^T x \leq b^T \bar{y} \\ & && \qquad \qquad \qquad \; x \succeq 0, x \in \mathbb{Z}^n \end{align}

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

\begin{align} f_s(x^{(k)})= \\ c_1^T x^{(k)} + & \text{minimize} && (b-A_1 x^{(k)})^T u \\ & \text{subject to} && A_2^T u \succeq c_2 \\ & && u \succeq 0, u \in \mathbb{R}^m \end{align}

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.

Data for the problem

The input data is from page 139, Integer programming by Garfinkel and Nemhauser, 1972.

In [15]:
# 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
Out[15]:
1000

How to implement the Benders decomposition algorithm in JuMP

There are two ways we can implement Benders decomposition in JuMP:

  • Classical approach: Adding the Benders cuts in a loop, and
  • Modern approach: Adding the Benders cuts as lazy constraints.

The classical approach might be inferior to the modern one, as the solver

  • might revisit previously eliminated solution, and
  • might discard the optimal solution to the original problem in favor of a better but ultimately infeasible solution to the relaxed one.

For more details on the comparison between the two approaches, see Paul Rubin's blog on Benders Decomposition.

Classical approach: adding the Benders cuts in a loop

Let's describe the master problem first. Note that there are no constraints, which we will added later using Benders decomposition.

In [16]:
# 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.

In [17]:
# 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]
################################################

In [ ]:

In [ ]:

In [ ]:

In [ ]:

In [ ]:

In [ ]:

In [ ]:

In [ ]:

In [ ]:

In [ ]:

In [ ]:

In [ ]:

In [ ]:

In [ ]:

In [ ]:

In [ ]:

In [ ]:

Modern approach: adding the Benders cuts as lazy constraints

What are lazy constraints?

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:

  • does not revisit previously eliminated solution, and
  • does not discard the optimal solution to the original problem in favor of a better but ultimately infeasible solution to the relaxed one.

How to add lazy constraints in Julia?

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.

In [ ]:
# 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
In [ ]:
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.

In [ ]:
addlazycallback(masterProblemModel, addBendersCut) # Telling the solver to use the 
# callback function
In [ ]:
JuMP.optimize(masterProblemModel)