diff --git a/notebooks/Shuvomoy - Benders decomposition.ipynb b/notebooks/Shuvomoy - Benders decomposition.ipynb new file mode 100644 index 0000000..a8c664f --- /dev/null +++ b/notebooks/Shuvomoy - Benders decomposition.ipynb @@ -0,0 +1,976 @@ +{ + "metadata": { + "language": "Julia", + "name": "", + "signature": "sha256:67099c03592aa05ccaa21fc95b9d604602d2abac82fcf313a8c03cabb1e4be41" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "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 nonnegative 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 arbitraily 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 feasbile 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" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + " # 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\n", + "\n" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 10, + "text": [ + "1000" + ] + } + ], + "prompt_number": 10 + }, + { + "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", + "collapsed": false, + "input": [ + "# Loading the necessary packages\n", + "#-------------------------------\n", + "using JuMP \n", + "using GLPKMathProgInterface\n", + "\n", + "# Master Problem Description\n", + "# --------------------------\n", + "\n", + "\n", + "masterProblemModel = Model(solver=GLPKSolverMIP()) \n", + "\n", + "# Variable Definition \n", + "# ----------------------------------------------------------------\n", + "@defVar(masterProblemModel, 0<= x[1:dimX]<= 1e6 , Int) \n", + "@defVar(masterProblemModel, t<=1e6)\n", + "\n", + "# Objective Setting\n", + "# -----------------\n", + "@setObjective(masterProblemModel, Max, t)\n", + "iC=1 # iC is the iteration counter \n", + "\n", + "print(masterProblemModel)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Max t\n", + "Subject to\n", + " 0 <= x[i] <= 1.0e6, integer, for all i in {1,2}\n", + " t <= 1.0e6\n" + ] + } + ], + "prompt_number": 11 + }, + { + "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", + "collapsed": false, + "input": [ + "# 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", + " statusMasterProblem = solve(masterProblemModel)\n", + " \n", + "if statusMasterProblem == :Infeasible\n", + " println(\"The problem is infeasible :-(\")\n", + " break\n", + "end\n", + " \n", + "if statusMasterProblem == :Unbounded \n", + " fmCurrent = M\n", + " xCurrent=M*ones(dimX)\n", + "end\n", + "\n", + "\n", + "if statusMasterProblem == :Optimal\n", + " fmCurrent = getValue(t)\n", + " xCurrent=Float64[]\n", + " for i in 1:dimX\n", + " push!(xCurrent,getValue(x[i]))\n", + " end\n", + "end\n", + "\n", + "println(\"Status of the master problem is\", statusMasterProblem, \n", + " \"\\nwith fmCurrent = \", fmCurrent, \n", + " \"\\nxCurrent = \", xCurrent)\n", + "\n", + " \n", + "subProblemModel = Model(solver=GLPKSolverLP())\n", + " \n", + "cSub=b-A1*xCurrent\n", + " \n", + "@defVar(subProblemModel, u[1:dimU]>=0)\n", + " \n", + "\n", + "@addConstraint(subProblemModel, constrRefSubProblem[j=1:size(A2,2)],sum{A2[i,j]*u[i] , i in [1:size(A2,1)]}>=c2[j])\n", + "# The second arguemnt of @addConstraint macro, constrRefSubProblem[j=1:size(A2,2)] means that the j-th constraint is\n", + "# referenced by constrRefSubProblem[j]. \n", + "\n", + " \n", + "@setObjective(subProblemModel, Min, dot(c1, xCurrent) + sum{cSub[i]*u[i], i in [1:dimU]})\n", + "\n", + "print(\"\\nThe current subproblem model is \\n\", subProblemModel)\n", + " \n", + "statusSubProblem = solve(subProblemModel) \n", + " \n", + "fsxCurrent = getObjectiveValue(subProblemModel) \n", + "\n", + "uCurrent = Float64[]\n", + " \n", + "for i in 1:dimU\n", + " push!(uCurrent, getValue(u[i]))\n", + "end\n", + "\n", + "\u03b3=dot(b,uCurrent)\n", + "\n", + " println(\"Status of the subproblem is \", statusSubProblem, \n", + " \"\\nwith fsxCurrent= \", fsxCurrent, \n", + " \"\\nand fmCurrent= \", fmCurrent) \n", + " \n", + " if statusSubProblem == :Optimal && 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 \", getDual(constrRefSubProblem))\n", + " println(\"################################################\\n\")\n", + " break\n", + " end \n", + " \n", + " if statusSubProblem == :Optimal && fsxCurrent < fmCurrent \n", + " println(\"\\nThere is a suboptimal vertex, add the corresponding constraint\")\n", + " cv= A1'*uCurrent - c1\n", + " @addConstraint(masterProblemModel, t+sum{cv[i]*x[i], i in 1:dimX} <= \u03b3 )\n", + " println(\"t + \", cv, \"\u1d40 x <= \", \u03b3)\n", + " end \n", + " \n", + " if statusSubProblem == :Unbounded \n", + " println(\"\\nThere is an extreme ray, adding the corresponding constraint\")\n", + " ce = A1'*uCurrent\n", + " @addConstraint(masterProblemModel, sum{ce[i]*x[i], i in 1:dimX} <= \u03b3)\n", + " println(ce, \"\u1d40 x <= \", \u03b3)\n", + " end\n", + " \n", + " iC=iC+1\n", + " \n", + "end" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n", + "-----------------------" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n", + "Iteration number = 1\n", + "-----------------------\n", + "\n", + "The current master problem is\n", + "Max t\n", + "Subject to\n", + " 0 <= x[i] <= 1.0e6, integer, for all i in {1,2}\n", + " t <= 1.0e6\n", + "Status of the master problem isOptimal\n", + "with fmCurrent = 1.0e6\n", + "xCurrent = [0.0,0.0]\n" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n", + "The current subproblem model is \n" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Min -2 u[1] - 3 u[2]\n", + "Subject to\n", + " u[1] - u[2] >= -2\n", + " -2 u[1] - u[2] >= -3\n", + " u[i] >= 0 for all i in {1,2}\n", + "Status of the subproblem is Optimal\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]\u1d40 x <= -7.666666666666667\n", + "\n", + "-----------------------\n", + "Iteration number = 2\n", + "-----------------------\n", + "\n", + "The current master problem is\n", + "Max t\n", + "Subject to\n", + " t - x[1] - 4 x[2] <= -7.666666666666667\n", + " 0 <= x[i] <= 1.0e6, integer, for all i in {1,2}\n", + " t <= 1.0e6\n", + "Status of the master problem isOptimal\n", + "with fmCurrent = 1.0e6\n", + "xCurrent = [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", + " u[1] - u[2] >= -2\n", + " -2 u[1] - u[2] >= -3\n", + " u[i] >= 0 for all i in {1,2}\n", + "Status of the subproblem is Optimal\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]\u1d40 x <= -0.0\n", + "\n", + "-----------------------\n", + "Iteration number = 3\n", + "-----------------------\n", + "\n", + "The current master problem is\n", + "Max t\n", + "Subject to\n", + " t - x[1] - 4 x[2] <= -7.666666666666667\n", + " t + x[1] + 4 x[2] <= 0\n", + " 0 <= x[i] <= 1.0e6, integer, for all i in {1,2}\n", + " t <= 1.0e6\n", + "Status of the master problem isOptimal\n", + "with fmCurrent = -4.0\n", + "xCurrent = [4.0,0.0]\n", + "\n", + "The current subproblem model is \n", + "Min -6 u[1] + u[2] - 4\n", + "Subject to\n", + " u[1] - u[2] >= -2\n", + " -2 u[1] - u[2] >= -3\n", + " u[i] >= 0 for all i in {1,2}\n", + "Status of the subproblem is Optimal\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]\u1d40 x <= -3.0\n", + "\n", + "-----------------------\n", + "Iteration number = 4\n", + "-----------------------\n", + "\n", + "The current master problem is\n", + "Max t\n", + "Subject to\n", + " t - x[1] - 4 x[2] <= -7.666666666666667\n", + " t + x[1] + 4 x[2] <= 0\n", + " t + 2.5 x[1] - 0.5 x[2] <= -3\n", + " 0 <= x[i] <= 1.0e6, integer, for all i in {1,2}\n", + " t <= 1.0e6\n", + "Status of the master problem isOptimal\n", + "with fmCurrent = -4.0\n", + "xCurrent = [0.0,1.0]\n", + "\n", + "The current subproblem model is \n", + "Min u[1] - 4\n", + "Subject to\n", + " u[1] - u[2] >= -2\n", + " -2 u[1] - u[2] >= -3\n", + " u[i] >= 0 for all i in {1,2}\n", + "Status of the subproblem is Optimal\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" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "The optimal v is " + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "constrRefSubProblem: 1 dimensions:\n", + "[1] = 0.0\n", + "[2] = 0.0\n", + "\n", + "################################################\n", + "\n" + ] + } + ], + "prompt_number": 12 + }, + { + "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 managable 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. " + ] + }, + { + "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 `@addConstraint` 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", + ">``setLazyCallback(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", + "collapsed": false, + "input": [ + " # 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(solver = GurobiSolver(Heuristics=0.0, Cuts = 0)) # If we want to add Benders lazy constraints \n", + "# in Gurobi, then we have to turn of Gurobi's own Cuts and Heuristics in the master problem\n", + "\n", + "# Variable Definition (Only CplexSolver() works properly for these)\n", + "# ----------------------------------------------------------------\n", + "@defVar(masterProblemModel, x[1:dimX] >=0 , Int) \n", + "@defVar(masterProblemModel, t)\n", + "\n", + "# ***************ALTERNATIVE VARIABLE DEFINITION FOR GUROBI************\n", + "#If we replace the two lines above with the follwoing:\n", + "#@defVar(masterProblemModel, 0<= x[1:dimX] <= 1e6 , Int)\n", + "#@defVar(masterProblemModel, t <= 1e6)\n", + "# then all the solvers give the expected solution\n", + "#**********************************************************************\n", + "\n", + "# Objective Setting\n", + "# -----------------\n", + "@setObjective(masterProblemModel, Max, t)\n", + "\n", + "print(masterProblemModel)\n", + "\n", + "stringOfBenderCuts=String[] # 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" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Max t\n", + "Subject to\n", + " x[i] >= 0, integer, for all i in {1,2}\n", + " t free\n" + ] + }, + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 13, + "text": [ + "0-element Array{String,1}" + ] + } + ], + "prompt_number": 13 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "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,getValue(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(solver=GurobiSolver(Presolve=0.0))\n", + " \n", + " cSub=b-A1*xCurrent\n", + " \n", + " @defVar(subProblemModel, u[1:dimU]>=0)\n", + " \n", + "\n", + " @addConstraint(subProblemModel, constrRefSubProblem[j=1:size(A2,2)], sum{A2[i,j]*u[i] , i in [1:size(A2,1)]}>=c2[j])\n", + "\n", + " \n", + " @setObjective(subProblemModel, Min, dot(c1, xCurrent) + sum{cSub[i]*u[i], i in [1:dimU]})\n", + " \n", + " println(\"The subproblem is being solved\")\n", + " \n", + " statusSubProblem = solve(subProblemModel) \n", + "\n", + " # We store the results achieved from the subproblem in conventional data structures \n", + " \n", + " fsxCurrent = getObjectiveValue(subProblemModel) \n", + "\n", + " uCurrent = Float64[]\n", + " for i in 1:dimU\n", + " push!(uCurrent, getValue(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", + " \u03b3=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 \", getDual(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", + " @addLazyConstraint(cb, t+sum{cv[i]*x[i], i in 1:dimX} <= \u03b3 )\n", + " println(\"t + \", cv, \"\u1d40 x <= \", \u03b3)\n", + " push!(stringOfBenderCuts, string(\"t+\", cv, \"'x <=\", \u03b3))\n", + " end\n", + " \n", + " if statusSubProblem == :Unbounded \n", + " println(\"\\nThere is an extreme ray, adding the corresponding constraint\")\n", + " ce = A1'*uCurrent\n", + " @addLazyConstraint(cb, sum{ce[i]*x[i], i in 1:dimX} <= \u03b3)\n", + " println(ce, \"x <= \", \u03b3)\n", + " push!(stringOfBenderCuts, string(ce, \"\u1d40 x <= \", \u03b3))\n", + " end\n", + " println(\"\\n\") \n", + " #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n", + " \n", + "end" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 14, + "text": [ + "addBendersCut (generic function with 1 method)" + ] + } + ], + "prompt_number": 14 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we tell the solver to use the callback function and solve the problem." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "setLazyCallback(masterProblemModel, addBendersCut) # Telling the solver to use the \n", + "# callback function" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 15, + "text": [ + "(addBendersCut,false)" + ] + } + ], + "prompt_number": 15 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "solve(masterProblemModel)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\r\n", + "User cuts applied: 2\r\n", + "\r\n", + "Root node processing (before b&c):\r\n", + " Real time = 0.27 sec. (0.04 ticks)\r\n", + "Sequential b&c:\r\n", + " Real time = 0.00 sec. (0.00 ticks)\r\n", + " ------------\r\n", + "Total (root+branch&cut) = 0.27 sec. (0.04 ticks)\r\n", + "Warning: Control callbacks may disable some MIP features.\r\n", + "Tried aggregator 1 time.\r\n", + "Reduced MIP has 0 rows, 3 columns, and 0 nonzeros.\r\n", + "Reduced MIP has 0 binaries, 2 generals, 0 SOSs, and 0 indicators.\r\n", + "Presolve time = 0.00 sec. (0.00 ticks)\r\n", + "Tried aggregator 1 time.\r\n", + "Reduced MIP has 0 rows, 3 columns, and 0 nonzeros.\r\n", + "Reduced MIP has 0 binaries, 2 generals, 0 SOSs, and 0 indicators.\r\n", + "Presolve time = 0.00 sec. (0.00 ticks)\r\n", + "MIP emphasis: balance optimality and feasibility.\r\n", + "MIP search method: traditional branch-and-cut.\r\n", + "Parallel mode: none, using 1 thread.\r\n", + "Root relaxation solution time = 0.00 sec. (0.00 ticks)\r\n", + "----------------------------" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n", + "ITERATION NUMBER = 1\n", + "---------------------------\n", + "\n", + "MASTERPROBLEM INFORMATION\n", + "-------------------------\n", + "The master problem that was solved was:\n", + "Max t\n", + "Subject to\n", + " x[i] >= 0, integer, for all i in {1,2}\n", + " t free\n", + "with 0 added lazy constraints\n", + "String[]\n", + "Current Value of x is: [0.0,0.0]\n", + "Current objective value of master problem, fmCurrent is: 1.0e20\n", + "\n", + "\n" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "The subproblem is being solved\n", + "Tried aggregator 1 time.\r\n", + "LP Presolve eliminated 1 rows and 1 columns.\r\n", + "Aggregator did 1 substitutions.\r\n", + "All rows and columns eliminated.\r\n", + "Presolve time = 0.00 sec. (0.00 ticks)\r\n", + "SUBPROBLEM INFORMATION\n", + "----------------------\n", + "The subproblem that was solved was: \n" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Min -2 u[1] - 3 u[2]\n", + "Subject to\n", + " u[1] - u[2] >= -2\n", + " -2 u[1] - u[2] >= -3\n", + " u[i] >= 0 for all i in {1,2}\n", + "Current status of the subproblem is Optimal\n", + "Current Value of u is: [0.3333333333333335,2.333333333333333]\n", + "Current Value of fs(xCurrent) is: -7.666666666666666\n", + "\n", + "\n", + "-------------------ADDING LAZY CONSTRAINT----------------\n", + "\n", + "There is a suboptimal vertex, add the corresponding constraint\n", + "t + [-0.9999999999999996,-4.0]\u1d40 x <= -7.666666666666666\n", + "\n", + "\n", + "----------------------------\n", + "ITERATION NUMBER = 2\n", + "---------------------------\n", + "\n", + "MASTERPROBLEM INFORMATION\n", + "-------------------------\n", + "The master problem that was solved was:\n", + "Max t\n", + "Subject to\n", + " x[i] >= 0, integer, for all i in {1,2}\n", + " t free\n", + "with 1 added lazy constraints\n", + "String[\"t+[-0.9999999999999996,-4.0]'x <=-7.666666666666666\"]\n", + "Current Value of x is: [0.0,2.5e19]\n", + "Current objective value of master problem, fmCurrent is: 1.0e20\n", + "\n", + "\n", + "The subproblem is being solved\n", + "Tried aggregator 1 time.\r\n", + "LP Presolve eliminated 2 rows and 2 columns.\r\n", + "All rows and columns eliminated.\r\n", + "Presolve time = 0.00 sec. (0.00 ticks)\r\n", + "SUBPROBLEM INFORMATION\n", + "----------------------\n", + "The subproblem that was solved was: \n", + "Min 7.5e19 u[1] + 7.5e19 u[2] - 1.0e20\n", + "Subject to\n", + " u[1] - u[2] >= -2\n", + " -2 u[1] - u[2] >= -3\n", + " u[i] >= 0 for all i in {1,2}\n", + "Current status of the subproblem is Optimal\n", + "Current Value of u is: [0.0,0.0]\n", + "Current Value of fs(xCurrent) is: -1.0e20\n", + "\n", + "\n", + "-------------------ADDING LAZY CONSTRAINT----------------\n", + "\n", + "There is a suboptimal vertex, add the corresponding constraint\n", + "t + [1.0,4.0]\u1d40 x <= -0.0\n", + "\n", + "\n", + "\r\n", + " Nodes Cuts/\r\n", + " Node Left Objective IInf Best Integer Best Bound ItCnt Gap\r\n", + "\r\n", + " 0 0 -3.8333 1 -3.8333 0 \r\n", + "----------------------------\n", + "ITERATION NUMBER = 3\n", + "---------------------------\n", + "\n", + "MASTERPROBLEM INFORMATION\n", + "-------------------------\n", + "The master problem that was solved was:\n", + "Max t\n", + "Subject to\n", + " x[i] >= 0, integer, for all i in {1,2}\n", + " t free\n", + "with 2 added lazy constraints\n", + "String[\"t+[-0.9999999999999996,-4.0]'x <=-7.666666666666666\",\"t+[1.0,4.0]'x <=-0.0\"]\n", + "Current Value of x is: [0.0,1.0]\n", + "Current objective value of master problem, fmCurrent is: -4.0\n", + "\n", + "\n", + "The subproblem is being solved\n", + "Tried aggregator 1 time.\r\n", + "LP Presolve eliminated 2 rows and 2 columns.\r\n", + "All rows and columns eliminated.\r\n", + "Presolve time = -0.00 sec. (0.00 ticks)\r\n", + "SUBPROBLEM INFORMATION\n", + "----------------------\n", + "The subproblem that was solved was: \n", + "Min u[1] - 4\n", + "Subject to\n", + " u[1] - u[2] >= -2\n", + " -2 u[1] - u[2] >= -3\n", + " u[i] >= 0 for all i in {1,2}\n", + "Current status of the subproblem is Optimal\n", + "Current Value of u is: [0.0,0.0]\n", + "Current Value of fs(xCurrent) is: -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 " + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "constrRefSubProblem: 1 dimensions:\n", + "[1] = -0.0\n", + "[2] = -0.0\n", + "\n", + "\n", + "\n" + ] + }, + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 16, + "text": [ + ":Optimal" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "* 0+ 0 -4.0000 -3.8333 0 4.17%\r\n", + " 0 1 -3.8333 1 -4.0000 -3.8333 0 4.17%\r\n", + "Elapsed time = 1.30 sec. (0.04 ticks, tree = 0.00 MB, solutions = 1)\r\n" + ] + } + ], + "prompt_number": 16 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [], + "language": "python", + "metadata": {}, + "outputs": [] + } + ], + "metadata": {} + } + ] +} \ No newline at end of file diff --git a/notebooks/Shuvomoy - Column generation.ipynb b/notebooks/Shuvomoy - Column generation.ipynb new file mode 100644 index 0000000..77789dd --- /dev/null +++ b/notebooks/Shuvomoy - Column generation.ipynb @@ -0,0 +1,938 @@ +{ + "metadata": { + "language": "Julia", + "name": "", + "signature": "sha256:700779758ccb2f45ddf8010aad45e697c4b31201d020bfd562bb84c3df669a5c" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Description:** This notebook describes how to implement column generation, which is a large scale optimization scheme, in JuMP. The cutting stock problem has been used as an illustrative example.\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 - column generation \n", + "--------------------------\n", + "
\n", + "\n", + "Implementing large scale optimization techniques such as column generation is really easy using JuMP. To explain how to implement column generation in JuMP, we consider the famous cutting stock problem. For more details about the problem, see pages 234-236 of Introduction to Linear Optimization by Bertsimas and Tsitsiklis." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Notation and notions:\n", + "\n", + "- Width of a large roll is $W$, and it needs to be cut into smaller width papers according to customer demand\n", + "- The set of indices of all feasible patterns is, $\\mathcal{J}=\\{1,2,\\ldots,n\\}$, where $n$ is a very large number\n", + "- A strict subset of $\\mathcal{J}$ that is considered in the master problem is $\\mathcal{J}'$\n", + "- The dummy index for a pattern is $j$ \n", + "- The index set of all possible paper-widths is, $\\mathcal{M}=\\{1,2,\\ldots,m\\}$\n", + "- The width of the paper with index $i$ is $w_i$\n", + "- The demand for the paper of width $w_i$ is $b_i$\n", + "- Number of smaller rolls of width $w_i$ produced by pattern $j$ is denoted by $a_{ij}$\n", + "- Number of large rolls cut according to pattern $j$ is denoted by $x_j$\n", + "\n", + "----------------------\n", + "## Original unabridged problem: \n", + "\n", + "\\begin{align}\n", + "&\\text{minimize} && \\sum_{j \\in \\mathcal{J}}{x_j} \\\\\n", + "&\\text{subject to} &&\\\\\n", + "& &&\\forall i \\in \\mathcal{M} \\quad \\sum_{j \\in \\mathcal{J}}{a_{ij} x_j}=b_i \\\\\n", + "& && \\forall j \\in \\mathcal{J} \\quad x_j \\geq 0 \\\\\n", + "\\end{align}\n", + "\n", + "----------------------\n", + "\n", + "## Structure of the decomposition\n", + "\n", + "**Master Problem:**\n", + "\\begin{align}\n", + "&\\text{minimize} && \\sum_{j \\in \\mathcal{J}'}{x_j} \\\\\n", + "&\\text{subject to} &&\\\\\n", + "& &&\\forall i \\in \\mathcal{M} \\quad \\sum_{j \\in \\mathcal{J}'}{a_{ij} x_j}=b_i \\\\\n", + "& && \\forall j \\in \\mathcal{J}' \\quad x_j \\geq 0 \\\\\n", + "\\end{align}\n", + "\n", + "**Subproblem:**\n", + "\\begin{align}\n", + "&\\text{minimize} && 1 - \\sum_{i \\in \\mathcal{M}} \\quad {p_i a_{i {j^*}}} \\\\\n", + "&\\text{subject to} &&\\\\\n", + "& && \\forall i \\in \\mathcal{M} \\quad a_{i {j^*}} \\geq 0, \\quad a_{ij^*} \\; \\text{integer} \\\\\n", + "& && \\sum_{i \\in \\mathcal{M}}{w_i a_{i{j^*}}} \\leq W\\\\\n", + "\\end{align}\n", + "\n", + "If objective value of the subproblem (min of the reduced cost vector) is greater than or equal to $0$, then it is optimal. Otherwise, add the column $(a_{i {j^*}})_{i \\in \\mathcal{M}}=A_{j*}$ and a corresponding new variable $x_{j*}$ is added to the master problem. The modified master problem is as follows:\n", + "\n", + "**Modified Master Problem**\n", + " \\begin{align}\n", + "&\\text{minimize} && \\sum_{j \\in \\mathcal{J}'}{x_j} + x_{j^*} \\\\\n", + "&\\text{subject to} &&\\\\\n", + "& &&\\forall i \\in \\mathcal{M} \\quad \\sum_{j \\in \\mathcal{J}'}{a_{ij} x_j}+a_{i j^*} x_{j^*}=b_i \\\\\n", + "& && \\forall j \\in \\mathcal{J}' \\quad x_j \\geq 0, x_j^* \\geq 0 \\\\\n", + "\\end{align}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Pseduocode\n", + "\n", + "- Input preliminary data for starting the problem\n", + "- Solve the master problem with the initial data\n", + "\\begin{align}\n", + "&\\text{minimize} && \\sum_{j \\in \\mathcal{J}'}{x_j} \\\\\n", + "&\\text{subject to} &&\\\\\n", + "& &&\\forall i \\in \\mathcal{M} \\quad \\sum_{j \\in \\mathcal{J}'}{a_{ij} x_j}=b_i \\\\\n", + "& && \\forall j \\in \\mathcal{J}' \\quad x_j \\geq 0 \\\\\n", + "\\end{align}\n", + "\n", + "- Collect the dual variables for the equality constraints and store them in an array $(p_i)_{i \\in \\mathcal{M}}$\n", + "\n", + "- Solve the sub problem \n", + "\\begin{align}\n", + "&\\text{minimize} && 1 - \\sum_{i \\in \\mathcal{M}} \\quad {p_i a_{i {j^*}}} \\\\\n", + "&\\text{subject to} &&\\\\\n", + "& && \\forall i \\in \\mathcal{M} \\quad a_{i {j^*}} \\geq 0, \\quad a_{ij^*} \\; \\text{integer} \\\\\n", + "& && \\sum_{i \\in \\mathcal{M}}{w_i a_{i{j^*}}} \\leq W\\\\\n", + "\\end{align}\n", + "- Flow control:
\n", + "\n", + "\n", + "while ( $\\text{optimal value of the subproblem} < 0$)
\n", + "> * Add the column $(a_{i {j^*}})_{i \\in \\mathcal{M}}=A_{j*}$ to $A$
\n", + "> * Add a corresponding new variable $x_{j*}$ to the list of variables
\n", + "> * Solve the modified master problem
\n", + " \\begin{align}\n", + "&\\text{minimize} && \\sum_{j \\in \\mathcal{J}'}{x_j} + x_{j^*} \\\\\n", + "&\\text{subject to} &&\\\\\n", + "& &&\\forall i \\in \\mathcal{M} \\quad \\sum_{j \\in \\mathcal{J}'}{a_{ij} x_j}+a_{i j^*} x_{j^*}=b_i \\\\\n", + "& && \\forall j \\in \\mathcal{J}' \\quad x_j \\geq 0 \\\\\n", + "& && \\qquad \\qquad \\; \\; x_{j^*} \\geq 0\n", + "\\end{align}\n", + "> * Collect the dual variables for the equality constraints and store them in an array $(p_i)_{i \\in \\mathcal{M}}$\n", + "> * Solve the sub problem as before
\n", + "> * Set $\\mathcal{J}':=\\mathcal{J}'\\cup \\{j^*\\}$
\n", + "\n", + " end while
\n", + "\n", + "\n", + "\n", + "- Display the results " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Master Problem Modification in JuMP\n", + "The problem modification can be done by using the already mentioned `@defVar` macro:\n", + "\n", + "$\n", + "@\\texttt{defVar}(m, l \\leq x_\\text{new} \\leq u, \\texttt{Int}, \\texttt{objective} = c_\\text{new}, \\texttt{inconstraints} = \\text{arrayConstrrefs}, \\texttt{coefficients} = \\text{arrayCoefficients}) \n", + "$\n", + "Here: \n", + "\n", + "- The name of the original model is $m$.\n", + "- The new variable to be added is $x_\\text{new}$ with lower bound $l$ and upper bound $u$.\n", + "- The type of the variable can be `Int`, `Bin`. For real variable the third argument is left vacant.\n", + "- The original objective, say $f_o(x)$ will become $f_o(x) + c_\\text{new} x_\\text{new}$ after modification\n", + "- The array $\\texttt{arrayConstrrefs}$ contain references to those constraints that need to be modified by inclusion of $x_\\text{new}$\n", + "- The array $\\texttt{arrayCoefficients}$ contain the coefficients that have to multiplied with $x_\\text{new}$ and then added to the constraints referened by $\\texttt{arrayConstrrefs}$ in an orderly manner. For example, if the $i$th element of $\\texttt{arrayConstrrefs}$ refers to a constraint $a_i^T x \\lesseqgtr b_i$, then after invoking the command, the constraint is modified as:\n", + "$a_i^T x +\\texttt{arrayCoefficients}[i] x_\\text{new} \\lesseqgtr b_i$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Implementing one iteration of the column generation algorithm\n", + "\n", + "To understand how the column generation is working in Julia, we implement one iteration of the column generation algorithm manually. The entire code is presented in the next section.\n" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + " # Uploading the packages:\n", + " # -----------------------\n", + "\n", + "using JuMP \n", + "\n", + "# We will use default solvers" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 1 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + " # Input preliminary data for starting the problem\n", + " # -----------------------------------------------\n", + "\n", + "W=110\n", + "cardinalityM=5\n", + "M=[1:cardinalityM]\n", + "A=eye(cardinalityM)\n", + "p=zeros(5)\n", + "b=[48; 35; 24; 10; 8]\n", + "w=[20; 45; 50; 55; 75]" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 2, + "text": [ + "5-element Array{Int64,1}:\n", + " 20\n", + " 45\n", + " 50\n", + " 55\n", + " 75" + ] + } + ], + "prompt_number": 2 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# Description of the master problem with the initial data\n", + "#----------------------\n", + "\n", + "cutstockMain = Model() # Model for the master problem\n", + "Jprime=[1:size(A,2)] # Intial number of variables\n", + "@defVar(cutstockMain, 0 <= x[Jprime] <= 1000000) # Defining the variables\n", + "@setObjective(cutstockMain, Min, sum{1*x[j],j in Jprime}) # Setting the objective\n", + "@addConstraint(cutstockMain, consRef[i=1:cardinalityM], sum{A[i,j]*x[j], j in Jprime}==b[i]) \n", + "# Here the second argument consRef[i=1:cardinalityM] means that the i-th constraint a\u1d62\u1d40x = b\u1d62 has the corresponding constraint reference\n", + "# consRef[i]\n", + "print(cutstockMain)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Min x[1] + x[2] + x[3] + x[4] + x[5]\n", + "Subject to\n", + " x[1] == 48\n", + " x[2] == 35\n", + " x[3] == 24\n", + " x[4] == 10\n", + " x[5] == 8\n", + " 0 <= x[i] <= 1.0e6 for all i in {1,2..4,5}\n" + ] + } + ], + "prompt_number": 4 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# Solving the master problem with the initial data\n", + "# ------------------------------------------------\n", + "solve(cutstockMain)\n", + "println(\"Current solution of the master problem is \", getValue(x))\n", + "println(\"Current objective value of the master problem is \", getObjectiveValue(cutstockMain))" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Current solution of the master problem is " + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "x: 1 dimensions, 5 entries:\n", + " [1] = 48.0\n", + " [2] = 35.0\n", + " [3] = 24.0\n", + " [4] = 10.0\n", + " [5] = 8.0" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n", + "Current objective value of the master problem is " + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "125.0\n" + ] + } + ], + "prompt_number": 5 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "#Collect the dual variables for the equality constraints and store them in an array p\n", + "for i in M\n", + " p[i] = getDual(consRef[i]) # These p[i] are the input data for the subproblem\n", + "end \n", + "println(\"The array storing the dual variables is \", p)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "The array storing the dual variables is " + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "[" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "1.0,1.0,1.0,1.0,1.0]\n" + ] + } + ], + "prompt_number": 6 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# Describe the sub problem\n", + "# ------------------------\n", + "cutstockSub=Model() # Model for the subproblem\n", + "@defVar(cutstockSub, 0 <= Ajstar[M] <= 1000000, Int )\n", + "@setObjective(cutstockSub, Min, 1-sum{p[i]*Ajstar[i],i in M})\n", + "@addConstraint(cutstockSub, sum{w[i]*Ajstar[i], i in M} <= W)\n", + "print(cutstockSub)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Min -Ajstar[1] - Ajstar[2] - Ajstar[3] - Ajstar[4] - Ajstar[5] + 1\n", + "Subject to\n", + " 20 Ajstar[1] + 45 Ajstar[2] + 50 Ajstar[3] + 55 Ajstar[4] + 75 Ajstar[5] <= 110\n", + " 0 <= Ajstar[i] <= 1.0e6, integer, for all i in {1,2..4,5}\n" + ] + } + ], + "prompt_number": 7 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# Solve the sub problem\n", + "# ---------------------\n", + "solve(cutstockSub)\n", + "minreducedCost=getObjectiveValue(cutstockSub)\n", + "println(\"The minimum component of the reduced cost vector is \", minreducedCost)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "The minimum component of the reduced cost vector is " + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "-4.0\n" + ] + } + ], + "prompt_number": 8 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The minimum component of the reduced cost vector is negative, so we have a suboptimal solution." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "if minreducedCost >= 0\n", + " println(\"We are done, current solution of the master problem is optimal\")\n", + "else\n", + " println(\"We have a cost reducing column \", getValue(Ajstar))\n", + "end\n" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "We have a cost reducing column " + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Ajstar: 1 dimensions, 5 entries:\n", + " [1] = 5.0\n", + " [2] = 0.0\n", + " [3] = 0.0\n", + " [4] = 0.0\n", + " [5] = 0.0\n" + ] + } + ], + "prompt_number": 9 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "typeof(Ajstar)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 10, + "text": [ + "JuMPDict{Variable,1} (constructor with 1 method)" + ] + } + ], + "prompt_number": 10 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now `Ajstar` is of type JuMPDict. To use it in the modified master problem, we have to store values from `Ajstar` in a column vector." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "Anew=Float64[] # This Anew correspond to the newly added column to the A matrix\n", + "for (i in 1:cardinalityM)\n", + " push!(Anew, getValue(Ajstar)[i])\n", + "end" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 11 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "When we add the cost reducing column `Anew` to the original matrix `A`, it also gives rise to a new variable `xNew` corresponding to `Anew`. Now we want to keep track of the new variables that are added by the subproblem. We do this by declaring an array of `Variable`s named `xNewArray`, which will contain all such newly added variables in the process of column generation. " + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "xNewArray=Variable[] # The newly added variables by flow control will be\n", + "# pushed to the new array of variables xNewArray" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 12, + "text": [ + "0-element Array{Variable,1}" + ] + } + ], + "prompt_number": 12 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here we just illustrate one iteration of the while loop manually, because, for now, we are interested to understand how JuMP is managing the flow control and modifying the master problem and the sub problem. \n", + "\n", + "Let's modify the master problem by adding the new column `Anew` to the old `A` matrix. Note that we do not have to rewrite the entire model." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# Modify the master problem by adding the new column Anew to the old A matrix\n", + "@defVar(\n", + "cutstockMain, # Model to be modified\n", + "0 <= xNew <= 1000000, # New variable to be added\n", + "objective=1, # cost coefficient of new varaible in the objective\n", + "inconstraints=consRef, # constraints to be modified\n", + "coefficients=Anew # the coefficients of the variable in those constraints\n", + ") \n", + "\n", + "# The line above adds the column (a\u1d62\u2c7c*)\u1d62=A\u2c7c* to A
\n", + "# and add a corresponding new variable x\u2c7c* to the list of variable\n", + "\n", + "push!(xNewArray, xNew) # Pushing the new variable in the array of new variables\n", + "print(cutstockMain)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Min x[1] + x[2] + x[3] + x[4] + x[5] + xNew\n", + "Subject to\n", + " x[1] + 5 xNew == 48\n", + " x[2] == 35\n", + " x[3] == 24\n", + " x[4] == 10\n", + " x[5] == 8\n", + " 0 <= x[i] <= 1.0e6 for all i in {1,2..4,5}\n", + " 0 <= xNew <= 1.0e6\n" + ] + } + ], + "prompt_number": 13 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Though we are showing only one iteration of the flow control, in the final code for sure we want to have a \n", + ">```\n", + "while ( some condition )\n", + "(\n", + "...\n", + ")\n", + "end\n", + "```\n", + "\n", + "block. \n", + "\n", + "Now if we do not do anything else in the final code, all the names of the newly added variables by the `while` loop will be the same: `xNew`! JuMP is intelligent enough to treat them as seperate varaibles, but it is not very human-friendly. It is more convenient if the newly added variables were given different names, which we can achieve by `setName(oldName, newName)` function.\n" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "setName(xNew, string(\"x[\",size(A,2)+1,\"]\")) # Changing the name of the variable \n", + "# otherwise all the newly added variables will have name xNew!\n", + "# size(A,2) gives the column number of A" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 14, + "text": [ + "\"x[6]\"" + ] + } + ], + "prompt_number": 14 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let us see if the name of the variable has changed as desired." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "print(cutstockMain) # Let us see if the name of the variables have changed as desired" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Min x[1] + x[2] + x[3] + x[4] + x[5] + x[6]\n", + "Subject to\n", + " x[1] + 5 x[6] == 48\n", + " x[2] == 35\n", + " x[3] == 24\n", + " x[4] == 10\n", + " x[5] == 8\n", + " 0 <= x[i] <= 1.0e6 for all i in {1,2..4,5}\n", + " 0 <= x[6] <= 1.0e6\n" + ] + } + ], + "prompt_number": 15 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Indeed it has! Now let's solve the modified master problem, and then collect the associated dual variables for the equality constraints and store them in the array `p`." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "statusControlFlow=solve(cutstockMain) # Solve the modified master problem\n", + "\n", + "getDual(consRef)\n", + "for i in M\n", + " p[i] = getDual(consRef)[i] \n", + "end \n", + "\n", + "println(p)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "[" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "0.2,1.0,1.0,1.0,1.0]\n" + ] + } + ], + "prompt_number": 16 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we solve the subproblem for the current solution of the master problem:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# Solving the modified sub problem \n", + "@defVar(cutstockSub, 0 <= Ajstar[M] <= 1000000, Int )\n", + "@setObjective(cutstockSub, Min, 1-sum{p[i]*Ajstar[i],i in M})\n", + "@addConstraint(cutstockSub, sum{w[i]*Ajstar[i], i in M} <= W)\n", + "print(cutstockSub) # Let's see what is the current subproblem looks like\n", + "solve(cutstockSub)\n", + "minReducedCost=getObjectiveValue(cutstockSub)\n", + "println(\"Current value of the minimum of the reduced cost vector is \", minReducedCost)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Min -0.2 Ajstar[1] - Ajstar[2] - Ajstar[3] - Ajstar[4] - Ajstar[5] + 1\n", + "Subject to\n", + " 20 Ajstar[1] + 45 Ajstar[2] + 50 Ajstar[3] + 55 Ajstar[4] + 75 Ajstar[5] <= 110\n", + " 20 Ajstar[1] + 45 Ajstar[2] + 50 Ajstar[3] + 55 Ajstar[4] + 75 Ajstar[5] <= 110\n", + " 0 <= Ajstar[i] <= 1.0e6, integer, for all i in {1,2..4,5}\n", + " 0 <= Ajstar[i] <= 1.0e6, integer, for all i in {1,2..4,5}\n" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Current value of the minimum of the reduced cost vector is " + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "-1.2000000000000002\n" + ] + } + ], + "prompt_number": 17 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The optimal value of the current subproblem is negative (which will be tested by the conditional statement of the while loop in the final code), giving us a cost reducing column to be added in the master problem. As before we have to store the column `Ajstar` in a column vector `Anew`." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "#Store the components of the solution of current subproblem into the column Anew \n", + "Anew=Float64[]\n", + "for (i in 1:cardinalityM)\n", + " push!(Anew, getValue(Ajstar)[i])\n", + "end\n", + "\n", + "println(\"New column to be added to A is: \", Anew)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "New column to be added to A is: " + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "[1.0,2.0,0.0,0.0,0.0]\n" + ] + } + ], + "prompt_number": 18 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Okay, we have understood how JuMP is working in the column generation process. The entire code of the cutting stock problem is given below:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Cutting stock problem code:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + " # Verfied to be working:\n", + "\n", + " # Uploading the packages:\n", + " # -----------------------\n", + "\n", + "using JuMP\n", + "using GLPKMathProgInterface\n", + "\n", + " # Input preliminary data for starting the problem\n", + " # ----------------------\n", + "\n", + "W=110\n", + "cardinalityM=5\n", + "M=[1:cardinalityM]\n", + "A=eye(cardinalityM)\n", + "p=zeros(5)\n", + "b=[48; 35; 24; 10; 8]\n", + "w=[20; 45; 50; 55; 75]\n", + "\n", + "@time begin # time measurement begins\n", + "\n", + " # Solve the master problem with the initial data\n", + " #----------------------\n", + "\n", + "cutstockMain = Model() # Model for the master problem\n", + "Jprime=[1:size(A,2)] # Intial number of variables\n", + "@defVar(cutstockMain, 0 <= x[Jprime] <= 1000000) # Defining the variables\n", + " \n", + "@setObjective(cutstockMain, Min, sum{1*x[j],j in Jprime}) # Setting the objective\n", + " \n", + "@addConstraint(cutstockMain, consRef[i=1:cardinalityM], sum{A[i,j]*x[j], j in Jprime}==b[i]) # Adding the constraints\n", + "# Here the second argument consRef[i=1:cardinalityM] means that the i-th constraint a\u1d62\u1d40x = b\u1d62 has \n", + "# the corresponding constraint reference consRef[i]\n", + "\n", + "solve(cutstockMain)\n", + " \n", + "#Collect the dual variables for the equality constraints and store them in an array p\n", + "getDual(consRef)\n", + "for i in M\n", + " p[i] = getDual(consRef)[i] # These p[i] are the input data for the subproblem\n", + "end \n", + "\n", + " # Solve the sub problem\n", + " # -------------------\n", + "\n", + "cutstockSub=Model() # Model for the subproblem\n", + "@defVar(cutstockSub, 0 <= Ajstar[M] <= 1000000, Int )\n", + "@setObjective(cutstockSub, Min, 1-sum{p[i]*Ajstar[i],i in M})\n", + "@addConstraint(cutstockSub, sum{w[i]*Ajstar[i], i in M} <= W)\n", + "solve(cutstockSub)\n", + "minReducedCost=getObjectiveValue(cutstockSub)\n", + "\n", + "Anew=Float64[] # This Anew correspond to the newly added column to the A matrix\n", + "for (i in 1:cardinalityM)\n", + " push!(Anew, getValue(Ajstar)[i])\n", + "end\n", + "\n", + "xNewArray=Variable[] # The newly added variables by flow control will be pushed to the new array of variables xNewArray\n", + "\n", + "k=1 # Counter for the while loop\n", + "\n", + " # Flow control\n", + " # ------------\n", + "\n", + "while (minReducedCost < 0) #while (current solution of the master problem is suboptimal, i.e., subproblem objective value < 0)\n", + " # Solve the master problem by adding the new column Anew to the old A matrix\n", + " @defVar(\n", + " cutstockMain, # Model to be modified\n", + " 0 <= xNew <= 1000000, # New variable to be added\n", + " objective=1, # cost coefficient of new varaible in the objective\n", + " inconstraints=consRef, # constraints to be modified\n", + " coefficients=Anew # the coefficients of the variable in those constraints\n", + " ) \n", + " # The line above adds the column (a\u1d62\u2c7c*)\u1d62=A\u2c7c* to A
\n", + " # and add a corresponding new variable x\u2c7c* to the list of variable\n", + " push!(xNewArray, xNew) # Pushing the new variable in the array of new variables\n", + " setName(xNew, string(\"x[\",size(A,2)+k,\"]\")) # Changing the name of the variable \n", + " # otherwise all the newly added variables will have name xNew!\n", + " k=k+1 # Increasing k by 1\n", + "statusControlFlow=solve(cutstockMain)\n", + "\n", + "#Collect the dual variables for the equality constraints and store them in an array p\n", + "getDual(consRef)\n", + "for i in M\n", + " p[i] = getDual(consRef)[i] \n", + "end \n", + " \n", + "# Solving the modified sub problem \n", + "@defVar(cutstockSub, 0 <= Ajstar[M] <= 1000000, Int )\n", + "@setObjective(cutstockSub, Min, 1-sum{p[i]*Ajstar[i],i in M})\n", + "@addConstraint(cutstockSub, sum{w[i]*Ajstar[i], i in M} <= W)\n", + "solve(cutstockSub)\n", + "minReducedCost=getObjectiveValue(cutstockSub)\n", + "\n", + "#Store the components of the solution of current subproblem into the column Anew \n", + "Anew=Float64[]\n", + "for (i in 1:cardinalityM)\n", + " push!(Anew, getValue(Ajstar)[i])\n", + "end\n", + " end # While loop ends\n", + " \n", + " end # time measurement ends\n", + "\n", + " # Print the results\n", + " # -----------------\n", + "\n", + "println(\"Objective value: \", getObjectiveValue(cutstockMain))\n", + "println(\"Current Solution is: \", getValue(x))\n", + "println(\"With \", length(xNewArray), \" variables added by flow control:\")\n", + "for i in 1:length(xNewArray)\n", + " println(\"[\",size(A,2)+i,\"] = \",getValue(xNewArray[i]))\n", + "end\n", + "println(\"Reduced cost of the current solution is \", getObjectiveValue(cutstockSub))" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "elapsed time: " + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "0.147638368 seconds (799460 bytes allocated)\n", + "Objective value: 46.25\n", + "Current Solution is: x: 1 dimensions, 5 entries:\n", + " [1] = 0.0\n", + " [2] = 0.0\n", + " [3] = 0.0\n", + " [4] = 0.0\n", + " [5] = 0.0\n", + "With 6 variables added by flow control:\n" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "[6] = 0.0\n", + "[7] = 17.5\n", + "[8] = 8.25\n", + "[9] = 5.0\n", + "[10] = 8.0\n", + "[11] = 7.499999999999999\n", + "Reduced cost of the current solution is 0.0\n" + ] + } + ], + "prompt_number": 20 + } + ], + "metadata": {} + } + ] +} \ No newline at end of file diff --git a/notebooks/Shuvomoy - Exploiting sparsity.ipynb b/notebooks/Shuvomoy - Exploiting sparsity.ipynb new file mode 100644 index 0000000..5570c83 --- /dev/null +++ b/notebooks/Shuvomoy - Exploiting sparsity.ipynb @@ -0,0 +1,1040 @@ +{ + "metadata": { + "language": "Julia", + "name": "", + "signature": "sha256:47332028aaa24f65ac4a8b4dc79424584346d0dd7d47ee1e1c587b99f84d5402" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Description:** This notebook describes how to make an optimization model more efficient by exploiting sparsity.\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 - exploiting sparsity \n", + "--------------------------\n", + "
\n", + "\n", + "## Introduction\n", + "\n", + "**What is a sparse data structure?**\n", + "\n", + "A sparse data structure is one that has a lot of zeros in it. If a matrix has many more zeros than nonzeros, then it is a sparse matrix. \n", + "\n", + "** Why do we need to exploit sparsity?** \n", + "\n", + "- Sparsity in the input data increases with the dimension.\n", + "\n", + "- Exploiting sparsity \n", + " * keeps the data size small\n", + " * saves memory\n", + " * reduces the running time \n", + " * improves the efficiency of the model\n", + "\n", + "\n", + "** How to exploit sparsity in Julia?**\n", + "\n", + "- Define `immutable type` and create a `dictionary` or an `array` of it." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## A test example: Transportation Problem\n", + "Consider a transportation problem which is going to be our test example:\n", + "\n", + "- ** Problem setup:** Some products have to transported from origin cities to destination cities \n", + "\n", + "- **Objective:** \n", + " * Minimize the total cost of shipment over all **relevant** routes\n", + " \n", + "- **Decision Variables** \n", + " * Find the optimum amount of every product to be shipped from one city to another\n", + "\n", + "- **Constraints**\n", + " * How much of a product a city can supply to other cities is fixed. \n", + " * The amount of any product demanded by a city is also fixed. \n", + " * The total amount of products shipped between every pair of different cities can not exceed a given limit. \n", + "\n", + "Suppose there are ten cities and three products in our problem. " + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "cities = \n", + "[\n", + ":GARY; :CLEV; :PITT; :FRA; :DET; :LAN; :WIN; :STL; :FRE; :LAF \n", + "]\n", + "\n", + "products = \n", + "[\n", + ":bands; :coils; :plate \n", + "]" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 1, + "text": [ + "3-element Array{Symbol,1}:\n", + " :bands\n", + " :coils\n", + " :plate" + ] + } + ], + "prompt_number": 1 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Defining `immutable` type\n", + "If we do not exploit sparsity, then number of ways we can ship the products from one city to other will be $3 \\times (^{10}P_2+10)=300$. \n", + "\n", + "Clearly many of them will be redundant, because of reasons like\n", + "- a product might not be needed by a city\n", + "- a product might not be produced by a city\n", + "etc.\n", + "\n", + "We just need to consider *relevant route*s, where a product can be shipped from one production city to the other demand city. So a *relevant route* can be defined by 3 features:\n", + "\n", + "- a product\n", + "- a city that produces that product\n", + "- a city that demands that product\n", + "\n", + "So, we define an `immutable` type `Route` as follows:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "immutable Route\n", + " p::Symbol # p stands for product\n", + " o::Symbol # o stands for origin\n", + " d::Symbol # d stands for destination\n", + "end" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 2 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here the datatype `Symbol` is a special type of immutable string. Then we create an array of only relevant routes." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "routesExample = \n", + "[\n", + " Route(:bands,:GARY,:FRA);\n", + " Route(:bands,:GARY,:DET);\n", + " Route(:bands,:GARY,:LAN);\n", + " Route(:bands,:GARY,:WIN);\n", + "] " + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 4, + "text": [ + "4-element Array{Route,1}:\n", + " Route(:bands,:GARY,:FRA)\n", + " Route(:bands,:GARY,:DET)\n", + " Route(:bands,:GARY,:LAN)\n", + " Route(:bands,:GARY,:WIN)" + ] + } + ], + "prompt_number": 4 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If we want to access $i$th element of the array by typing in `routes[i]`. When we want to access the product name associated with the $it$h elelment of the array, we can do so by typing `routes[i].p`.\n" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "routesExample[2] # Will give the second route " + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 14, + "text": [ + "Route(:bands,:GARY,:DET)" + ] + } + ], + "prompt_number": 14 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "routesExample[4].d # Will give the demand city of the 4th route" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 5, + "text": [ + ":WIN" + ] + } + ], + "prompt_number": 5 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Creating new arrays efficiently from existing arrays\n", + "\n", + "Often we need to create new arrays, where the elements of them are extracted from some already existing array conditionally. \n", + "\n", + "Consider the immutable type `Supply`" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "immutable Supply\n", + " p::Symbol # p stands for product name\n", + " o::Symbol # o stands for the origin city\n", + "end" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 8 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We want to create an array `suppliesExample`, that contains all relevant product-city pairs, where the particular product is produced in that city. Clearly we can construct this array by plucking each product and corresponding city producing it from `routesExample`. This is how we do it efficiently: \n", + "\n", + "- Create a emty array of type `Supply`\n", + "- Add elements to this array by \n", + "\n", + " * selecting the product and origin from the elements of `routes` \n", + " * *push*ing them one by one in `supplies`" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "suppliesExample = Supply[] # Creates a 0 element array of immutable type Supply\n", + "for r in routesExample # For every element of the route route\n", + " push!(suppliesExample, Supply(r.p, r.o)) # pick the product and origin city and push it in supplies\n", + "end" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 9 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "suppliesExample" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 10, + "text": [ + "4-element Array{Supply,1}:\n", + " Supply(:bands,:GARY)\n", + " Supply(:bands,:GARY)\n", + " Supply(:bands,:GARY)\n", + " Supply(:bands,:GARY)" + ] + } + ], + "prompt_number": 10 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Dictionary\n", + "**What is a dictionary?**\n", + "A dictionary is a datatype which can be useful in exploiting sparsity.\n", + "\n", + "**Why is it needed?**\n", + " Often we might be interested to index a variable by a composite data type, rather than a number. \n", + " \n", + "For example, for the transportation problem in consideration, it would be more convenient to index the decision variables in the routes that are present. Let\n", + "\n", + "\\begin{align}\n", + "R=\\{(p,o,d) \\in P \\times C \\times C: \\text{product } p \\text{ has to be transported from city } o \\text{ to city } d\\}\n", + "\\end{align}\n", + "\n", + "be the set of all the routes that are relevant for the problem. So, we can define our decision variables as below:\n", + "\n", + "\\begin{align}\n", + "\\forall (p,o,d) \\in R \\left(x_{(p,o,d)}= \\text{the amount of a product } p \\text{ that is trasported from city } o \\text{ city d} \\right)\n", + "\\end{align}\n", + "\n", + "From a data structure point of view, $ \\left(x_{(p,o,d)}\\right)_{(p,o,d) \\in R}$ is a dictionary which\n", + "\n", + "- takes $(p,o,d) \\in R$ as its **key** and \n", + "- has the **value** the optimum amount of the product $p$ to be shipped from city $o$ to city $d$." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Efficiently Constructing Dictionary of `immutable` types\n", + "\n", + "Suppose we want to create a dictionary called `costRoutes`. Every element of the dictionary `costRoutes` contains the value of shipping cost along a particular route belonging to the array `routesExample`. So, \n", + "\n", + "- the **key** to an element belonging to the dictionary is a specific route belonging to the array `routesExample`\n", + "- the **value** is the cost for that shipment. \n", + "\n", + "Suppose the values of the costs are stored in an array named `costCofExample`. " + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "costCofExample = [100; 200; 300; 50.0]" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 11, + "text": [ + "4-element Array{Float64,1}:\n", + " 100.0\n", + " 200.0\n", + " 300.0\n", + " 50.0" + ] + } + ], + "prompt_number": 11 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We create the dictionary `costRoutes` similar to an array: \n", + "\n", + "- we create an empty dictionary, and then\n", + "- push the elements of the dictionary (the key and the corresponding value) one by one.\n", + "\n", + "``\n", + "costRoutesExample=Dict{Route, Float64}()# Create an empty dictionary where the key is Route and the value is Float64\n", + "for i in [1:length(routes)]\n", + " push!(costRoutesExample, routes[i], costCof[i]) # routes[i] is the key, and costCof[i] is the value\n", + "end\n", + "``\n", + "\n" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "costRoutesExample=Dict{Route, Float64}()# Create an empty dictionary \n", + "# where the key is Route and the value is Float64" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 12, + "text": [ + "Dict{Route,Float64} with 0 entries" + ] + } + ], + "prompt_number": 12 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "for i in [1:length(routesExample)]\n", + " push!(costRoutesExample, routesExample[i], costCofExample[i]) \n", + " # routesExample[i] is the key, and costCofExample[i] is the value\n", + "end" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 13 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "costRoutesExample" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 14, + "text": [ + "Dict{Route,Float64} with 4 entries:\n", + " Route(:bands,:GARY,:WIN) => 50.0\n", + " Route(:bands,:GARY,:FRA) => 100.0\n", + " Route(:bands,:GARY,:LAN) => 300.0\n", + " Route(:bands,:GARY,:DET) => 200.0" + ] + } + ], + "prompt_number": 14 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "After the dictionary is initialized, we can access the cost associated with some route `routes[i]` by typing in `costRoutes[routes[i]]`" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "costRoutesExample[routesExample[4]]" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 15, + "text": [ + "50.0" + ] + } + ], + "prompt_number": 15 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Or we can input the description of the route itself:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "routesExample[3]" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 16, + "text": [ + "Route(:bands,:GARY,:LAN)" + ] + } + ], + "prompt_number": 16 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "costRoutesExample[Route(:bands,:GARY,:LAN)]" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 17, + "text": [ + "300.0" + ] + } + ], + "prompt_number": 17 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Mathematical representation of the transportation problem\n", + "\n", + "\n", + "The problem is a classic transportation problem. We will consider the sparse representation of the problem. Let\n", + "\\begin{align}\n", + "&C=\\text{Set of all cities} \\\\\n", + "&P=\\text{Set of all products} \\\\\n", + "&R=\\{(p,o,d) \\in P \\times C \\times C: \\text{product } p \\text{ has to be transported from city } o \\text{ to city } d\\} \\\\\n", + "&\\forall (p,o,d) \\in R \\quad \\left(c_{(p,o,d)} = \\text{cost of transporting some product } p \\text{ from city } o \\text{ to city } d \\right)\\\\\n", + "& \\forall p \\in P \\quad \\left(O_p = \\text{Set of all origin cities where a product } p \\text{ is produced}\\right) \\\\\n", + "& \\forall p \\in P \\quad \\left(D_p = \\text{Set of all destination cities where a product } p \\text{\n", + "has to be delivered}\\right) \\\\\n", + "& \\forall p \\in P \\; \\forall o \\in O_p \\quad \\left(s_{(p,o)} = \\text{the total amount of the product } p \\text{ that can be supplied by city } o\\right) \\\\\n", + "& \\forall p \\in P \\; \\forall d \\in D_p \\quad \\left(d_{(p,d)} = \\text{the total amount of the product } p \\text{ that is demanded by city } d\\right) \\\\\n", + "& \\text{Total amount of shipped product between each pair of cities cannot exceed } \\gamma\n", + "\\end{align}\n", + "\n", + "The decision variable for this problem is \n", + "\\begin{align}\n", + "\\forall (p,o,d) \\in R \\quad \\left(x_{(p,o,d)}= \\text{the amount of a product } p \\text{ that is trasported from city } o \\text{ city d} \\right)\n", + "\\end{align}\n", + "The optimization problem can be described as below:\n", + "\\begin{align}\n", + "&\\text{minimize} && \\sum_{(p,o,d) \\in R}{c_{(p,o,d)} x_{(p,o,d)}} \\\\\n", + "&\\text{subject to} && \\\\\n", + "& && \\forall p \\in P \\; \\forall o \\in O_p \\quad \\left(\\sum_{d \\in D_p} x_{(p,o,d)}=s_{(p,o)}\\right) \\\\\n", + "& && \\text{ : The amount of any product a city can supply to other cities is fixed} \\\\\n", + "& && \\forall p \\in P \\; \\forall d \\in D_p \\quad \\left(\\sum_{o \\in O_p}{x_{(p,o,d)}}=d_{(p,d)}\\right) \\\\\n", + "& && \\text{ : The amount of any product demanded by a city is also fixed.} \\\\\n", + "& && \\forall o \\in C \\; \\forall d \\in C \\setminus \\{o\\} \\quad \n", + " \\left(\\sum_{(p,\\bar{o},\\bar{d}) \\in R \\; : \\; o=\\bar{o} \\wedge d=\\bar{d}}{x_{(p,\\bar{o},\\bar{d})}}\\leq \\gamma\\right) \\\\\n", + "& && \\text{ :The total amount of products shipped between every pair of different cities can not exceed a given limit.}\n", + " \\end{align}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Mapping of the mathematical symbols to JuMP\n", + "\n", + "In the data file, the symbols in the model above are mapped as follows:\n", + "\n", + "| Mathematical Symbol | In the code | Comment |\n", + "|---------------------- :|------------: |-----------:|\n", + "| $C$ | `cities` | `cities` is an array of `Symbol`s |\n", + "| $P$ | `products` | `products` is an array of `Symbol`s |\n", + "| $R$ | `routes` | `routes` is an array of immutable type `Route` |\n", + "| $(O_p)_{p \\in P}$ | `orig` | `orig` is a dictionary |\n", + "| $(D_p)_{p \\in P}$ | `dest` | `dest` is a dictionary |\n", + "| $((s_{(p,o)})_{o \\in O_p})_{p \\in P}$ | `suppliedAmount`| `suppliedAmount` is a dictionary |\n", + "| $((d_{(p,d)})_{d \\in D_p})_{p \\in P}$ | `demandedAmount` | `demandedAmount` is a dictionary with key `Customer` and value `Float64` |\n", + "| $(x_{(p,o,d)})_{(p,o,d) \\in R}$ | `trans` | `trans` is a dictionary and the variable in the problem |\n", + "| $(c_{(p,o,d)})_{(p,o,d) \\in R}$ | `costRoutes` | `costRoutes` is a dictionary |\n", + "| $\\gamma$ | `capacity` | It is of type `Float64` |\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Data file" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "#\n", + "cities = \n", + "[\n", + ":GARY; :CLEV; :PITT; :FRA; :DET; :LAN; :WIN; :STL; :FRE; :LAF \n", + "]\n", + "\n", + "products = \n", + "[\n", + ":bands; :coils; :plate \n", + "]\n", + "\n", + "capacity = 625 \n", + "\n", + "immutable Route\n", + " p::Symbol # p stands for product\n", + " o::Symbol # o stands for origin\n", + " d::Symbol # d stands for destination\n", + "end\n", + "\n", + "routes = \n", + "[\n", + " Route(:bands,:GARY,:FRA);\n", + " Route(:bands,:GARY,:DET);\n", + " Route(:bands,:GARY,:LAN);\n", + " Route(:bands,:GARY,:WIN);\n", + " Route(:bands,:GARY,:STL);\n", + " Route(:bands,:GARY,:FRE);\n", + " Route(:bands,:GARY,:LAF);\n", + " Route(:bands,:CLEV,:FRA);\n", + " Route(:bands,:CLEV,:DET);\n", + " Route(:bands,:CLEV,:LAN);\n", + " Route(:bands,:CLEV,:WIN);\n", + " Route(:bands,:CLEV,:STL);\n", + " Route(:bands,:CLEV,:FRE);\n", + " Route(:bands,:CLEV,:LAF); \n", + " Route(:bands,:PITT,:FRA);\n", + " Route(:bands,:PITT,:DET);\n", + " Route(:bands,:PITT,:LAN);\n", + " Route(:bands,:PITT,:WIN);\n", + " Route(:bands,:PITT,:STL);\n", + " Route(:bands,:PITT,:FRE);\n", + " Route(:bands,:PITT,:LAF); \n", + "\n", + " Route(:coils,:GARY,:FRA);\n", + " Route(:coils,:GARY,:DET);\n", + " Route(:coils,:GARY,:LAN);\n", + " Route(:coils,:GARY,:WIN);\n", + " Route(:coils,:GARY,:STL);\n", + " Route(:coils,:GARY,:FRE);\n", + " Route(:coils,:GARY,:LAF); \n", + "\n", + " Route(:coils,:CLEV,:FRA);\n", + " Route(:coils,:CLEV,:DET);\n", + " Route(:coils,:CLEV,:LAN);\n", + " Route(:coils,:CLEV,:WIN);\n", + " Route(:coils,:CLEV,:STL);\n", + " Route(:coils,:CLEV,:FRE);\n", + " Route(:coils,:CLEV,:LAF); \n", + "\n", + " Route(:coils,:PITT,:FRA);\n", + " Route(:coils,:PITT,:DET);\n", + " Route(:coils,:PITT,:LAN);\n", + " Route(:coils,:PITT,:WIN);\n", + " Route(:coils,:PITT,:STL);\n", + " Route(:coils,:PITT,:FRE);\n", + " Route(:coils,:PITT,:LAF); \n", + "\n", + " Route(:plate,:GARY,:FRA);\n", + " Route(:plate,:GARY,:DET);\n", + " Route(:plate,:GARY,:LAN);\n", + " Route(:plate,:GARY,:WIN);\n", + " Route(:plate,:GARY,:STL);\n", + " Route(:plate,:GARY,:FRE);\n", + " Route(:plate,:GARY,:LAF); \n", + "\n", + " Route(:plate,:CLEV,:FRA);\n", + " Route(:plate,:CLEV,:DET);\n", + " Route(:plate,:CLEV,:LAN);\n", + " Route(:plate,:CLEV,:WIN);\n", + " Route(:plate,:CLEV,:STL);\n", + " Route(:plate,:CLEV,:FRE);\n", + " Route(:plate,:CLEV,:LAF); \n", + "\n", + " Route(:plate,:PITT,:FRA);\n", + " Route(:plate,:PITT,:DET);\n", + " Route(:plate,:PITT,:LAN);\n", + " Route(:plate,:PITT,:WIN);\n", + " Route(:plate,:PITT,:STL);\n", + " Route(:plate,:PITT,:FRE);\n", + " Route(:plate,:PITT,:LAF); \n", + "]\n", + "\n", + "immutable Supply\n", + " p::Symbol\n", + " o::Symbol\n", + "end\n", + "\n", + "#Creating the array supplies\n", + "# --------\n", + "supplies = Supply[] # Creates a 0 element array of immutable type Supply\n", + "for r in routes\n", + " push!(supplies, Supply(r.p, r.o))\n", + "end\n", + "\n", + "# Creating suppliedAmount dictionary\n", + "# --------------\n", + "#It might be better to create this as a dictionary, where the key is the \n", + "# element of the array supplies and the value is the corresponding supplied\n", + "#amount\n", + "suppliedAmount = Dict{Supply, Float64}()\n", + "for s in supplies\n", + " if s.p == :bands && s.o == :CLEV\n", + " push!(suppliedAmount, s, 700)\n", + " elseif s.p == :bands && s.o==:GARY\n", + " push!(suppliedAmount, s, 400)\n", + " elseif s.p == :bands && s.o==:PITT\n", + " push!(suppliedAmount, s, 800)\n", + " elseif s.p == :coils && s.o==:GARY\n", + " push!(suppliedAmount, s, 800)\n", + " elseif s.p == :coils && s.o==:CLEV\n", + " push!(suppliedAmount, s, 1600)\n", + " elseif s.p == :coils && s.o == :PITT\n", + " push!(suppliedAmount, s, 1800)\n", + " elseif s.p == :plate && s.o == :GARY\n", + " push!(suppliedAmount, s, 200)\n", + " elseif s.p == :plate && s.o == :CLEV\n", + " push!(suppliedAmount, s, 300)\n", + " elseif s.p == :plate && s.o == :PITT\n", + " push!(suppliedAmount, s, 300)\n", + " else\n", + " push!(suppliedAmount, -1) \n", + "end #if\n", + "end #for\n", + "\n", + "immutable Customer\n", + " p::Symbol\n", + " d::Symbol\n", + "end\n", + "\n", + "# Creating customers array, which is an array of custom immutable Customer \n", + "# ---------\n", + "customers = Customer[]\n", + "for r in routes\n", + " push!(customers, Customer(r.p, r.d))\n", + "end\n", + "\n", + "# Let's create the demandedAmount dictionary now\n", + "# --------------\n", + "demandedAmount = Dict{Customer, Float64}()\n", + "for c in customers\n", + "#1\n", + " if c.p==:bands && c.d==:FRA\n", + " push!(demandedAmount, c, 300) \n", + "#2\n", + " elseif c.p==:coils && c.d==:FRA\n", + " push!(demandedAmount, c, 500) \n", + "#3\n", + " elseif c.p==:plate && c.d==:FRA\n", + " push!(demandedAmount, c, 100) \n", + "#4\n", + " elseif c.p==:bands && c.d==:DET\n", + " push!(demandedAmount, c, 300) \n", + "#5\n", + " elseif c.p==:coils && c.d==:DET\n", + " push!(demandedAmount, c, 750) \n", + "#6\n", + " elseif c.p==:plate && c.d==:DET \n", + " push!(demandedAmount, c, 100) \n", + "#7\n", + " elseif c.p==:bands && c.d==:LAN\n", + " push!(demandedAmount, c, 100) \n", + "#8\n", + " elseif c.p==:coils && c.d==:LAN\n", + " push!(demandedAmount, c, 400) \n", + "#9\n", + " elseif c.p==:plate && c.d==:LAN\n", + " push!(demandedAmount, c, 0) \n", + "#10\n", + " elseif c.p==:bands && c.d==:WIN\n", + " push!(demandedAmount, c, 75) \n", + "#11\n", + " elseif c.p==:coils && c.d==:WIN\n", + " push!(demandedAmount, c, 250) \n", + "#12\n", + " elseif c.p==:plate && c.d==:WIN\n", + " push!(demandedAmount, c, 50)\n", + "#13\n", + " elseif c.p==:bands && c.d==:STL\n", + " push!(demandedAmount, c, 650)\n", + "#14\n", + " elseif c.p==:coils && c.d==:STL\n", + " push!(demandedAmount, c, 950)\n", + "#15\n", + " elseif c.p==:plate && c.d==:STL\n", + " push!(demandedAmount, c, 200) \n", + "#16\n", + " elseif c.p==:bands && c.d==:FRE\n", + " push!(demandedAmount, c, 225)\n", + "#17\n", + " elseif c.p==:coils && c.d==:FRE\n", + " push!(demandedAmount, c, 850)\n", + "#18\n", + " elseif c.p==:plate && c.d==:FRE\n", + " push!(demandedAmount, c, 100)\n", + "#19\n", + " elseif c.p==:bands && c.d==:LAF\n", + " push!(demandedAmount, c, 250)\n", + "#20\n", + " elseif c.p==:coils && c.d==:LAF\n", + " push!(demandedAmount, c, 500)\n", + "#21 \n", + " elseif c.p==:plate && c.d==:LAF\n", + " push!(demandedAmount, c, 250)\n", + " else\n", + " push!(demandedAmount, -1)\n", + " println(\"Something is wrong index\", i, \"with customers\", c)\n", + " end\n", + "end\n", + "\n", + "costCof = [30,10, 8,10,11,71, 6,22, 7,10, 7,21,82,13,19,11,12,10,25,83,15,\n", + " 39,14,11,14,16,82, 8,27, 9,12, 9, 26,95, 17, 24,14,17,13, 28,99, 20,\n", + " 41,15,12,16,17,86, 8,29, 9,13, 9,28,99,18,26,14,17,13,31,104,20]\n", + "\n", + "# Creating costRoutes dictionary which contains\n", + "# ----------\n", + "costRoutes=Dict{Route, Float64}()\n", + "for i in [1:length(routes)]\n", + " push!(costRoutes, routes[i], costCof[i])\n", + "end\n", + "\n", + "# Creating orig, which takes the product as the input and gives the set of origins of that product\n", + "# ----\n", + "orig = Dict{Symbol, Array}()\n", + "for i in [1:length(products)]\n", + " dummy_array = Symbol[]\n", + " for j in [1:length(routes)]\n", + " #println(i, j, products[i] == routes[j].p)\n", + " if products[i] == routes[j].p\n", + " push!(dummy_array, routes[j].o)\n", + " #println(orig[products[i]])\n", + " else \n", + " #println(\"Oops, something is not right\")\n", + " end #if\n", + " end #for \n", + " push!(orig,products[i],unique(dummy_array))\n", + " end #for\n", + "\n", + "# Creating dest, which takes the product as the input and gives the set of destinations of that product\n", + "# ----\n", + "dest = Dict{Symbol, Array}()\n", + "for i in [1:length(products)]\n", + " dummy_array = Symbol[]\n", + " for j in [1:length(routes)]\n", + " #println(i, j, products[i] == routes[j].p)\n", + " if products[i] == routes[j].p\n", + " push!(dummy_array, routes[j].d)\n", + " #println(orig[products[i]])\n", + " else \n", + " #println(\"Oops, something is not right\")\n", + " end #if\n", + " end #for \n", + " push!(dest,products[i],unique(dummy_array))\n", + "end #for\n", + "\n" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 2 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Model File" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "using JuMP # Need to say it whenever we use JuMP\n", + "#using CPLEX # Telling Julia what solver we are using \n", + "using GLPKMathProgInterface\n", + "\n", + "#MODEL CONSTRUCTION\n", + "#------------------\n", + "# transpModel = Model(solver=CplexSolver()) # Name of the model object\n", + "transpModel = Model(solver=GLPKSolverLP())\n", + "\n", + "#VARIABLES\n", + "#---------\n", + "@defVar(transpModel, trans[routes] >= 0) \n", + "\n", + "#OBJECTIVE\n", + "#---------\n", + "\n", + "@setObjective(transpModel, Min, sum{costRoutes[l]*trans[l], l in routes})\n", + "\n", + "#Constraints\n", + "#-----------\n", + "\n", + "\n", + "# First Constraint\n", + "# ----------------\n", + "\n", + "for pr in products\n", + " for org in orig[pr]\n", + " @addConstraint(transpModel, sum{trans[Route(pr, org, de)], de in dest[pr]}\n", + " ==\n", + " suppliedAmount[Supply(pr,org)])\n", + " end\n", + "end\n", + "\n", + "#Second Constraint\n", + "#-----------------\n", + "\n", + "for pr in products\n", + " for de in dest[pr]\n", + " @addConstraint(transpModel, sum{trans[Route(pr, org, de)], org in orig[pr]}\n", + " ==\n", + " demandedAmount[Customer(pr,de)])\n", + " end\n", + "end\n", + "\n", + "# Final constraint:\n", + "# ----------------\n", + "for org in cities\n", + " for de in cities\n", + " if org!=de\n", + " @addConstraint(transpModel, \n", + " sum{\n", + " trans[r], r in routes\n", + " ; r.o == org && r.d==de # This will be used as an filtering condition\n", + " }\n", + " <=\n", + " capacity)\n", + " else\n", + " continue\n", + " end\n", + " end\n", + "end\n", + "\n", + "status=solve(transpModel)\n", + "\n", + "println(\"The optimal objective value is: \", getObjectiveValue(transpModel))\n", + "\n", + "println(\"The optimal solution is, trans= \\n\", getValue(trans))" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "The optimal objective value is: " + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "199500.0\n", + "The optimal solution is, trans= \n" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "trans: 1 dimensions, 63 entries:\n", + " [Route(:bands,:CLEV,:LAN)] = 0.0\n", + " [Route(:bands,:GARY,:FRA)] = 0.0\n", + " [Route(:bands,:CLEV,:FRE)] = 0.0\n", + " [Route(:plate,:CLEV,:DET)] = 100.0\n", + " [Route(:coils,:PITT,:WIN)] = 0.0\n", + " [Route(:plate,:CLEV,:LAN)] = -0.0\n", + " [Route(:bands,:PITT,:FRE)] = 225.0\n", + " [Route(:plate,:GARY,:WIN)] = 0.0\n", + " [Route(:coils,:GARY,:LAF)] = 150.0\n", + " [Route(:coils,:CLEV,:STL)] = 300.0\n", + " [Route(:plate,:CLEV,:STL)] = -0.0\n", + " [Route(:bands,:PITT,:LAN)] = 100.0\n", + " [Route(:coils,:GARY,:FRE)] = 625.0\n", + " [Route(:plate,:CLEV,:FRA)] = 0.0\n", + " [Route(:coils,:CLEV,:FRA)] = 0.0\n", + " [Route(:bands,:PITT,:LAF)] = 250.0\n", + " [Route(:coils,:CLEV,:FRE)] = 50.0\n", + " [Route(:bands,:GARY,:DET)] = 0.0\n", + " [Route(:plate,:GARY,:FRA)] = 0.0\n", + " [Route(:coils,:CLEV,:LAF)] = 175.0\n", + " [Route(:bands,:PITT,:WIN)] = 0.0\n", + " [Route(:plate,:PITT,:FRE)] = 0.0\n", + " [Route(:coils,:CLEV,:DET)] = 425.0\n", + " [Route(:plate,:GARY,:FRE)] = 0.0\n", + " [Route(:bands,:GARY,:FRE)] = 0.0\n", + " [Route(:bands,:PITT,:DET)] = 200.0\n", + " [Route(:bands,:GARY,:LAN)] = 0.0\n", + " [Route(:bands,:PITT,:FRA)] = 25.0\n", + " [Route(:coils,:PITT,:LAF)] = 175.0\n", + " [Route(:plate,:CLEV,:LAF)] = 50.0\n", + " [Route(:bands,:CLEV,:STL)] = 250.0\n", + " [Route(:coils,:GARY,:WIN)] = 0.0\n", + " [Route(:plate,:GARY,:STL)] = 200.0\n", + " [Route(:bands,:CLEV,:WIN)] = 75.0\n", + " [Route(:coils,:PITT,:FRA)] = 500.0\n", + " [Route(:bands,:CLEV,:FRA)] = 275.0\n", + " [Route(:plate,:PITT,:FRA)] = 100.0\n", + " [Route(:plate,:GARY,:LAF)] = 0.0\n", + " [Route(:coils,:CLEV,:LAN)] = 400.0\n", + " [Route(:plate,:PITT,:STL)] = 0.0\n", + " [Route(:coils,:CLEV,:WIN)] = 250.0\n", + " [Route(:coils,:PITT,:DET)] = 325.0\n", + " [Route(:coils,:GARY,:LAN)] = 0.0\n", + " [Route(:plate,:PITT,:LAF)] = 200.0\n", + " [Route(:plate,:PITT,:DET)] = 0.0\n", + " [Route(:plate,:CLEV,:WIN)] = 50.0\n", + " [Route(:coils,:GARY,:DET)] = 0.0\n", + " [Route(:bands,:PITT,:STL)] = 0.0\n", + " [Route(:plate,:CLEV,:FRE)] = 100.0\n", + " [Route(:bands,:GARY,:WIN)] = 0.0\n", + " [Route(:plate,:PITT,:LAN)] = 0.0\n", + " [Route(:coils,:GARY,:FRA)] = 0.0\n", + " [Route(:plate,:GARY,:LAN)] = 0.0\n", + " [Route(:plate,:PITT,:WIN)] = 0.0\n", + " [Route(:bands,:CLEV,:DET)] = 100.0\n", + " [Route(:coils,:PITT,:STL)] = 625.0\n", + " [Route(:bands,:GARY,:LAF)] = 0.0\n", + " [Route(:bands,:CLEV,:LAF)] = 0.0\n", + " [Route(:plate,:GARY,:DET)] = 0.0\n", + " [Route(:coils,:PITT,:LAN)] = 0.0\n", + " [Route(:coils,:PITT,:FRE)] = 175.0\n", + " [Route(:coils,:GARY,:STL)] = 25.0\n", + " [Route(:bands,:GARY,:STL)] = 400.0" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n" + ] + } + ], + "prompt_number": 8 + } + ], + "metadata": {} + } + ] +} \ No newline at end of file diff --git a/notebooks/Shuvomoy - Getting started with JuMP.ipynb b/notebooks/Shuvomoy - Getting started with JuMP.ipynb new file mode 100644 index 0000000..0cb85f4 --- /dev/null +++ b/notebooks/Shuvomoy - Getting started with JuMP.ipynb @@ -0,0 +1,1601 @@ +{ + "metadata": { + "language": "Julia", + "name": "", + "signature": "sha256:95d8d8e322ebd34d3bef9e9d130b7f2b6909feb2b0c0669a6b936075c2beae59" + }, + "nbformat": 3, + "nbformat_minor": 0, + "worksheets": [ + { + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Description:** This notebook is an introduction to JuMP. The topics described are as follows:\n", + "- Installing Julia and JuMP\n", + "- Representing vectors in Julia\n", + "- Structure of a JuMP model\n", + "- Solving general purpose linear programming problem\n", + "- Solving general purpose integer programming problem\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 - getting started\n", + "--------------------------\n", + "\n", + "## Julia\n", + "\n", + "\n", + "### What is Julia?\n", + "\n", + "Julia is a new programming language.\n", + "\n", + "### Why Julia?\n", + "\n", + "- Free and open-source\n", + "- Syntax similar to MATLAB, speed similar to C\n", + "\n", + "\n", + "## JuMP\n", + "\n", + "### What is JuMP?\n", + "\n", + "JuMP is a modelling language for mathematical optimization [1]. It is embedded in Julia.\n", + "\n", + "### Why JuMP?\n", + "\n", + "- Very user-friendly\n", + "- Speed similar to special purpose commercial modelling language like AMPL\n", + "- Solver independent code: same code will run for both commercial and open-source solvers\n", + "- Very easy to implement solver callback and problem modification\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Installation\n", + "\n", + "## Installing Julia and IJulia\n", + "\n", + "### Installing Julia\n", + "\n", + "Go to the [download page](http://julialang.org/downloads/ \"downloadJulia\") of Julia, download the appropriate installer and run it. Now you can start an interactive Julia session! \n", + "\n", + "\n", + "\n", + "### Installing IJulia\n", + "\n", + "IJulia will allow us to create powerful graphical notebooks, which is very convenient for the tutorial. We can insert codes, text, mathematical formulas and multimedia etc in the same notebook. To install IJulia: \n", + "\n", + "- Install Anaconda from http://continuum.io/downloads. If you are on Windows, then while running the Anaconda installer please check the options *Add Anaconda to the System Path* and also *Register Anaconda as default Python version of the system*.\n", + "\n", + "- Now start a Julia interactive session. Type the following `Pkg.add(\"IJulia\")`\n", + "\n", + "- To open a new notebook\n", + " * In the Julia interactive session, run `Using IJulia` and then `notebook()`. It will open the Home for IJulia in your web browser. \n", + " \n", + " \n", + " \n", + " * The directory location can be checked by the command `pwd()`. If you want to change the directory to something else, the before running `Using IJulia` and `notebook()`, run `cd(path to your preferred directory)`, e.g., `cd(\"E:\\\\Dropbox\\\\Julia_Workspaces\")`. \n", + " * Click on *New Notebook* on the top right corner. In the new notebook, you can execute any Julia command pressing SHIFT+ENTER.\n", + " \n", + " \n", + " \n", + "### Installing JuMP\n", + " \n", + "At first add the JuMP package by running the following code in the notebook:\n" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "Pkg.add(\"JuMP\")" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stderr", + "text": [ + "INFO: Nothing to be done\n" + ] + } + ], + "prompt_number": 2 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We need to install a solver package. Let's install the opensource solvers GLPK, Cbc and Clp by typing in `Pkg.add(\"GLPKMathProgInterface\")`, `Pkg.add(\"Cbc\")` and `Pkg.add(\"Clp\")` respectively. Let's add the Julia package associated with CPLEX by typing in `Pkg.add(\"CPLEX\")`. The other choices are `\"CPLEX\"`, `\"Cbc\"`, `\"Clp\"`, `\"Gurobi\"` and `\"MOSEK\"`. \n", + "\n", + "It should be noted that, in order to use commercial solvers such as CPLEX, Gurobi and Mosek in JuMP, we will require working installations of them with appropriate licences. Both Gurobi and Mosek are free for academic use. CPLEX is free for faculty members and graduate teaching assistants. " + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "Pkg.add(\"GLPKMathProgInterface\")" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stderr", + "text": [ + "INFO: Nothing to be done\n" + ] + } + ], + "prompt_number": 4 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "Pkg.add(\"Cbc\")" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stderr", + "text": [ + "INFO: Nothing to be done\n" + ] + } + ], + "prompt_number": 6 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "Pkg.add(\"Clp\")" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stderr", + "text": [ + "INFO: Nothing to be done\n" + ] + } + ], + "prompt_number": 8 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "Pkg.add(\"CPLEX\")" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stderr", + "text": [ + "INFO: Nothing to be done\n" + ] + } + ], + "prompt_number": 11 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "Pkg.add(\"Gurobi\")" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stderr", + "text": [ + "INFO: Nothing to be done\n" + ] + } + ], + "prompt_number": 14 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If you have not updated your Julia packages in a while, a good idea might be updating them." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "Pkg.update()" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stderr", + "text": [ + "INFO: Updating METADATA...\n" + ] + }, + { + "output_type": "stream", + "stream": "stderr", + "text": [ + "INFO: Computing changes...\n" + ] + }, + { + "output_type": "stream", + "stream": "stderr", + "text": [ + "INFO: No packages to install, update or remove\n" + ] + } + ], + "prompt_number": 12 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "println(\"Hello World!\")" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Hello World!" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n" + ] + } + ], + "prompt_number": 13 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# The very first example\n", + "\n", + "At first let us try to solve a very simple and trivial optimization problem using JuMP to check if everything is working properly. \n", + "\n", + "\n", + "\\begin{align}\n", + "\\text{minimize} \\qquad & x+y \\\\\n", + " \\text{subject to} \\quad \\quad & x+y \\leq 1 \\\\\n", + " \\qquad \\qquad & x \\geq 0, y \\geq 0 \\\\\n", + " \\qquad \\qquad & x,y \\in \\mathbb{R}\n", + "\\end{align}\n", + "\n", + "\n", + "Here is the JuMP code to solve the mentioned problem:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "using JuMP # Need to say it whenever we use JuMP\n", + "\n", + "using GLPKMathProgInterface # Loading the GLPK module for using its solver\n", + "\n", + "\n", + "#MODEL CONSTRUCTION\n", + "#--------------------\n", + "\n", + "myModel = Model(solver=GLPKSolverLP()) \n", + "# Name of the model object. All constraints and variables of an optimization problem are associated \n", + "# with a particular model object. The name of the model object does not have to be myModel, it can be yourModel too! The argument of Model,\n", + "# solver=GLPKsolverLP() means that to solve the optimization problem we will use GLPK solver.\n", + "\n", + "#VARIABLES\n", + "#---------\n", + "\n", + "# A variable is modeled using @defVar(name of the model object, variable name and bound, variable type)\n", + "# Bound can be lower bound, upper bound or both. If no variable type is defined, then it is treated as \n", + "#real. For binary variable write Bin and for integer use Int.\n", + "\n", + "@defVar(myModel, x >= 0) # Models x >=0\n", + "\n", + "# Some possible variations:\n", + "# @defVar(myModel, x, Binary) # No bound on x present, but x is a binary variable now\n", + "# @defVar(myModel, x <= 10) # This one defines a variable with lower bound x <= 10\n", + "# @defVar(myModel, 0 <= x <= 10, Int) # This one has both lower and upper bound, and x is an integer\n", + "\n", + "@defVar(myModel, y >= 0) # Models y >= 0\n", + "\n", + "#OBJECTIVE\n", + "#---------\n", + "\n", + "@setObjective(myModel, Min, x + y) # Sets the objective to be minimized. For maximization use Max\n", + "\n", + "#CONSTRAINTS\n", + "#-----------\n", + "\n", + "@addConstraint(myModel, x + y <= 1) # Adds the constraint x + y <= 1\n", + "\n", + "#THE MODEL IN A HUMAN-READABLE FORMAT\n", + "#------------------------------------\n", + "println(\"The optimization problem to be solved is:\")\n", + "print(myModel) # Shows the model constructed in a human-readable form\n", + "\n", + "#SOLVE IT AND DISPLAY THE RESULTS\n", + "#--------------------------------\n", + "status = solve(myModel) # solves the model \n", + "\n", + "println(\"Objective value: \", getObjectiveValue(myModel)) # getObjectiveValue(model_name) gives the optimum objective value\n", + "println(\"x = \", getValue(x)) # getValue(decision_variable) will give the optimum value of the associated decision variable\n", + "println(\"y = \", getValue(y))" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "The optimization problem to be solved is:" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n", + "Min x + y\n", + "Subject to\n", + " x + y <= 1\n", + " x >= 0\n", + " y >= 0\n" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Objective value: " + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "0.0\n", + "x = 0.0\n", + "y = 0.0\n" + ] + } + ], + "prompt_number": 3 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This was certainly not the most exciting optimization problem to solve. This was for test purpose only. However, before going into the structure of a JuMP model, let us learn how to represent vectors in Julia.\n", + " \n", + "# Representing vectors in Julia\n", + " \n", + " * A column vector, $y=(y_1, y_2, \\ldots, y_n)= \\begin{pmatrix}\n", + " y_1 \\\\\n", + " y_2 \\\\\n", + " . \\\\\n", + " . \\\\\n", + " y_n\n", + " \\end{pmatrix} \\in \\mathbb{R}^n$ will be written in Julia as `[y[1];y[2];...;y[n]]`. \n", + " \n", + " For example to create column vector $\\begin{pmatrix}\n", + " 3 \\\\\n", + " 2.4 \\\\\n", + " 9.1 \\\\\n", + " \\end{pmatrix}$ use: `[3; 2.4; 9.1]`.\n" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "[3; 2.4; 9.1] # Column vector" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 19, + "text": [ + "3-element Array{Float64,1}:\n", + " 3.0\n", + " 2.4\n", + " 9.1" + ] + } + ], + "prompt_number": 19 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "* A row vector, $z=(z_1 \\; z_2 \\; \\ldots \\; z_n) \\in \\mathbb{R}^{1 \\times n}$ will be written in Julia as `[z[1] y[2]...z[n]]`. \n", + " \n", + " For example to create row vector $(1.2 \\; 3.5 \\; 8.21)$ use: `[1.2 3.5 8.21]`.\n", + " " + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "[1.2 3.5 8.21] # Row vector" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 11, + "text": [ + "1x3 Array{Float64,2}:\n", + " 1.2 3.5 8.21" + ] + } + ], + "prompt_number": 11 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "* To create a $m \\times n$ matrix $A = \\begin{pmatrix}\n", + " A_{11} & A_{12} & A_{13} & \\ldots &A_{1n} \\\\\n", + " \\ldots & \\ldots & \\ldots & \\ldots & \\ldots \\\\\n", + " A_{m1} & A_{m2} & A_{m3} & \\ldots & A_{mn}\n", + " \\end{pmatrix}$ write:\n", + " \n", + " `[A[1,1] A[1,2] A[1,3]... A[1,n];`\n", + " ` ... ; `\n", + " `A[m,1] A[m,2] ... A[m,n]]`.\n", + " \n", + " So the matrix $A = \\begin{pmatrix}\n", + " 1 & 1 & 9 & 5 \\\\\n", + " 3 & 5 & 0 & 8 \\\\\n", + " 2 & 0 & 6 & 13\n", + " \\end{pmatrix}$ is represented in Julia as: \n", + " \n", + " `A= [\n", + " 1 1 9 5;\n", + " 3 5 0 8;\n", + " 2 0 6 13\n", + " ]`" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# Generating a matrix\n", + "A= [\n", + " 1 1 9 5;\n", + " 3 5 0 8;\n", + " 2 0 6 13\n", + " ]" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 12, + "text": [ + "3x4 Array{Int64,2}:\n", + " 1 1 9 5\n", + " 3 5 0 8\n", + " 2 0 6 13" + ] + } + ], + "prompt_number": 12 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " $A_{ij}$ can be accessed by `A[i,j]` ,the $i$th row of the matrix A is represented by `A[i,:]`, the $j$th column of the matrix A is represented by `A[:,j]`. \n", + " \n", + " The size of a matrix $A$ can be determined by running the command `size(A)`. If we write `numRows, numCols = size(A)`, then `numRows` and `numCols` will contain the total number of rows and columns of A respectively." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "numRows, numCols = size(A)\n", + "println(\n", + "\"A has \", numRows, \" rows and \", numCols, \" columns \\n\",\n", + "\"A[3,3] is \", A[3,3], \"\\n\",\n", + "\"The 3rd column of A is \", A[:,3], \"\\n\",\n", + "\"The 2nd row of A is \", A[2,:]\n", + ")\n" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "A has " + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "3 rows and 4 columns \n", + "A[3,3] is 6\n", + "The 3rd column of A is [9,0,6]\n", + "The 2nd row of A is " + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "[3 5 0 8]\n" + ] + } + ], + "prompt_number": 13 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Suppose $x,y \\in \\mathbb{R}^n$. Then $x^T y =\\sum_{i=1}^{n} {x_i y_i}$ is written as `dot(x,y)`. " + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "y=[1; 2; 3; 4]\n", + "x=[5; 6; 7; 8]\n", + "xTy=dot(x,y)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 14, + "text": [ + "70" + ] + } + ], + "prompt_number": 14 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Structure of a JuMP model\n", + "\n", + "Any JuMP model that describes an optimization problem must have four parts: \n", + "\n", + "- **Model Object**, \n", + "- **Variables**, \n", + "- **Objective**, \n", + "- **Constraints**.\n", + "\n", + "## Model\n", + "\n", + "Any instance of an optimization problem corresponds to a model object. This model object is associated with all the variables, constraints and objective of the instance. It is constructed using `modelName = Model(solver=`*solver of our preference*`)`. If no solver is specified, then `ClpSolver()` and/or `CbcSolver()` will be used by default. Here `modelName` is any valid name. We will limit ourselves in the open source solvers such as:\n", + "\n", + "* Linear Programming Solver: `ClpSolver(), GLPKSolverLP()`\n", + "* Mixed Integer Programming Solver: `GLPKSolverMIP() CbcSolver()`\n" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "using JuMP\n", + "myModel = Model() # ClpSolver() and/or CbcSolver() will be used based on the problem" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "latex": [ + "$$ \\begin{alignat*}{1}\\min\\quad & 0\\\\\n", + "\\text{Subject to} \\quad\\end{alignat*}\n", + " $$" + ], + "metadata": {}, + "output_type": "pyout", + "prompt_number": 2, + "text": [ + "Feasibility problem with:\n", + " * 0 linear constraints\n", + " * 0 variables\n", + "Solver set to Default" + ] + } + ], + "prompt_number": 2 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n", + "##Variables\n", + "Variables are defined using `@defVar` macro, which takes upto three input arguments. The *first* arguement is the name of the model. Then the *second* argument contains the name of the variable, and a bound on the variable if it exists. The *third* argument is not needed if the variable is real. When the variable is binary or integer, then `Bin` or `Int`, respectively, is used in place of the third argument.\n", + "\n", + "\n", + "\n", + "###Examples of Variables\n", + "Suppose the model object is `myModel`. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "- To describe a variable $z \\in \\mathbb{R}$ such that $0 \\leq z \\leq 10$\n", + "write" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "@defVar(myModel, 0 <= z <= 10)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "latex": [ + "$$ z $$" + ], + "metadata": {}, + "output_type": "pyout", + "prompt_number": 4, + "text": [ + "z" + ] + } + ], + "prompt_number": 4 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "- Now consider a decision variable $x \\in \\mathbb{R}^n$, and it has a bound $l \\preceq x \\preceq u$, where naturally $l, u \\in \\mathbb{R}^n$. For that we write
\n" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# INPUT DATA, CHANGE THEM TO YOUR REQUIREMENT\n", + "#-------------------------------------------\n", + "n = 10\n", + "l = [1; 2; 3; 4; 5; 6; 7; 8; 9; 10]\n", + "u = [10; 11; 12; 13; 14; 15; 16; 17; 18; 19]" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 16, + "text": [ + "10-element Array{Int64,1}:\n", + " 10\n", + " 11\n", + " 12\n", + " 13\n", + " 14\n", + " 15\n", + " 16\n", + " 17\n", + " 18\n", + " 19" + ] + } + ], + "prompt_number": 16 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# VARIABLE DEFINITION\n", + "# ------------------- \n", + "@defVar(myModel, l[i] <= x[i=1:n] <= u [i])" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "latex": [ + "$$ .. \\leq x_{i} \\leq .. \\quad\\forall i \\in \\{1,2,\\dots,9,10\\} $$" + ], + "metadata": {}, + "output_type": "pyout", + "prompt_number": 17, + "text": [ + ".. <= x[i] <= .. for all i in {1,2..9,10}" + ] + } + ], + "prompt_number": 17 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "- Suppose we have decision variables $x \\in \\mathbb{R}^n$, $y \\in \\mathbb{Z}^m$ and $z \\in \\mathbb \\{0,1\\}^p$ such that $x \\succeq 0$, $a \\preceq y \\preceq b$. Here $a, b \\in \\mathbb{Z}^m$. To express this in JuMP we write" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# INPUT DATA, CHANGE THEM TO YOUR REQUIREMENT\n", + "#-------------------------------------------\n", + "n = 4 # dimension of x\n", + "m = 3 # dimension of y\n", + "p = 2 # dimensin of z\n", + "a = [0; 1; 2]\n", + "b = [3; 4; 7]" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 18, + "text": [ + "3-element Array{Int64,1}:\n", + " 3\n", + " 4\n", + " 7" + ] + } + ], + "prompt_number": 18 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "# VARIABLE DEFINITION\n", + "# -------------------\n", + "@defVar(myModel, x[i=1:n] >= 0)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "latex": [ + "$$ x_{i} \\geq 0 \\quad\\forall i \\in \\{1,2,3,4\\} $$" + ], + "metadata": {}, + "output_type": "pyout", + "prompt_number": 19, + "text": [ + "x[i] >= 0 for all i in {1,2,3,4}" + ] + } + ], + "prompt_number": 19 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "@defVar(myModel, a[i] <= y[i=1:m] <= b[i], Int)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "latex": [ + "$$ .. \\leq y_{i} \\leq .., \\in \\mathbb{Z}, \\quad\\forall i \\in \\{1,2,3\\} $$" + ], + "metadata": {}, + "output_type": "pyout", + "prompt_number": 20, + "text": [ + ".. <= y[i] <= .., integer, for all i in {1,2,3}" + ] + } + ], + "prompt_number": 20 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "@defVar(myModel, z[i=1:p], Bin)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "latex": [ + "$$ z_{i} \\in \\{0,1\\} \\quad\\forall i \\in \\{1,2\\} $$" + ], + "metadata": {}, + "output_type": "pyout", + "prompt_number": 21, + "text": [ + "z[i] in {0,1} for all i in {1,2}" + ] + } + ], + "prompt_number": 21 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Constraints\n", + "Constraints are added by using `@addConstraint` macro. The first argument is the model object the constraint is associated with, the second argument is the reference to that constraint and the third argument is the constraint description. The constraint reference comes handy when we want to manipulate the constraint later or access the dual variables associated with it. If no constraint reference is needed, then the second argument is the constraint description.
" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "###Examples of Constraints\n", + "Let's give some examples on writing constraints in JuMP. Suppose the model name is `yourModel`." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "yourModel = Model()" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "latex": [ + "$$ \\begin{alignat*}{1}\\min\\quad & 0\\\\\n", + "\\text{Subject to} \\quad\\end{alignat*}\n", + " $$" + ], + "metadata": {}, + "output_type": "pyout", + "prompt_number": 5, + "text": [ + "Feasibility problem with:\n", + " * 0 linear constraints\n", + " * 0 variables\n", + "Solver set to Default" + ] + } + ], + "prompt_number": 5 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "1. Consider variables $x, y \\in \\mathbb{R}$ which are coupled by the constraints $5 x +3 y \\leq 5$. We write this as
\n", + "`@addConstraint(yourModel, 5*x + 3*y <= 5)`
\n", + "Naturally, `x` and `y` have to be defined first using `@defVar` macro." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "@defVar(yourModel, x)\n", + "@defVar(yourModel, y)\n", + "@addConstraint(yourModel, 5*x + 3*y <= 5)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "latex": [ + "$$ 5 x + 3 y \\leq 5 $$" + ], + "metadata": {}, + "output_type": "pyout", + "prompt_number": 6, + "text": [ + "5 x + 3 y <= 5" + ] + } + ], + "prompt_number": 6 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + " Here no constraint reference is given. Now suppose we want to get the dual value of some constraint after solving the problem, then we would need a constraint reference to assign to the constraint first. Let's call the constraint reference as `conRef1` (it could be any valid name). Then the same constraint have to be written as:
\n", + "`@addConstraint(yourModel, conRef1, 6*x + 4*y >= 5)`
\n", + "When we would need the dual value after solving the problem we just write `println(getDual(conRef1))`.
" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "@addConstraint(yourModel, conRef1, 6*x + 4*y >= 5) " + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "latex": [ + "$$ 6 x + 4 y \\geq 5 $$" + ], + "metadata": {}, + "output_type": "pyout", + "prompt_number": 7, + "text": [ + "6 x + 4 y >= 5" + ] + } + ], + "prompt_number": 7 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "2. Consider a variable $x \\in \\mathbb{R}^4$, a coefficient vector $a=(1, -3, 5, -7)$ We want to write a constraint of the form $\\sum_{i=1}^4{a_i x_i} \\leq 3$. In JuMP we write:
\n" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "a = [1; -3; 5; 7] \n", + "@defVar(yourModel, w[1:4])\n", + "@addConstraint(yourModel, sum{a[i]*w[i],i=1:4} <= 3)" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "latex": [ + "$$ w_{1} - 3 w_{2} + 5 w_{3} + 7 w_{4} \\leq 3 $$" + ], + "metadata": {}, + "output_type": "pyout", + "prompt_number": 8, + "text": [ + "w[1] - 3 w[2] + 5 w[3] + 7 w[4] <= 3" + ] + } + ], + "prompt_number": 8 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##Objective\n", + "Objective is set using `@setObjective` macro. It has three arguments. The first argument is as usual the model object. The second one is either `Max` if we want to maximize the objective function, or `Min` when we want to minimize. The last argument is the description of the objective which has similar syntax to that of constraint definition.\n", + "\n", + "###Example of objective\n", + "For the previous model, consider the decision variable $w \\in \\mathbb{R}^4$ and cost vector $c = (2, 3 , 4, 5)$. We want to minimize $c^T w$. In JuMP we would write:" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "c = [2; 3; 4; 5] \n", + "@setObjective(yourModel, Min, sum{c[i]*w[i],i=1:4})" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "latex": [ + "$$ 2 w_{1} + 3 w_{2} + 4 w_{3} + 5 w_{4} $$" + ], + "metadata": {}, + "output_type": "pyout", + "prompt_number": 10, + "text": [ + "2 w[1] + 3 w[2] + 4 w[3] + 5 w[4]" + ] + } + ], + "prompt_number": 10 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##Solving a standard form Linear Programming \n", + "problem\n", + "Let us try to write the JuMP code for the following standard form optimization problem:\n", + "\n", + "\n", + "\\begin{align}\n", + "& \\text{minimize} && c^T x \\\\\n", + "& \\text{subject to} && A x = b \\\\\n", + " & && x \\succeq 0 \\\\\n", + " & && x \\in \\mathbb{R}^n\n", + "\\end{align}\n", + "\n", + "\n", + "where, $n = 4$, $c=(1, 3, 5, 2)$, $A = \\begin{pmatrix}\n", + " 1 & 1 & 9 & 5 \\\\\n", + " 3 & 5 & 0 & 8 \\\\\n", + " 2 & 0 & 6 & 13\n", + " \\end{pmatrix}$ and $b=(7, 3, 5)$. The symbol $\\succeq$ ($\\preceq$) stands for element-wise greater (less) than or equal to.\n", + "\n", + "### Entering different parts of the code one by one\n", + "Let us input differnt parts of the JuMP code one by one and see the corresponding outputs to detect if everything is okay. Of course we could input the whole code at once." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "using JuMP # Need to say it whenever we use JuMP\n", + "\n", + "using GLPKMathProgInterface # Loading the package for using the GLPK solver\n" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 72 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "#MODEL CONSTRUCTION\n", + "#------------------\n", + "\n", + "sfLpModel = Model(solver=GLPKSolverLP()) # Name of the model object" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "latex": [ + "$$ \\begin{alignat*}{1}\\min\\quad & 0\\\\\n", + "\\text{Subject to} \\quad\\end{alignat*}\n", + " $$" + ], + "metadata": {}, + "output_type": "pyout", + "prompt_number": 74, + "text": [ + "Feasibility problem with:\n", + " * 0 linear constraints\n", + " * 0 variables\n", + "Solver set to GLPK" + ] + } + ], + "prompt_number": 74 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "#INPUT DATA\n", + "#----------\n", + "\n", + "c = [1; 3; 5; 2] \n", + "\n", + "A= [\n", + " 1 1 9 5;\n", + " 3 5 0 8;\n", + " 2 0 6 13\n", + " ]\n", + "\n", + "b = [7; 3; 5] \n", + "\n", + "m, n = size (A) # m = number of rows of A, n = number of columns of A" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 75, + "text": [ + "(3,4)" + ] + } + ], + "prompt_number": 75 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "#VARIABLES\n", + "#---------\n", + "\n", + "@defVar(sfLpModel, x[1:n] >= 0) # Models x >=0\n" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "latex": [ + "$$ x_{i} \\geq 0 \\quad\\forall i \\in \\{1,2,3,4\\} $$" + ], + "metadata": {}, + "output_type": "pyout", + "prompt_number": 76, + "text": [ + "x[i] >= 0 for all i in {1,2,3,4}" + ] + } + ], + "prompt_number": 76 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "#CONSTRAINTS\n", + "#-----------\n", + "\n", + "for i=1:m # for all rows do the following\n", + " @addConstraint(sfLpModel, sum{A[i,j]*x[j] , j=1:n} == b[i]) # the ith row \n", + " # of A*x is equal to the ith component of b\n", + "end # end of the for loop" + ], + "language": "python", + "metadata": {}, + "outputs": [], + "prompt_number": 77 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "#OBJECTIVE\n", + "#---------\n", + "\n", + "@setObjective(sfLpModel, Min, sum{c[j]*x[j], j=1:n}) # minimize c'x " + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "latex": [ + "$$ x_{1} + 3 x_{2} + 5 x_{3} + 2 x_{4} $$" + ], + "metadata": {}, + "output_type": "pyout", + "prompt_number": 78, + "text": [ + "x[1] + 3 x[2] + 5 x[3] + 2 x[4]" + ] + } + ], + "prompt_number": 78 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "#THE MODEL IN A HUMAN-READABLE FORMAT\n", + "#------------------------------------\n", + "\n", + "println(\"The optimization problem to be solved is:\")\n", + "print(sfLpModel) # Shows the model constructed in a human-readable form" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "The optimization problem to be solved is:" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n", + "Min x[1] + 3 x[2] + 5 x[3] + 2 x[4]\n", + "Subject to\n", + " x[1] + x[2] + 9 x[3] + 5 x[4] == 7\n", + " 3 x[1] + 5 x[2] + 8 x[4] == 3\n", + " 2 x[1] + 6 x[3] + 13 x[4] == 5\n", + " x[i] >= 0 for all i in {1,2,3,4}\n" + ] + } + ], + "prompt_number": 79 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "status = solve(sfLpModel) # solves the model " + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 80, + "text": [ + ":Optimal" + ] + } + ], + "prompt_number": 80 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "#SOLVE IT AND DISPLAY THE RESULTS\n", + "#--------------------------------\n", + "\n", + "println(\"Objective value: \", getObjectiveValue(sfLpModel)) # getObjectiveValue(model_name) gives the optimum objective value\n", + "\n", + "println(\"Optimal solution is x = \\n\", getValue(x)) # getValue(decision_variable) will give the optimum value \n", + " # of the associated decision variable\n" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Objective value: " + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "4.923076923076923\n" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Optimal solution is x = \n" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "[1] = 0.4230769230769228\n", + "[2] = 0.34615384615384626\n", + "[3] = 0.6923076923076923\n", + "[4] = 0.0\n", + "\n" + ] + } + ], + "prompt_number": 81 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### The whole code" + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "using JuMP \n", + "\n", + "using GLPKMathProgInterface \n", + "\n", + "#MODEL CONSTRUCTION\n", + "#------------------\n", + "\n", + "sfLpModel = Model(solver=GLPKSolverLP())\n", + "\n", + "#INPUT DATA\n", + "#----------\n", + "\n", + "c = [1; 3; 5; 2] \n", + "\n", + "A= [\n", + " 1 1 9 5;\n", + " 3 5 0 8;\n", + " 2 0 6 13\n", + " ]\n", + "\n", + "b = [7; 3; 5] \n", + "\n", + "m, n = size (A) # m = number of rows of A, n = number of columns of A\n", + "\n", + "#VARIABLES\n", + "#---------\n", + "\n", + "@defVar(sfLpModel, x[1:n] >= 0) # Models x >=0\n", + "\n", + "#CONSTRAINTS\n", + "#-----------\n", + "\n", + "for i=1:m # for all rows do the following\n", + " @addConstraint(sfLpModel, sum{A[i,j]*x[j] , j=1:n} == b[i]) # the ith row \n", + " # of A*x is equal to the ith component of b\n", + "end # end of the for loop\n", + "\n", + "#OBJECTIVE\n", + "#---------\n", + "\n", + "@setObjective(sfLpModel, Min, sum{c[j]*x[j], j=1:n}) # minimize c'x \n", + "\n", + "\n", + "\n", + "#THE MODEL IN A HUMAN-READABLE FORMAT\n", + "#------------------------------------\n", + "\n", + "println(\"The optimization problem to be solved is:\")\n", + "print(sfLpModel) # Shows the model constructed in a human-readable form\n", + "\n", + "@time begin\n", + "status = solve(sfLpModel) # solves the model \n", + "end\n", + "#SOLVE IT AND DISPLAY THE RESULTS\n", + "#-------------------------------- \n", + "\n", + "println(\"Objective value: \", getObjectiveValue(sfLpModel)) # getObjectiveValue(model_name) gives the optimum objective value\n", + "\n", + "println(\"Optimal solution is x = \\n\", getValue(x)) # getValue(decision_variable) will give the optimum value \n", + " # of the associated decision variable\n" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "The optimization problem to be solved is:" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Min x[1] + 3 x[2] + 5 x[3] + 2 x[4]\n", + "Subject to\n", + " x[1] + x[2] + 9 x[3] + 5 x[4] == 7\n", + " 3 x[1] + 5 x[2] + 8 x[4] == 3\n", + " 2 x[1] + 6 x[3] + 13 x[4] == 5\n", + " x[i] >= 0 for all i in {1,2,3,4}\n", + "elapsed time: 0.000363692 seconds (6792 bytes allocated)\n", + "Objective value: 4.923076923076923\n", + "Optimal solution is x = \n", + "x: 1 dimensions:\n", + "[1] = 0.4230769230769228\n", + "[2] = 0.34615384615384626\n", + "[3] = 0.6923076923076923\n", + "[4] = 0.0\n", + "\n" + ] + } + ], + "prompt_number": 5 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Solving a standard form Mixed Integer Programming Problem\n", + "\n", + "Let us try to write the JuMP code for the following standard form optimization problem:\n", + "\n", + "$\n", + "\\begin{align}\n", + "& \\text{minimize} && c^T x + d^T y\\\\\n", + "& \\text{subject to} && A x + B y= f \\\\\n", + " & && x \\succeq 0, y \\succeq 0 \\\\\n", + " & && x \\in \\mathbb{R}^n, y \\in \\mathbb{Z}^p\n", + "\\end{align}\n", + "$\n", + "\n", + "Here, $A \\in \\mathbb{R}^{m \\times n}, B \\in \\mathbb{R}^{m \\times p}, c \\in \\mathbb{R}^n, d \\in \\mathbb{R}^p, f \\in \\mathbb{R}^m$. The data were randomly generated. The symbol $\\succeq$ ($\\preceq$) stands for element-wise greater (less) than or equal to." + ] + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "n = 5\n", + "p = 4\n", + "m = 3\n", + "A=\n", + "[0.7511 -0.1357 0.7955 -0.4567 0.1356\n", + "-0.6670 -0.3326 0.1657 -0.5519 -0.9367\n", + " 1.5894 -0.1302 -0.4313 -0.4875 0.4179]\n", + "\n", + "B=\n", + "[-0.09520 -0.28056 -1.33978 0.6506\n", + " -0.8581 -0.3518 1.2788 1.5114\n", + " -0.5925 1.3477 0.1589 0.03495]\n", + "\n", + "c=[0.3468,0.8687,0.1200,0.5024,0.2884]\n", + "\n", + "d=[0.2017,0.2712,0.4997,0.9238]\n", + "\n", + "f = [0.1716,0.3610,0.0705]\n" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "metadata": {}, + "output_type": "pyout", + "prompt_number": 27, + "text": [ + "3-element Array{Float64,1}:\n", + " 0.1716\n", + " 0.361 \n", + " 0.0705" + ] + } + ], + "prompt_number": 27 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "using JuMP\n", + "using GLPKMathProgInterface\n", + "\n", + "sfMipModel = Model(solver = GLPKSolverMIP())\n", + "\n", + "@defVar(sfMipModel, x[1:n] >=0)\n", + "@defVar(sfMipModel, y[1:p] >= 0, Int)\n", + "\n", + "@setObjective(sfMipModel, Min, sum{c[i] * x[i], i in 1:n}+sum{d[i]*y[i], i in 1:p})\n", + "\n", + "for i in [1:m]\n", + " @addConstraint(sfMipModel, sum{A[i,j]*x[j], j in [1:n]}+ sum{B[i,j]*y[j], j in [1:p]} == f[i])\n", + "end\n", + "\n", + "print(sfMipModel, \"\\n\")\n", + "statusMipModel = solve(sfMipModel)\n", + "print(\"Status of the problem is \", statusMipModel, \"\\n\")\n", + "\n", + "if statusMipModel == :Optimal\n", + " print (\"Optimal objective value = \", getObjectiveValue(sfMipModel), \"\\nOptimal x = \", getValue(x), \"\\nOptimal y = \", getValue(y))\n", + "end" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Min 0.3468 x[1] + 0.8687 x[2] + 0.12 x[3] + 0.5024 x[4] + 0.2884 x[5] + 0.2017 y[1] + 0.2712 y[2] + 0.4997 y[3] + 0.9238 y[4]\n", + "Subject to\n", + " 0.7511 x[1] - 0.1357 x[2] + 0.7955 x[3] - 0.4567 x[4] + 0.1356 x[5] - 0.0952 y[1] - 0.28056 y[2] - 1.33978 y[3] + 0.6506 y[4] == 0.1716\n", + " -0.667 x[1] - 0.3326 x[2] + 0.1657 x[3] - 0.5519 x[4] - 0.9367 x[5] - 0.8581 y[1] - 0.3518 y[2] + 1.2788 y[3] + 1.5114 y[4] == 0.361\n", + " 1.5894 x[1] - 0.1302 x[2] - 0.4313 x[3] - 0.4875 x[4] + 0.4179 x[5] - 0.5925 y[1] + 1.3477 y[2] + 0.1589 y[3] + 0.03495 y[4] == 0.0705\n", + " x[i] >= 0 for all i in {1,2..4,5}\n", + " y[i] >= 0, integer, for all i in {1,2,3,4}\n" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "\n" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "Status of the problem is Optimal\n", + "Optimal objective value = 1.070277955983598\n", + "Optimal x = " + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "x: 1 dimensions:\n", + "[1] = 0.0654906726037024\n", + "[2] = 0.0\n", + "[3] = 1.6298599862096397\n", + "[4] = 0.0\n", + "[5] = 1.221506908389311\n", + "\n", + "Optimal y = y: 1 dimensions:\n", + "[1] = 0.0\n", + "[2] = 0.0\n", + "[3] = 1.0\n", + "[4] = 0.0\n" + ] + } + ], + "prompt_number": 28 + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Reference\n", + "\n", + "[1] M. Lubin and I. Dunning, \u201cComputing in Operations Research using Julia\u201d, INFORMS Journal on Computing, to appear, 2014. [arXiv:1312.1431](http://arxiv.org/abs/1312.1431)" + ] + } + ], + "metadata": {} + } + ] +} \ No newline at end of file