Benders decomposition (via callbacks)

Originally Contributed by: Mathieu Besançon

This notebook describes how to implement the Benders decomposition in JuMP using lazy constraints. We keep the same notation and problem form as the first notebook Benders decomposition.

\[\begin{aligned} & \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{aligned}\]

where $b \in \mathbb{R}^m$, $A_1 \in \mathbb{R}^{m \times n}$, $A_2 \in \mathbb{R}^{m \times p}$ and $\mathbb{Z}$ is the set of integers.

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.

For a detailed explanation on the Benders decomposition algorithm, see the introduction notebook.

Lazy constraints

Some optimization solvers allow users to interact with them during the solution process by providing user-defined functions which are triggered under certain conditions.

The generic term for these functions is callback.

In integer optimization, the main callback types are lazy constraints, user-cuts and heuristic solutions. See the Callbacks section for an introduction on how to use them.

Some callbacks define a new constraint which is only activated when necessary, i.e., when a current solution does not respect them. It can avoid building an optimization model with too many constraints up-front.

This is the case for Benders decomposition, since the sub-problem defines an exponential number of primal vertices and therefore dual cuts.

A detailed explanation on the distinction between user-cuts and lazy constraints is also available on Paul Rubin's blog. He also describes this approach to Benders Decomposition.

We use the data from the original notebook and change the solution algorithm to leverage lazy constraints:

  • Step 1 (Initialization)
  • Step 2 (defining the subproblem model)
  • Step 3 (registering the lazy constraint of the subproblem)

Input data

The input data is from page 139, Integer programming by Garfinkel and Nemhauser[1].

c1 = [-1; -4]
c2 = [-2; -3]

dim_x = length(c1)
dim_u = length(c2)

b = [-2; -3]

A1 = [
    1 -3
    -1 -3
]
A2 = [
    1 -2
    -1 -1
]

M = 1000;

Implementation

using JuMP
import GLPK
import Test
Info

This tutorial uses the MathOptInterface API. By default, JuMP exports the MOI symbol as an alias for the MathOptInterface.jl package. We recommend making this more explicit in your code by adding the following lines:

import MathOptInterface
const MOI = MathOptInterface

Subproblem creation

function build_subproblem()
    sub_problem_model = Model(GLPK.Optimizer)
    @variable(sub_problem_model, u[1:dim_u] >= 0)
    @constraint(
        sub_problem_model,
        constr_ref_subproblem[j = 1:size(A2, 2)],
        A2[:, j]' * u >= c2[j],
    )
    return (sub_problem_model, u)
end
build_subproblem (generic function with 1 method)

Master Problem Description

master_problem_model = Model(GLPK.Optimizer)
@variable(master_problem_model, 0 <= x[1:dim_x] <= 1e6, Int)
@variable(master_problem_model, t <= 1e6)
@objective(master_problem_model, Max, t)
print(master_problem_model)
Max t
Subject to
 x[1] ≥ 0.0
 x[2] ≥ 0.0
 x[1] ≤ 1.0e6
 x[2] ≤ 1.0e6
 t ≤ 1.0e6
 x[1] integer
 x[2] integer
(sub_problem_model, u) = build_subproblem();

Track the calls to the callback

iter_num = 0
0

Define lazy constraints

function benders_lazy_constraint_callback(cb_data)
    global iter_num
    iter_num += 1
    println("Iteration number = ", iter_num)

    x_current = callback_value.(Ref(cb_data), x)
    fm_current = callback_value(cb_data, t)

    c_sub = b - A1 * x_current
    @objective(sub_problem_model, Min, c1' * x_current + c_sub' * u)
    optimize!(sub_problem_model)

    print("\nThe current subproblem model is \n")
    print(sub_problem_model)

    t_status_sub = termination_status(sub_problem_model)
    p_status_sub = primal_status(sub_problem_model)

    fs_x_current = objective_value(sub_problem_model)

    u_current = value.(u)

    γ = b' * u_current

    if p_status_sub == FEASIBLE_POINT && fs_x_current ≈ fm_current # we are done
        @info("No additional constraint from the subproblem")
    end

    if p_status_sub == FEASIBLE_POINT && fs_x_current < fm_current
        println(
            "\nThere is a suboptimal vertex, add the corresponding constraint",
        )
        cv = A1' * u_current - c1
        new_optimality_cons = @build_constraint(t + cv' * x <= γ)
        MOI.submit(
            master_problem_model,
            MOI.LazyConstraint(cb_data),
            new_optimality_cons,
        )
    end

    if t_status_sub == INFEASIBLE_OR_UNBOUNDED
        println(
            "\nThere is an  extreme ray, adding the corresponding constraint",
        )
        ce = A1' * u_current
        new_feasibility_cons = @build_constraint(dot(ce, x) <= γ)
        MOI.submit(
            master_problem_model,
            MOI.LazyConstraint(cb_data),
            new_feasibility_cons,
        )
    end
end

MOI.set(
    master_problem_model,
    MOI.LazyConstraintCallback(),
    benders_lazy_constraint_callback,
)

optimize!(master_problem_model)

t_status = termination_status(master_problem_model)
p_status = primal_status(master_problem_model)

if p_status == INFEASIBLE_POINT
    println("The problem is infeasible :-(")
end

if t_status == INFEASIBLE_OR_UNBOUNDED
    fm_current = M
    x_current = M * ones(dim_x)
end

if p_status == FEASIBLE_POINT
    fm_current = value(t)
    x_current = value.(x)
end


println(
    "Status of the master problem is ",
    t_status,
    "\nwith fm_current = ",
    fm_current,
    "\nx_current = ",
    x_current,
)
Iteration number = 1

The current subproblem model is
Min -2 u[1] - 3 u[2]
Subject to
 constr_ref_subproblem[1] : u[1] - u[2] ≥ -2.0
 constr_ref_subproblem[2] : -2 u[1] - u[2] ≥ -3.0
 u[1] ≥ 0.0
 u[2] ≥ 0.0

There is a suboptimal vertex, add the corresponding constraint
Iteration number = 2

The current subproblem model is
Min 750003.75 u[1] + 750002.75 u[2] - 1.0000076666666666e6
Subject to
 constr_ref_subproblem[1] : u[1] - u[2] ≥ -2.0
 constr_ref_subproblem[2] : -2 u[1] - u[2] ≥ -3.0
 u[1] ≥ 0.0
 u[2] ≥ 0.0

There is a suboptimal vertex, add the corresponding constraint
Iteration number = 3

The current subproblem model is
Min 0.875 u[1] - 0.125 u[2] - 3.8333333333333335
Subject to
 constr_ref_subproblem[1] : u[1] - u[2] ≥ -2.0
 constr_ref_subproblem[2] : -2 u[1] - u[2] ≥ -3.0
 u[1] ≥ 0.0
 u[2] ≥ 0.0

There is a suboptimal vertex, add the corresponding constraint
Iteration number = 4

The current subproblem model is
Min -3.8333333333333335
Subject to
 constr_ref_subproblem[1] : u[1] - u[2] ≥ -2.0
 constr_ref_subproblem[2] : -2 u[1] - u[2] ≥ -3.0
 u[1] ≥ 0.0
 u[2] ≥ 0.0
[ Info: No additional constraint from the subproblem
Iteration number = 5

The current subproblem model is
Min -0.875 u[1] + 0.125 u[2] - 3.8333333333333335
Subject to
 constr_ref_subproblem[1] : u[1] - u[2] ≥ -2.0
 constr_ref_subproblem[2] : -2 u[1] - u[2] ≥ -3.0
 u[1] ≥ 0.0
 u[2] ≥ 0.0

There is a suboptimal vertex, add the corresponding constraint
Iteration number = 6

The current subproblem model is
Min u[2] - 5
Subject to
 constr_ref_subproblem[1] : u[1] - u[2] ≥ -2.0
 constr_ref_subproblem[2] : -2 u[1] - u[2] ≥ -3.0
 u[1] ≥ 0.0
 u[2] ≥ 0.0
[ Info: No additional constraint from the subproblem
Iteration number = 7

The current subproblem model is
Min u[1] - 4
Subject to
 constr_ref_subproblem[1] : u[1] - u[2] ≥ -2.0
 constr_ref_subproblem[2] : -2 u[1] - u[2] ≥ -3.0
 u[1] ≥ 0.0
 u[2] ≥ 0.0
[ Info: No additional constraint from the subproblem
Status of the master problem is OPTIMAL
with fm_current = -4.0
x_current = [0.0, 1.0]

References

  1. Garfinkel, R. & Nemhauser, G. L. Integer programming. (Wiley, 1972).

Tip

This tutorial was generated using Literate.jl. View the source .jl file on GitHub.