{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "title: Benders Decomposition\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Originally Contributed by**: Shuvomoy Das Gupta" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This notebook describes how to implement Benders decomposition in JuMP, which is a large scale optimization scheme. \n", "We only discuss the classical approach (using loops) here.\n", "The approach using lazy constraints is showed in the correponding notebook." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To illustrate an implementation of the Benders decomposition in JuMP,\n", "we apply it to the following general mixed integer problem:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\\begin{align*}\n", "& \\text{maximize} \\quad &&c_1^T x+c_2^T v \\\\\n", "& \\text{subject to} \\quad &&A_1 x+ A_2 v \\preceq b \\\\\n", "& &&x \\succeq 0, x \\in \\mathbb{Z}^n \\\\\n", "& &&v \\succeq 0, v \\in \\mathbb{R}^p \\\\\n", "\\end{align*}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "where $b \\in \\mathbb{R}^m$, $A_1 \\in \\mathbb{R}^{m \\times n}$, $A_2 \\in \\mathbb{R}^{m \\times p}$ and\n", "\\mathbb{Z} is the set of integers. \n", "Here the symbol $\\succeq$ ($\\preceq$) stands for element-wise greater (less) than or equal to. \n", "Any mixed integer programming problem can be written in the form above." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We want to write the Benders decomposition algorithm for the problem above. \n", "Consider the polyhedron $\\{u \\in \\mathbb{R}^m| A_2^T u \\succeq 0, u \\succeq 0\\}$. \n", "Assume the set of vertices and extreme rays of the polyhedron is denoted by $P$ and $Q$ respectively. \n", "Assume on the $k$th iteration the subset of vertices of the polyhedron mentioned is denoted by \n", "$T(k)$ and the subset of extreme rays are denoted by $Q(k)$, \n", "which will be generated by the Benders decomposition problem below." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Benders decomposition algorithm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Step 1 (Initialization)**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We start with $T(1)=Q(1)=\\emptyset$. \n", "Let $f_m^{(1)}$ be arbitrarily large and $x^{(1)}$ be any non-negative integer vector and go to Step 2." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Step 2 (Solving the master problem)**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solve the master problem, $f_\\text{m}^{(k)}$ =" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\\begin{align*}\n", "&\\text{maximize} &&\\quad t \\\\\n", "&\\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} \\\\\n", "& && \\forall \\bar{y} \\in Q(k) \\qquad (A_1 ^T \\bar{y})^T x \\leq b^T \\bar{y} \\\\\n", "& && \\qquad \\qquad \\qquad \\; x \\succeq 0, x \\in \\mathbb{Z}^n\n", "\\end{align*}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let the maximizer corresponding to the objective value $f_\\text{m}^{(k)}$ be denoted by $x^{(k)}$. \n", "Now there are three possibilities:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- If $f_\\text{m}^{(k)}=-\\infty$, i.e., the master problem is infeasible, \n", "then the original proble is infeasible and sadly, we are done.\n", "\n", "- If $f_\\text{m}^{(k)}=\\infty$, i.e. the master problem is unbounded above, \n", "then we take $f_\\text{m}^{(k)}$ to be arbitrarily large and $x^{(k)}$ to be a corresponding feasible solution. \n", "Go to Step 3\n", "\n", "- If $f_\\text{m}^{(k)}$ is finite, then we collect $t^{(k)}$ and $x^{(k)}$ and go to Step 3." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Step 3 (Solving the subproblem and add Benders cut when needed)**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Solve the subproblem, $f_s(x^{(k)})$ =" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\\begin{align*}\n", " c_1^T x^{(k)} + & \\text{minimize} && (b-A_1 x^{(k)})^T u \\\\\n", " & \\text{subject to} && A_2^T u \\succeq c_2 \\\\\n", " & && u \\succeq 0, u \\in \\mathbb{R}^m\n", "\\end{align*}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let the minimizer corresponding to the objective value $f_s(x^{(k)})$ be denoted by $u^{(k)}$. There are three possibilities:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- If $f_s(x^{(k)}) = \\infty$, the original problem is either infeasible or unbounded. \n", "We quit from Benders algorithm and use special purpose algorithm to find a feasible solution if there exists one.\n", "\n", "- If $f_s(x^{(k)}) = - \\infty$, we arrive at an extreme ray $y^{(k)}$. \n", "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 \n", "i.e., $Q(k+1):= Q(k) \\cup \\{y^{(k)}\\}$. Take $k:=k+1$ and go to Step 3. \n", "\n", "- If $f_s(x^{(k)})$ is finite, then \n", "\n", " * If $f_s(x^{(k)})=f_m^{(k)}$ we arrive at the optimal solution. \n", "The optimum objective value of the original problem is $f_s(x^{(k)})=f_m^{(k)}$, \n", "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!\n", "\n", " * If $f_s(x^{(k)}) < f_m^{(k)}$ we get an suboptimal vertex $u^{(k)}$. \n", "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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For a more general approach to Bender's Decomposition you can have a look at \n", "[Mathieu Besançon's blog](https://matbesancon.github.io/post/2019-05-08-simple-benders/)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data for the problem" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The input data is from page 139, Integer programming by Garfinkel and Nemhauser[[1]](#c1)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "c1 = [-1; -4]\n", "c2 = [-2; -3]\n", "\n", "dim_x = length(c1)\n", "dim_u = length(c2)\n", "\n", "b = [-2; -3]\n", "\n", "A1 = [1 -3;\n", " -1 -3]\n", "A2 = [1 -2;\n", " -1 -1]\n", "\n", "M = 1000;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## How to implement the Benders decomposition algorithm in JuMP" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There are two ways we can implement Benders decomposition in JuMP: \n", "- *Classical approach:* Adding the Benders cuts in a loop, and\n", "- *Modern approach:* Adding the Benders cuts as lazy constraints." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The classical approach might be inferior to the modern one, as the solver\n", "- might revisit previously eliminated solution, and\n", "- might discard the optimal solution to the original problem in favor of a better but \n", "ultimately infeasible solution to the relaxed one." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For more details on the comparison between the two approaches, see [Paul Rubin's blog on Benders Decomposition](http://orinanobworld.blogspot.ca/2011/10/benders-decomposition-then-and-now.html)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Classical Approach: Adding the Benders Cuts in a Loop" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's describe the master problem first. Note that there are no constraints, which we will added later using Benders decomposition." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Max t\n", "Subject to\n", " x[1] ≥ 0.0\n", " x[2] ≥ 0.0\n", " x[1] ≤ 1.0e6\n", " x[2] ≤ 1.0e6\n", " t ≤ 1.0e6\n", " x[1] integer\n", " x[2] integer\n" ] } ], "source": [ "# Loading the necessary packages\n", "#-------------------------------\n", "using JuMP \n", "using GLPK\n", "using LinearAlgebra\n", "using Test\n", "\n", "# Master Problem Description\n", "# --------------------------\n", "\n", "master_problem_model = Model(GLPK.Optimizer)\n", "\n", "# Variable Definition \n", "# ----------------------------------------------------------------\n", "@variable(master_problem_model, 0 <= x[1:dim_x] <= 1e6, Int) \n", "@variable(master_problem_model, t <= 1e6)\n", "\n", "# Objective Setting\n", "# -----------------\n", "@objective(master_problem_model, Max, t)\n", "global iter_num = 1\n", "\n", "print(master_problem_model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here is the loop that checks the status of the master problem and the subproblem and \n", "then adds necessary Benders cuts accordingly." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "-----------------------\n", "Iteration number = 1\n", "-----------------------\n", "\n", "The current master problem is\n", "Max t\n", "Subject to\n", " x[1] ≥ 0.0\n", " x[2] ≥ 0.0\n", " x[1] ≤ 1.0e6\n", " x[2] ≤ 1.0e6\n", " t ≤ 1.0e6\n", " x[1] integer\n", " x[2] integer\n", "Status of the master problem is OPTIMAL\n", "with fm_current = 1.0e6\n", "x_current = [0.0, 0.0]\n", "\n", "The current subproblem model is \n", "Min -2 u[1] - 3 u[2]\n", "Subject to\n", " constr_ref_subproblem[1] : u[1] - u[2] ≥ -2.0\n", " constr_ref_subproblem[2] : -2 u[1] - u[2] ≥ -3.0\n", " u[1] ≥ 0.0\n", " u[2] ≥ 0.0\n", "Status of the subproblem is OPTIMAL\n", "with fs_x_current = -7.666666666666667\n", "and fm_current = 1.0e6\n", "\n", "There is a suboptimal vertex, add the corresponding constraint\n", "t + [-1.0, -4.0]ᵀ x <= -7.666666666666667\n", "\n", "-----------------------\n", "Iteration number = 2\n", "-----------------------\n", "\n", "The current master problem is\n", "Max t\n", "Subject to\n", " -x[1] - 4 x[2] + t ≤ -7.666666666666667\n", " x[1] ≥ 0.0\n", " x[2] ≥ 0.0\n", " x[1] ≤ 1.0e6\n", " x[2] ≤ 1.0e6\n", " t ≤ 1.0e6\n", " x[1] integer\n", " x[2] integer\n", "Status of the master problem is OPTIMAL\n", "with fm_current = 1.0e6\n", "x_current = [0.0, 250002.0]\n", "\n", "The current subproblem model is \n", "Min 750004 u[1] + 750003 u[2] - 1.000008e6\n", "Subject to\n", " constr_ref_subproblem[1] : u[1] - u[2] ≥ -2.0\n", " constr_ref_subproblem[2] : -2 u[1] - u[2] ≥ -3.0\n", " u[1] ≥ 0.0\n", " u[2] ≥ 0.0\n", "Status of the subproblem is OPTIMAL\n", "with fs_x_current = -1.000008e6\n", "and fm_current = 1.0e6\n", "\n", "There is a suboptimal vertex, add the corresponding constraint\n", "t + [1.0, 4.0]ᵀ x <= 0.0\n", "\n", "-----------------------\n", "Iteration number = 3\n", "-----------------------\n", "\n", "The current master problem is\n", "Max t\n", "Subject to\n", " -x[1] - 4 x[2] + t ≤ -7.666666666666667\n", " x[1] + 4 x[2] + t ≤ 0.0\n", " x[1] ≥ 0.0\n", " x[2] ≥ 0.0\n", " x[1] ≤ 1.0e6\n", " x[2] ≤ 1.0e6\n", " t ≤ 1.0e6\n", " x[1] integer\n", " x[2] integer\n", "Status of the master problem is OPTIMAL\n", "with fm_current = -4.0\n", "x_current = [4.0, 0.0]\n", "\n", "The current subproblem model is \n", "Min -6 u[1] + u[2] - 4\n", "Subject to\n", " constr_ref_subproblem[1] : u[1] - u[2] ≥ -2.0\n", " constr_ref_subproblem[2] : -2 u[1] - u[2] ≥ -3.0\n", " u[1] ≥ 0.0\n", " u[2] ≥ 0.0\n", "Status of the subproblem is OPTIMAL\n", "with fs_x_current = -13.0\n", "and fm_current = -4.0\n", "\n", "There is a suboptimal vertex, add the corresponding constraint\n", "t + [2.5, -0.5]ᵀ x <= -3.0\n", "\n", "-----------------------\n", "Iteration number = 4\n", "-----------------------\n", "\n", "The current master problem is\n", "Max t\n", "Subject to\n", " -x[1] - 4 x[2] + t ≤ -7.666666666666667\n", " x[1] + 4 x[2] + t ≤ 0.0\n", " 2.5 x[1] - 0.5 x[2] + t ≤ -3.0\n", " x[1] ≥ 0.0\n", " x[2] ≥ 0.0\n", " x[1] ≤ 1.0e6\n", " x[2] ≤ 1.0e6\n", " t ≤ 1.0e6\n", " x[1] integer\n", " x[2] integer\n", "Status of the master problem is OPTIMAL\n", "with fm_current = -4.0\n", "x_current = [0.0, 1.0]\n", "\n", "The current subproblem model is \n", "Min u[1] - 4\n", "Subject to\n", " constr_ref_subproblem[1] : u[1] - u[2] ≥ -2.0\n", " constr_ref_subproblem[2] : -2 u[1] - u[2] ≥ -3.0\n", " u[1] ≥ 0.0\n", " u[2] ≥ 0.0\n", "Status of the subproblem is OPTIMAL\n", "with fs_x_current = -4.0\n", "and fm_current = -4.0\n", "\n", "################################################\n", "Optimal solution of the original problem found\n", "The optimal objective value t is -4.0\n", "The optimal x is [0.0, 1.0]\n", "The optimal v is [0.0, 0.0]\n", "################################################\n", "\n" ] }, { "data": { "text/plain": [ "\u001b[32m\u001b[1mTest Passed\u001b[22m\u001b[39m" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "iter_num = 1\n", "\n", "while(true)\n", " println(\"\\n-----------------------\")\n", " println(\"Iteration number = \", iter_num)\n", " println(\"-----------------------\\n\")\n", " println(\"The current master problem is\")\n", " print(master_problem_model)\n", " \n", " optimize!(master_problem_model)\n", " \n", " t_status = termination_status(master_problem_model)\n", " p_status = primal_status(master_problem_model)\n", " \n", " if p_status == MOI.INFEASIBLE_POINT\n", " println(\"The problem is infeasible :-(\")\n", " break\n", " end\n", "\n", " (fm_current, x_current) = if t_status == MOI.INFEASIBLE_OR_UNBOUNDED\n", " (M, M * ones(dim_x)) \n", " elseif p_status == MOI.FEASIBLE_POINT\n", " (value(t), value.(x))\n", " else\n", " error(\"Unexpected status: $((t_status, p_status))\")\n", " end\n", "\n", " println(\"Status of the master problem is \", t_status, \n", " \"\\nwith fm_current = \", fm_current, \n", " \"\\nx_current = \", x_current)\n", "\n", " sub_problem_model = Model(GLPK.Optimizer)\n", "\n", " c_sub = b - A1 * x_current\n", "\n", " local u = @variable(sub_problem_model, u[1:dim_u] >= 0)\n", "\n", " @constraint(sub_problem_model, constr_ref_subproblem[j = 1:size(A2, 2)], dot(A2[:, j], u) >= c2[j])\n", " # The second argument of @constraint macro,\n", " # constr_ref_subproblem[j=1:size(A2,2)] means that the j-th constraint is\n", " # referenced by constr_ref_subproblem[j].\n", " \n", " @objective(sub_problem_model, Min, dot(c1, x_current) + dot(c_sub, u))\n", "\n", " print(\"\\nThe current subproblem model is \\n\", sub_problem_model)\n", "\n", " optimize!(sub_problem_model)\n", "\n", " t_status_sub = termination_status(sub_problem_model)\n", " p_status_sub = primal_status(sub_problem_model)\n", "\n", " fs_x_current = objective_value(sub_problem_model) \n", "\n", " u_current = value.(u)\n", "\n", " γ = dot(b, u_current)\n", "\n", " println(\"Status of the subproblem is \", t_status_sub, \n", " \"\\nwith fs_x_current = \", fs_x_current, \n", " \"\\nand fm_current = \", fm_current) \n", " \n", " if p_status_sub == MOI.FEASIBLE_POINT && fs_x_current == fm_current # we are done\n", " println(\"\\n################################################\")\n", " println(\"Optimal solution of the original problem found\")\n", " println(\"The optimal objective value t is \", fm_current)\n", " println(\"The optimal x is \", x_current)\n", " println(\"The optimal v is \", dual.(constr_ref_subproblem))\n", " println(\"################################################\\n\")\n", " break\n", " end \n", " \n", " if p_status_sub == MOI.FEASIBLE_POINT && fs_x_current < fm_current\n", " println(\"\\nThere is a suboptimal vertex, add the corresponding constraint\")\n", " cv = A1' * u_current - c1\n", " @constraint(master_problem_model, t + dot(cv, x) <= γ)\n", " println(\"t + \", cv, \"ᵀ x <= \", γ)\n", " end \n", " \n", " if t_status_sub == MOI.INFEASIBLE_OR_UNBOUNDED\n", " println(\"\\nThere is an extreme ray, adding the corresponding constraint\")\n", " ce = A1'* u_current\n", " @constraint(master_problem_model, dot(ce, x) <= γ)\n", " println(ce, \"ᵀ x <= \", γ)\n", " end\n", " \n", " global iter_num += 1\n", "end\n", "\n", "@test value(t) ≈ -4" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### References\n", "\n", "1. Garfinkel, R. & Nemhauser, G. L. Integer programming. (Wiley, 1972)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 1.5.1", "language": "julia", "name": "julia-1.5" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.5.1" } }, "nbformat": 4, "nbformat_minor": 2 }