Originally Contributed by: Shuvomoy Das Gupta
This notebook describes how to implement Benders decomposition in JuMP, which is a large scale optimization scheme. We only discuss the classical approach (using loops) here. The approach using lazy constraints is showed in the correponding notebook.
To illustrate an implementation of the Benders decomposition in JuMP, we apply it to the following general mixed integer problem:
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.
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 2.
Step 2 (Solving the master problem)
Solve the master problem, $f_\text{m}^{(k)}$ =
Let the maximizer corresponding to the objective value $f_\text{m}^{(k)}$ be denoted by $x^{(k)}$. Now there are three possibilities:
then the original proble is infeasible and sadly, we are done.
then we take $f_\text{m}^{(k)}$ to be arbitrarily large and $x^{(k)}$ to be a corresponding feasible solution. Go to Step 3
Step 3 (Solving the subproblem and add Benders cut when needed)
Solve the subproblem, $f_s(x^{(k)})$ =
Let the minimizer corresponding to the objective value $f_s(x^{(k)})$ be denoted by $u^{(k)}$. There are three possibilities:
We quit from Benders algorithm and use special purpose algorithm to find a feasible solution if there exists one.
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.
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!
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.
For a more general approach to Bender's Decomposition you can have a look at Mathieu Besançon's blog.
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;
There are two ways we can implement Benders decomposition in JuMP:
The classical approach might be inferior to the modern one, as the solver
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.
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 LinearAlgebra
using Test
# Master Problem Description
# --------------------------
master_problem_model = Model(GLPK.Optimizer)
# Variable Definition
# ----------------------------------------------------------------
@variable(master_problem_model, 0 <= x[1:dim_x] <= 1e6, Int)
@variable(master_problem_model, t <= 1e6)
# Objective Setting
# -----------------
@objective(master_problem_model, Max, t)
global iter_num = 1
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
Here is the loop that checks the status of the master problem and the subproblem and then adds necessary Benders cuts accordingly.
iter_num = 1
while(true)
println("\n-----------------------")
println("Iteration number = ", iter_num)
println("-----------------------\n")
println("The current master problem is")
print(master_problem_model)
optimize!(master_problem_model)
t_status = termination_status(master_problem_model)
p_status = primal_status(master_problem_model)
if p_status == MOI.INFEASIBLE_POINT
println("The problem is infeasible :-(")
break
end
(fm_current, x_current) = if t_status == MOI.INFEASIBLE_OR_UNBOUNDED
(M, M * ones(dim_x))
elseif p_status == MOI.FEASIBLE_POINT
(value(t), value.(x))
else
error("Unexpected status: $((t_status, p_status))")
end
println("Status of the master problem is ", t_status,
"\nwith fm_current = ", fm_current,
"\nx_current = ", x_current)
sub_problem_model = Model(GLPK.Optimizer)
c_sub = b - A1 * x_current
local u = @variable(sub_problem_model, u[1:dim_u] >= 0)
@constraint(sub_problem_model, constr_ref_subproblem[j = 1:size(A2, 2)], dot(A2[:, j], u) >= c2[j])
# The second argument of @constraint macro,
# constr_ref_subproblem[j=1:size(A2,2)] means that the j-th constraint is
# referenced by constr_ref_subproblem[j].
@objective(sub_problem_model, Min, dot(c1, x_current) + dot(c_sub, u))
print("\nThe current subproblem model is \n", sub_problem_model)
optimize!(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)
γ = dot(b, u_current)
println("Status of the subproblem is ", t_status_sub,
"\nwith fs_x_current = ", fs_x_current,
"\nand fm_current = ", fm_current)
if p_status_sub == MOI.FEASIBLE_POINT && fs_x_current == fm_current # we are done
println("\n################################################")
println("Optimal solution of the original problem found")
println("The optimal objective value t is ", fm_current)
println("The optimal x is ", x_current)
println("The optimal v is ", dual.(constr_ref_subproblem))
println("################################################\n")
break
end
if p_status_sub == MOI.FEASIBLE_POINT && fs_x_current < fm_current
println("\nThere is a suboptimal vertex, add the corresponding constraint")
cv = A1' * u_current - c1
@constraint(master_problem_model, t + dot(cv, x) <= γ)
println("t + ", cv, "ᵀ x <= ", γ)
end
if t_status_sub == MOI.INFEASIBLE_OR_UNBOUNDED
println("\nThere is an extreme ray, adding the corresponding constraint")
ce = A1'* u_current
@constraint(master_problem_model, dot(ce, x) <= γ)
println(ce, "ᵀ x <= ", γ)
end
global iter_num += 1
end
@test value(t) ≈ -4
----------------------- Iteration number = 1 ----------------------- The current master problem is 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 Status of the master problem is OPTIMAL with fm_current = 1.0e6 x_current = [0.0, 0.0] 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 Status of the subproblem is OPTIMAL with fs_x_current = -7.666666666666667 and fm_current = 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 Max t Subject to -x[1] - 4 x[2] + t ≤ -7.666666666666667 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 Status of the master problem is OPTIMAL with fm_current = 1.0e6 x_current = [0.0, 250002.0] The current subproblem model is Min 750004 u[1] + 750003 u[2] - 1.000008e6 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 Status of the subproblem is OPTIMAL with fs_x_current = -1.000008e6 and fm_current = 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 Max t Subject to -x[1] - 4 x[2] + t ≤ -7.666666666666667 x[1] + 4 x[2] + t ≤ 0.0 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 Status of the master problem is OPTIMAL with fm_current = -4.0 x_current = [4.0, 0.0] The current subproblem model is Min -6 u[1] + u[2] - 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 Status of the subproblem is OPTIMAL with fs_x_current = -13.0 and fm_current = -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 Max t Subject to -x[1] - 4 x[2] + t ≤ -7.666666666666667 x[1] + 4 x[2] + t ≤ 0.0 2.5 x[1] - 0.5 x[2] + t ≤ -3.0 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 Status of the master problem is OPTIMAL with fm_current = -4.0 x_current = [0.0, 1.0] 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 Status of the subproblem is OPTIMAL with fs_x_current = -4.0 and fm_current = -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] ################################################
Test Passed