{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Disclaimer\n", "This notebook is only working under the versions:\n", "\n", "- JuMP 0.19 (unreleased, but currently in master)\n", "\n", "- MathOptInterface 0.4.1\n", "\n", "- GLPK 0.6.0\n", "\n", "## Disclaimer 2\n", "\n", "The second part od this notebook (using Lazy Constraints) is not working" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**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.\n", "\n", "**Author:** [Shuvomoy Das Gupta](http://scg.utoronto.ca/~shuvomoy.dasgupta/)\n", "\n", "**License:** \"Creative
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.\n", "\n", "# Using Julia+JuMP for optimization - Benders decomposition\n", "\n", "\n", "--------------------------\n", "
\n", "\n", "To illustrate implementation of solver callback in JuMP, we consider applying Benders decomposition to the following general mixed integer problem:\n", "\n", "\\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}\n", "\n", "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.\n", "\n", "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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Benders decomposition algorithm\n", "---------------------------------------\n", "\n", "**Step 1 (Initialization)**
\n", "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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "** Step 2 (Solving the master problem)**
\n", "Solve the master problem:\n", "\n", "\n", "\\begin{align}\n", "& f_\\text{m}^{(k)}= \\\\\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}\n", "\n", "Let the maximizer corresponding to the objective value $f_\\text{m}^{(k)}$ be denoted by $x^{(k)}$. Now there are three possibilities:\n", "\n", "- If $f_\\text{m}^{(k)}=-\\infty$, i.e., the master problem is infeasible, 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, then we take $f_\\text{m}^{(k)}$ to be arbitrarily large and $x^{(k)}$ to be a corresponding feasible solution. 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.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "** Step 3 (Solving the subproblem and add Benders cut when needed) **
\n", "Solve the subproblem\n", "\n", "\\begin{align}\n", " f_s(x^{(k)})= \\\\\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}\n", "\n", "Let the minimizer corresponding to the objective value $f_s(x^{(k)})$ be denoted by $u^{(k)}$. There are three possibilities:\n", "\n", "- 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.\n", "\n", "- 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.\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. 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!\n", " \n", " * 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.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data for the problem\n", "\n", "The input data is from page 139, Integer programming by Garfinkel and Nemhauser, 1972." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "1000" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Data for the problem\n", "#---------------------\n", "c1=[-1;-4]\n", "c2=[-2; -3]\n", "dimX=length(c1)\n", "dimU=length(c2)\n", "b=[-2;-3]\n", "A1=[\n", " 1 -3;\n", " -1 -3\n", " ]\n", "A2=[\n", " 1 -2;\n", " -1 -1\n", " ]\n", "M=1000" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## How to implement the Benders decomposition algorithm in JuMP\n", "\n", "There are two ways we can implement Benders decomposition in JuMP:\n", "\n", "- *Classical approach:* Adding the Benders cuts in a loop, and\n", "- *Modern approach:* Adding the Benders cuts as lazy constraints.\n", "\n", "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 ultimately infeasible solution to the relaxed one.\n", "\n", "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": 16, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "A JuMP Model" ] } ], "source": [ "# Loading the necessary packages\n", "#-------------------------------\n", "using JuMP \n", "using GLPK\n", "using MathOptInterface\n", "const MOI = MathOptInterface\n", "\n", "# Master Problem Description\n", "# --------------------------\n", "\n", "masterProblemModel = Model(optimizer=GLPK.GLPKOptimizerMIP()) \n", "\n", "# Variable Definition \n", "# ----------------------------------------------------------------\n", "@variable(masterProblemModel, 0<= x[1:dimX]<= 1e6 , Int) \n", "@variable(masterProblemModel, t<=1e6)\n", "\n", "# Objective Setting\n", "# -----------------\n", "@objective(masterProblemModel, Max, t)\n", "iC=1 # iC is the iteration counter \n", "\n", "print(masterProblemModel)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Here is the loop that checks the status of the master problem and the subproblem and then adds necessary Benders cuts accordingly." ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "-----------------------\n", "Iteration number = 1\n", "-----------------------\n", "\n", "The current master problem is\n", "A JuMP ModelStatus of the master problem isSuccess\n", "with fmCurrent = 1.0e6\n", "xCurrent = [0.0, 0.0]\n", "\n", "The current subproblem model is \n", "A JuMP ModelStatus of the subproblem is Success\n", "with fsxCurrent= -7.666666666666667\n", "and fmCurrent= 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", "A JuMP ModelStatus of the master problem isSuccess\n", "with fmCurrent = 1.0e6\n", "xCurrent = [0.0, 250002.0]\n", "\n", "The current subproblem model is \n", "A JuMP ModelStatus of the subproblem is Success\n", "with fsxCurrent= -1.000008e6\n", "and fmCurrent= 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", "A JuMP ModelStatus of the master problem isSuccess\n", "with fmCurrent = -4.0\n", "xCurrent = [4.0, 0.0]\n", "\n", "The current subproblem model is \n", "A JuMP ModelStatus of the subproblem is Success\n", "with fsxCurrent= -13.0\n", "and fmCurrent= -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", "A JuMP ModelStatus of the master problem isSuccess\n", "with fmCurrent = -4.0\n", "xCurrent = [0.0, 1.0]\n", "\n", "The current subproblem model is \n", "A JuMP ModelStatus of the subproblem is Success\n", "with fsxCurrent= -4.0\n", "and fmCurrent= -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" ] } ], "source": [ "# Trying the entire benders decomposition in one step\n", "while(true)\n", " println(\"\\n-----------------------\")\n", " println(\"Iteration number = \", iC)\n", " println(\"-----------------------\\n\")\n", " println(\"The current master problem is\")\n", " print(masterProblemModel)\n", " \n", " JuMP.optimize(masterProblemModel)\n", " t_status = JuMP.terminationstatus(masterProblemModel)# == MOI.Success\n", " p_status = JuMP.primalstatus(masterProblemModel)# == MOI.FeasiblePoint\n", " \n", " if p_status == MOI.InfeasiblePoint\n", " println(\"The problem is infeasible :-(\")\n", " break\n", " end\n", "\n", " if t_status == MOI.InfeasibleOrUnbounded\n", " fmCurrent = M\n", " xCurrent=M*ones(dimX)\n", " end\n", "\n", "\n", " if p_status == MOI.FeasiblePoint\n", " fmCurrent = JuMP.resultvalue(t)\n", " xCurrent=Float64[]\n", " for i in 1:dimX\n", " push!(xCurrent, JuMP.resultvalue(x[i]))\n", " end\n", " end\n", "\n", " println(\"Status of the master problem is\", t_status, \n", " \"\\nwith fmCurrent = \", fmCurrent, \n", " \"\\nxCurrent = \", xCurrent)\n", "\n", "\n", " subProblemModel = Model(optimizer=GLPK.GLPKOptimizerLP())\n", "\n", " cSub=b-A1*xCurrent\n", "\n", " @variable(subProblemModel, u[1:dimU]>=0)\n", "\n", "\n", " @constraint(subProblemModel, constrRefSubProblem[j=1:size(A2,2)],sum(A2[i,j]*u[i] for i in 1:size(A2,1))>=c2[j])\n", " # The second argument of @constraint macro, constrRefSubProblem[j=1:size(A2,2)] means that the j-th constraint is\n", " # referenced by constrRefSubProblem[j]. \n", "\n", " \n", " @objective(subProblemModel, Min, dot(c1, xCurrent) + sum(cSub[i]*u[i] for i in 1:dimU))\n", "\n", " print(\"\\nThe current subproblem model is \\n\", subProblemModel)\n", "\n", " JuMP.optimize(subProblemModel) \n", " t_status_sub = JuMP.terminationstatus(subProblemModel)# == MOI.Success\n", " p_status_sub = JuMP.primalstatus(subProblemModel)# == MOI.FeasiblePoint\n", "\n", " fsxCurrent = JuMP.objectivevalue(subProblemModel) \n", "\n", " uCurrent = Float64[]\n", "\n", " for i in 1:dimU\n", " push!(uCurrent, JuMP.resultvalue(u[i]))\n", " end\n", "\n", " γ=dot(b,uCurrent)\n", "\n", " println(\"Status of the subproblem is \", t_status_sub, \n", " \"\\nwith fsxCurrent= \", fsxCurrent, \n", " \"\\nand fmCurrent= \", fmCurrent) \n", " \n", " if p_status_sub == MOI.FeasiblePoint && fsxCurrent == fmCurrent # we are done\n", " println(\"\\n################################################\")\n", " println(\"Optimal solution of the original problem found\")\n", " println(\"The optimal objective value t is \", fmCurrent)\n", " println(\"The optimal x is \", xCurrent)\n", " println(\"The optimal v is \", JuMP.resultdual.(constrRefSubProblem))\n", " println(\"################################################\\n\")\n", " break\n", " end \n", " \n", " if p_status_sub == MOI.FeasiblePoint && fsxCurrent < fmCurrent \n", " println(\"\\nThere is a suboptimal vertex, add the corresponding constraint\")\n", " cv= A1'*uCurrent - c1\n", " @constraint(masterProblemModel, t+sum(cv[i]*x[i] for i in 1:dimX) <= γ )\n", " println(\"t + \", cv, \"ᵀ x <= \", γ)\n", " end \n", " \n", " if t_status_sub == MOI.InfeasibleOrUnbounded\n", " println(\"\\nThere is an extreme ray, adding the corresponding constraint\")\n", " ce = A1'*uCurrent\n", " @constraint(masterProblemModel, sum(ce[i]*x[i] for i in 1:dimX) <= γ)\n", " println(ce, \"ᵀ x <= \", γ)\n", " end\n", " \n", " iC=iC+1\n", " sleep(0.5)\n", " \n", "end" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Modern approach: adding the Benders cuts as lazy constraints" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### What are lazy constraints?\n", "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. \n", "\n", "Lazy constraints have the following advantages:\n", "- does not revisit previously eliminated solution, and\n", "- does not discard the optimal solution to the original problem in favor of a better but ultimately infeasible solution to the relaxed one." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### How to add lazy constraints in Julia?\n", "\n", "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.\n", "\n", "When we add the lazy constraints we have to use the `@addLazyConstraint` macro as follows:\n", "\n", ">``@addLazyConstraint(cb, description of the constraint)\n", "``.\n", "\n", "Description of the constraint will be of the same manner as in adding a normal constraint in JuMP using `@constraint` macro. \n", "\n", "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:\n", "\n", ">``addLazyCallback(name of the master problem model, addBendersCut)\n", "``\n", "\n", ", which tells JuMP to add the lazy constraints generated by the function `addBendersCut` to the master problem." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ " # Loading the necessary packages\n", " #-------------------------------\n", "using JuMP \n", "#using CPLEX\n", "using Gurobi\n", "\n", " # Master Problem Description\n", " # --------------------------\n", "\n", "# Model name\n", " \n", "#masterProblemModel = Model(solver = CplexSolver())\n", "masterProblemModel = Model(optimizer = GurobiOptimizer(Heuristics=0, Cuts = 0, OutputFlag = 0)) # If we want to add Benders lazy constraints \n", "# in Gurobi, then we have to turn off Gurobi's own Cuts and Heuristics in the master problem\n", "\n", "# Variable Definition (Only CplexSolver() works properly for these)\n", "# ----------------------------------------------------------------\n", "#@variable(masterProblemModel, x[1:dimX] >=0 , Int) \n", "#@variable(masterProblemModel, t)\n", "\n", "# ***************ALTERNATIVE VARIABLE DEFINITION FOR GUROBI************\n", "#If we replace the two lines above with the follwoing:\n", "@variable(masterProblemModel, 0<= x[1:dimX] <= 1e6 , Int)\n", "@variable(masterProblemModel, t <= 1e6)\n", "# then all the solvers give the expected solution\n", "#**********************************************************************\n", "\n", "# Objective Setting\n", "# -----------------\n", "@objective(masterProblemModel, Max, t)\n", "\n", "print(masterProblemModel)\n", "\n", "stringOfBenderCuts=AbstractString[] # this is an array of strings which will contain all the\n", "# Benders cuts added to be displayed later\n", "\n", "# There are no constraints when we start, so we will add all the constraints in the\n", "# form of Benders cuts lazily" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "function addBendersCut(cb)\n", " #***************************************************************************\n", " # First we store the master problem solution in conventional data structures\n", " println(\"----------------------------\")\n", " println(\"ITERATION NUMBER = \", length(stringOfBenderCuts)+1)\n", " println(\"---------------------------\\n\")\n", " \n", " fmCurrent = getvalue(t)\n", " xCurrent=Float64[]\n", " for i in 1:dimX\n", " push!(xCurrent, JuMP.resultvalue(x[i]))\n", " end\n", " \n", " # Display the current solution of the master problem\n", " println(\"MASTERPROBLEM INFORMATION\")\n", " println(\"-------------------------\")\n", " println(\"The master problem that was solved was:\")\n", " print(masterProblemModel)\n", " println(\"with \", length(stringOfBenderCuts), \" added lazy constraints\")\n", " println(stringOfBenderCuts)\n", " println(\"Current Value of x is: \", xCurrent)\n", " println(\"Current objective value of master problem, fmCurrent is: \", fmCurrent)\n", " println(\"\\n\")\n", " \n", " #************************************************************************\n", " \n", " # ========================================================================\n", " # Now we solve the subproblem\n", "\n", " # subProblemModel=Model(solver=CplexSolver())\n", "\n", " subProblemModel = Model(optimizer=GurobiOptimizer(Presolve=0, OutputFlag = 0))\n", " \n", " cSub=b-A1*xCurrent\n", " \n", " @variable(subProblemModel, u[1:dimU]>=0)\n", " \n", "\n", " @constraint(subProblemModel, constrRefSubProblem[j=1:size(A2,2)], sum(A2[i,j]*u[i] for i in 1:size(A2,1))>=c2[j])\n", "\n", " \n", " @objective(subProblemModel, Min, dot(c1, xCurrent) + sum(cSub[i]*u[i] for i in 1:dimU))\n", " \n", " println(\"The subproblem is being solved\")\n", " \n", " statusSubProblem = JuMP.optimize(subProblemModel) \n", "\n", " # We store the results achieved from the subproblem in conventional data structures \n", " \n", " fsxCurrent = JuMP.objectivevalue(subProblemModel) \n", "\n", " uCurrent = Float64[]\n", " for i in 1:dimU\n", " push!(uCurrent, JuMP.resultvalue(u[i]))\n", " end\n", " \n", " # Display the solution corresponding to the subproblem\n", " \n", " println(\"SUBPROBLEM INFORMATION\")\n", " println(\"----------------------\")\n", " println(\"The subproblem that was solved was: \")\n", " print(subProblemModel)\n", " println(\"Current status of the subproblem is \", statusSubProblem)\n", " println(\"Current Value of u is: \", uCurrent) # JuMP will return an extreme ray\n", " # automatically (if the solver supports it), so we do not need to change the syntax\n", " println(\"Current Value of fs(xCurrent) is: \", fsxCurrent)\n", " println(\"\\n\")\n", " \n", " # ==========================================================================\n", " \n", " \n", " \n", " # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n", " # Now we check the status of the algorithm and add Benders cut when necessary\n", " γ=dot(b,uCurrent)\n", " \n", "\n", " \n", " if statusSubProblem == :Optimal && fsxCurrent==fmCurrent # we are done\n", " println(\"OPTIMAL SOLUTION OF THE ORIGINAL PROBLEM FOUND :-)\")\n", " println(\"The optimal objective value t is \", fmCurrent)\n", " println(\"The optimal x is \", xCurrent)\n", " println(\"The optimal v is \", JuMP.resultdual(constrRefSubProblem))\n", " println(\"\\n\")\n", " return\n", " end \n", "\n", " println(\"-------------------ADDING LAZY CONSTRAINT----------------\") \n", " if statusSubProblem == :Optimal && fsxCurrent < fmCurrent \n", " println(\"\\nThere is a suboptimal vertex, add the corresponding constraint\")\n", " cv= A1'*uCurrent - c1\n", " @lazyconstraint(cb, t+sum(cv[i]*x[i] for i in 1:dimX) <= γ)\n", " println(\"t + \", cv, \"ᵀ x <= \", γ)\n", " push!(stringOfBenderCuts, string(\"t+\", cv, \"'x <=\", γ))\n", " end\n", " \n", " if statusSubProblem == :Unbounded \n", " println(\"\\nThere is an extreme ray, adding the corresponding constraint\")\n", " ce = A1'*uCurrent\n", " @lazyconstraint(cb, sum(ce[i]*x[i] for i in 1:dimX) <= γ)\n", " println(ce, \"x <= \", γ)\n", " push!(stringOfBenderCuts, string(ce, \"ᵀ x <= \", γ))\n", " end\n", " println(\"\\n\") \n", " #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n", " \n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we tell the solver to use the callback function and solve the problem." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "addlazycallback(masterProblemModel, addBendersCut) # Telling the solver to use the \n", "# callback function" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false, "scrolled": false }, "outputs": [], "source": [ "JuMP.optimize(masterProblemModel)" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 0.6.0", "language": "julia", "name": "julia-0.6" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.6.0" } }, "nbformat": 4, "nbformat_minor": 0 }