{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# What is JuMP?\n", "\n", "JuMP is an _modeling language_ for optimization problems, writen in julia. \n", "\n", "When solving an optimization problem, typically you start with something like this:\n", "\n", "![alt text](img/pooling_problem.pdf)\n", "\n", "To solve this, you might use a _solver_: a software implementation of an optimization algorithm. They typically want the problem specified in a much more opaque way:\n", "\n", "java\n", "import ilog.concert.*;\n", "import ilog.cplex.*;\n", "public class Example {\n", " public static void main(String[] args) {\n", " try {\n", " IloCplex cplex = new IloCplex();\n", "double[] lb = {0.0, 0.0, 0.0};\n", "double[] ub = {40.0, Double.MAX_VALUE, Double.MAX_VALUE}; IloNumVar[] x = cplex.numVarArray(3, lb, ub);\n", " double[] objvals = {1.0, 2.0, 3.0};\n", "6\n", "cplex.addMaximize(cplex.scalProd(x, objvals));\n", " cplex.addLe(cplex.sum(cplex.prod(-1.0, x),\n", " cplex.prod( 1.0, x),\n", " cplex.prod( 1.0, x)), 20.0);\n", " cplex.addLe(cplex.sum(cplex.prod( 1.0, x),\n", " cplex.prod(-3.0, x),\n", " cplex.prod( 1.0, x)), 30.0);\n", "if ( cplex.solve() ) {\n", "cplex.out().println(\"Solution status = \" + cplex.getStatus()); cplex.out().println(\"Solution value = \" + cplex.getObjValue());\n", " double[] val = cplex.getValues(x);\n", " int ncols = cplex.getNcols();\n", " for (int j = 0; j < ncols; ++j)\n", " cplex.out().println(\"Column: \" + j + \" Value = \" + val[j]);\n", " }\n", " cplex.end();\n", " }\n", " catch (IloException e) {\n", " System.err.println(\"Concert exception '\" + e + \"' caught\");\n", "} }\n", "}\n", "\n", "\n", "For larger, more complex problem, programming like this is:\n", "* Time-consuming\n", "* Difficult\n", "* Hard to maintain/extend\n", "* Error-prone\n", "\n", "A modeling language (like JuMP) let's you code an optimization problem in a more natural way. It does the translation to the low-level solver format for you.\n", "\n", "There are a number of modeling languages out there. Why JuMP?\n", "\n", "* User-friendly\n", "* Matches performance of competitors\n", "* Solver-independent\n", "* Easy to extend and take advantage of advanced features\n", "\n", "In this session, we will focus on using JuMP for linear optimization. In the next class, we'll see how you can use if for discrete (i.e. integer) optimization, and for general nonlinear optimization.\n", "\n", "# Installing JuMP\n", "\n", "To install JuMP, just run" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[1m\u001b[34mINFO: Nothing to be done\n", "\u001b[0m\u001b[1m\u001b[34mINFO: METADATA is out-of-date — you may not have the latest version of JuMP\n", "\u001b[0m\u001b[1m\u001b[34mINFO: Use Pkg.update() to get the latest versions of your packages\n", "\u001b[0m" ] } ], "source": [ "Pkg.add(\"JuMP\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We already did this in the preassignment. To actually solve a problem, we will also need to install a solver package. There are 15+ options exposed in julia, each with support for different problem classes, different performance profiles, licensing requirements, etc. For the preassignment, we installed Gurobi, a best-of-breed linear/integer programming solver with a generous academic license." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[1m\u001b[34mINFO: Nothing to be done\n", "\u001b[0m\u001b[1m\u001b[34mINFO: METADATA is out-of-date — you may not have the latest version of Gurobi\n", "\u001b[0m\u001b[1m\u001b[34mINFO: Use Pkg.update() to get the latest versions of your packages\n", "\u001b[0m" ] } ], "source": [ "Pkg.add(\"Gurobi\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# A first example\n", "Let's see how we translate a simple, 2 variable LP to JuMP code.\n", "\n", "\n", "\\begin{align*}\n", "\\max_{x,y} \\quad& x + 2y \\\\\n", "\\text{s.t.}\\quad& x + y \\leq 1 \\\\\n", "& x, y \\geq 0.\n", "\\end{align*}\n", "\n", "\n", "First, we load the JuMP and Gurobi libraries." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [], "source": [ "using JuMP, Gurobi" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we construct a model object. This is a container for everything in our optimization problem: variables, constraints, solver options, etc." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "\\begin{alignat*}{1}\\min\\quad & 0\\\\\n", "\\text{Subject to} \\quad\\end{alignat*}\n", "" ], "text/plain": [ "Feasibility problem with:\n", " * 0 linear constraints\n", " * 0 variables\n", "Solver is Gurobi" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model = Model(solver=GurobiSolver())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we define the two decision variables in our optimization problem. We will use the @variable macro (a fancy function, essentially). The first argument is the model object to attach the variable to, and the second specifies the variable name and any bounds." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "$$y$$" ], "text/plain": [ "y" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@variable(model, x >= 0)\n", "@variable(model, y >= 0)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "\\begin{alignat*}{1}\\min\\quad & 0\\\\\n", "\\text{Subject to} \\quad & x \\geq 0\\\\\n", " & y \\geq 0\\\\\n", "\\end{alignat*}\n", "" ], "text/plain": [ "Feasibility problem with:\n", " * 0 linear constraints\n", " * 2 variables\n", "Solver is Gurobi" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now add the single constraint of our problem using the @constraint macro. We write it algebraically, exactly as we see it above." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "\\begin{alignat*}{1}\\min\\quad & 0\\\\\n", "\\text{Subject to} \\quad & x + y \\leq 1\\\\\n", " & x \\geq 0\\\\\n", " & y \\geq 0\\\\\n", "\\end{alignat*}\n", "" ], "text/plain": [ "Feasibility problem with:\n", " * 1 linear constraint\n", " * 2 variables\n", "Solver is Gurobi" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@constraint(model, x + y <= 1)\n", "model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We specify the objective function with the @objective macro." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "$$x + 2 y$$" ], "text/plain": [ "x + 2 y" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@objective(model, Max, x + 2y)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "\\begin{alignat*}{1}\\max\\quad & x + 2 y\\\\\n", "\\text{Subject to} \\quad & x + y \\leq 1\\\\\n", " & x \\geq 0\\\\\n", " & y \\geq 0\\\\\n", "\\end{alignat*}\n", "" ], "text/plain": [ "Maximization problem with:\n", " * 1 linear constraint\n", " * 2 variables\n", "Solver is Gurobi" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To solve the optimization problem, call the solve function." ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimize a model with 1 rows, 2 columns and 2 nonzeros\n", "Coefficient statistics:\n", " Matrix range [1e+00, 1e+00]\n", " Objective range [1e+00, 2e+00]\n", " Bounds range [0e+00, 0e+00]\n", " RHS range [1e+00, 1e+00]\n", "Presolve removed 1 rows and 2 columns\n", "Presolve time: 0.00s\n", "Presolve: All rows and columns removed\n", "Iteration Objective Primal Inf. Dual Inf. Time\n", " 0 2.0000000e+00 0.000000e+00 0.000000e+00 0s\n", "\n", "Solved in 0 iterations and 0.01 seconds\n", "Optimal objective 2.000000000e+00\n" ] }, { "data": { "text/plain": [ ":Optimal" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solve(model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now inspect the solution values and optimal cost." ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "getvalue(x) = 0.0\n", "getvalue(y) = 1.0\n", "getobjectivevalue(model) = 2.0\n" ] }, { "data": { "text/plain": [ "2.0" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@show getvalue(x)\n", "@show getvalue(y)\n", "@show getobjectivevalue(model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Exercise\n", "\n", "Code and solve the following optimization problem:\n", "\n", "\n", "\\begin{align*}\n", "\\min_{x,y} \\quad& 3x - y \\\\\n", "\\text{s.t.}\\quad& x + 2y \\geq 1 \\\\\n", "& x \\geq 0 \\\\\n", "& 0 \\leq y \\leq 1.\n", "\\end{align*}\n", "" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimize a model with 1 rows, 2 columns and 2 nonzeros\n", "Coefficient statistics:\n", " Matrix range [1e+00, 2e+00]\n", " Objective range [1e+00, 3e+00]\n", " Bounds range [1e+00, 1e+00]\n", " RHS range [1e+00, 1e+00]\n", "Presolve removed 1 rows and 2 columns\n", "Presolve time: 0.00s\n", "Presolve: All rows and columns removed\n", "Iteration Objective Primal Inf. Dual Inf. Time\n", " 0 -1.0000000e+00 0.000000e+00 0.000000e+00 0s\n", "\n", "Solved in 0 iterations and 0.00 seconds\n", "Optimal objective -1.000000000e+00\n" ] }, { "data": { "text/plain": [ ":Optimal" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "model = Model(solver=GurobiSolver())\n", "@variable(model, x >= 0)\n", "@variable(model, 0 <= y <= 1)\n", "@constraint(model, x + 2y >= 1)\n", "@objective(model, Min, 3x - y)\n", "\n", "stat = solve(model)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Airline Network Revenue Management\n", "\n", " \n", "\n", "In the airline network revenue management problem we are trying to decide how many tickets for each origin-destination (O-D) pair to sell at each price level. The goal is to maximize revenue, and we cannot sell more tickets than there is demand for, or space on the planes for.\n", "\n", "## Three Flight Problem\n", "\n", "We'll start with a toy problem that has three origin-destination pairs, with two price classes for each pair. The three origin-destination pairs are BOS-MDW, MDW-SFO, or BOS-SFO via MDW. BOS stands for Boston, MDW is Chicago-Midway, and SFO is San Francisco. Each O-D pair has a \"regular\" and \"discount\" fare class. The data we will use is summarized as follows:\n", "\n", "\n", "PLANE CAPACITY: 166\n", "\n", "BOS-MDW\n", " Regular Discount\n", "Price: 428 190\n", "Demand: 80 120\n", "\n", "BOS-SFO\n", " Regular Discount\n", "Price: 642 224\n", "Demand: 75 100\n", "\n", "MDW-SFO\n", " Regular Discount\n", "Price: 512 190\n", "Demand: 60 110\n", "" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "\\begin{alignat*}{1}\\min\\quad & 0\\\\\n", "\\text{Subject to} \\quad\\end{alignat*}\n", "" ], "text/plain": [ "Feasibility problem with:\n", " * 0 linear constraints\n", " * 0 variables\n", "Solver is Gurobi" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nrm = Model(solver=GurobiSolver())" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "\\begin{alignat*}{1}\\min\\quad & 0\\\\\n", "\\text{Subject to} \\quad & 0 \\leq BOStoMDW_R \\leq 80\\\\\n", " & 0 \\leq BOStoMDW_D \\leq 120\\\\\n", " & 0 \\leq BOStoSFO_R \\leq 75\\\\\n", " & 0 \\leq BOStoSFO_D \\leq 100\\\\\n", " & 0 \\leq MDWtoSFO_R \\leq 60\\\\\n", " & 0 \\leq MDWtoSFO_D \\leq 110\\\\\n", "\\end{alignat*}\n", "" ], "text/plain": [ "Feasibility problem with:\n", " * 0 linear constraints\n", " * 6 variables\n", "Solver is Gurobi" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@variables(nrm, begin \n", " 0 <= BOStoMDW_R <= 80\n", " 0 <= BOStoMDW_D <= 120\n", " 0 <= BOStoSFO_R <= 75\n", " 0 <= BOStoSFO_D <= 100\n", " 0 <= MDWtoSFO_R <= 60\n", " 0 <= MDWtoSFO_D <= 110\n", "end)\n", "nrm" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "\\begin{alignat*}{1}\\max\\quad & 428 BOStoMDW_R + 190 BOStoMDW_D + 642 BOStoSFO_R + 224 BOStoSFO_D + 512 MDWtoSFO_R + 190 MDWtoSFO_D\\\\\n", "\\text{Subject to} \\quad & 0 \\leq BOStoMDW_R \\leq 80\\\\\n", " & 0 \\leq BOStoMDW_D \\leq 120\\\\\n", " & 0 \\leq BOStoSFO_R \\leq 75\\\\\n", " & 0 \\leq BOStoSFO_D \\leq 100\\\\\n", " & 0 \\leq MDWtoSFO_R \\leq 60\\\\\n", " & 0 \\leq MDWtoSFO_D \\leq 110\\\\\n", "\\end{alignat*}\n", "" ], "text/plain": [ "Maximization problem with:\n", " * 0 linear constraints\n", " * 6 variables\n", "Solver is Gurobi" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@objective(nrm, Max, 428BOStoMDW_R + 190BOStoMDW_D +\n", " 642BOStoSFO_R + 224BOStoSFO_D +\n", " 512MDWtoSFO_R + 190MDWtoSFO_D)\n", "nrm" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "\\begin{alignat*}{1}\\max\\quad & 428 BOStoMDW_R + 190 BOStoMDW_D + 642 BOStoSFO_R + 224 BOStoSFO_D + 512 MDWtoSFO_R + 190 MDWtoSFO_D\\\\\n", "\\text{Subject to} \\quad & BOStoMDW_R + BOStoMDW_D + BOStoSFO_R + BOStoSFO_D \\leq 166\\\\\n", " & MDWtoSFO_R + MDWtoSFO_D + BOStoSFO_R + BOStoSFO_D \\leq 166\\\\\n", " & 0 \\leq BOStoMDW_R \\leq 80\\\\\n", " & 0 \\leq BOStoMDW_D \\leq 120\\\\\n", " & 0 \\leq BOStoSFO_R \\leq 75\\\\\n", " & 0 \\leq BOStoSFO_D \\leq 100\\\\\n", " & 0 \\leq MDWtoSFO_R \\leq 60\\\\\n", " & 0 \\leq MDWtoSFO_D \\leq 110\\\\\n", "\\end{alignat*}\n", "" ], "text/plain": [ "Maximization problem with:\n", " * 2 linear constraints\n", " * 6 variables\n", "Solver is Gurobi" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@constraint(nrm, BOStoMDW_R + BOStoMDW_D + \n", " BOStoSFO_R + BOStoSFO_D <= 166)\n", "@constraint(nrm, MDWtoSFO_R + MDWtoSFO_D + \n", " BOStoSFO_R + BOStoSFO_D <= 166)\n", "nrm" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimize a model with 2 rows, 6 columns and 8 nonzeros\n", "Coefficient statistics:\n", " Matrix range [1e+00, 1e+00]\n", " Objective range [2e+02, 6e+02]\n", " Bounds range [6e+01, 1e+02]\n", " RHS range [2e+02, 2e+02]\n", "Presolve time: 0.00s\n", "Presolved: 2 rows, 6 columns, 8 nonzeros\n", "\n", "Iteration Objective Primal Inf. Dual Inf. Time\n", " 0 1.7921000e+05 3.880000e+02 0.000000e+00 0s\n", " 3 1.2109000e+05 0.000000e+00 0.000000e+00 0s\n", "\n", "Solved in 3 iterations and 0.00 seconds\n", "Optimal objective 1.210900000e+05\n" ] }, { "data": { "text/plain": [ ":Optimal" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "status = solve(nrm)\n", "status # Should be :Optimal" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "getvalue(BOStoMDW_R) = 80.0\n", "getvalue(BOStoMDW_D) = 11.0\n", "getobjectivevalue(nrm) = 121090.0\n" ] }, { "data": { "text/plain": [ "121090.0" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@show getvalue(BOStoMDW_R)\n", "@show getvalue(BOStoMDW_D)\n", "@show getobjectivevalue(nrm)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "At this point, an exercise might be to extend this model by, say, adding another airport in this same fashion. I won't assign that, though, because it's a little tedious. It also doesn't scale well to problems with many airports. Instead, we can use JuMP's collections and summation notation to make compact, clear, and maintainable models for larger, more complex problems.\n", "\n", "First, we would like to construct a _collection of variables_ all at once. This is a very common idiom; for example, you might have a variable named x that is indexed from 1 to 10:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "$$x_{i} \\geq 0 \\quad\\forall i \\in \\{1,2,\\dots,9,10\\}$$" ], "text/plain": [ "10-element Array{JuMP.Variable,1}:\n", " x \n", " x \n", " x \n", " x \n", " x \n", " x \n", " x \n", " x \n", " x \n", " x" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m = Model()\n", "@variable(m, x[1:10] >= 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The index sets are specified inside the [...] block. You can create multidimensional containers by specifying multiple index sets, separated by commas:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "$$y_{i,j} \\leq 1 \\quad\\forall i \\in \\{1,2,\\dots,9,10\\}, j \\in \\{red,blue\\}$$" ], "text/plain": [ "y[i,j] ≤ 1 ∀ i ∈ {1,2,…,9,10}, j ∈ {red,blue}" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@variable(m, y[1:10,[\"red\",\"blue\"]] <= 1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For more complicated expressions, you can name the indices for the index sets and use them in the rest of the variable definition:" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "$$\\dots \\leq z_{i,j} \\leq \\dots \\quad\\forall i \\in \\{1,2,\\dots,9,10\\}, j \\in \\{\\dots\\}$$" ], "text/plain": [ "… ≤ z[i,j] ≤ … ∀ i ∈ {1,2,…,9,10}, j ∈ {…}" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ub = rand(10)\n", "@variable(m, i <= z[i=1:10,j=(i+1):10] <= ub[j])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To specify conditions on the indexing, you can add conditionals inside the [...] block, separated by a semicolon:" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "$$w_{i,c} \\geq 0 \\quad\\forall i \\in \\{1,2,\\dots,9,10\\}, c \\in \\{red,blue\\} s.t. iseven(i) || c == \"red\"$$" ], "text/plain": [ "w[i,c] ≥ 0 ∀ i ∈ {1,2,…,9,10}, c ∈ {red,blue} s.t. iseven(i) || c == \"red\"" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@variable(m, w[i=1:10, c=[\"red\",\"blue\"]; iseven(i) || c == \"red\"] >= 0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we can programatically create arrays of variables, we would like to be able to use them to full-effect in the constraints of our problem. That is, we want a way to express multi-dimensional summations, with conditionals. To do this, we use the sum(...) construction. The first argument is the ''inner loop'' of the summation, the index sets are specified after a for, and any conditionals are stated following an if (similar to variable definition, but with a slightly different syntax)." ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "$$x_{1} + x_{2} + x_{3} + x_{4} + x_{5} + x_{6} + x_{7} + x_{8} + x_{9} + x_{10} \\leq 1$$" ], "text/plain": [ "x + x + x + x + x + x + x + x + x + x ≤ 1" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@constraint(m, sum(x[i] for i in 1:10) <= 1)" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "$$2 y_{1,red} + 3 y_{1,blue} + 2 y_{2,red} + 3 y_{2,blue} + 2 y_{3,red} + 3 y_{3,blue} + 2 y_{4,red} + 3 y_{4,blue} + 2 y_{5,red} + 3 y_{5,blue} + 2 y_{6,red} + 3 y_{6,blue} + 2 y_{7,red} + 3 y_{7,blue} + 2 y_{8,red} + 3 y_{8,blue} + 2 y_{9,red} + 3 y_{9,blue} + 2 y_{10,red} + 3 y_{10,blue} = 1$$" ], "text/plain": [ "2 y[1,red] + 3 y[1,blue] + 2 y[2,red] + 3 y[2,blue] + 2 y[3,red] + 3 y[3,blue] + 2 y[4,red] + 3 y[4,blue] + 2 y[5,red] + 3 y[5,blue] + 2 y[6,red] + 3 y[6,blue] + 2 y[7,red] + 3 y[7,blue] + 2 y[8,red] + 3 y[8,blue] + 2 y[9,red] + 3 y[9,blue] + 2 y[10,red] + 3 y[10,blue] = 1" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "coef = Dict(\"red\" => 2, \"blue\" => 3)\n", "@constraint(m, sum(coef[c]*y[i,c] for i in 1:10, c in [\"red\",\"blue\"]) == 1)" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "$$2 z_{1,2} + 3 z_{1,3} + 4 z_{1,4} + 5 z_{1,5} + 6 z_{1,6} + 7 z_{1,7} + 8 z_{1,8} + 9 z_{1,9} + 10 z_{1,10} + 6 z_{2,3} + 8 z_{2,4} + 10 z_{2,5} + 12 z_{2,6} + 14 z_{2,7} + 16 z_{2,8} + 18 z_{2,9} + 20 z_{2,10} + 12 z_{3,4} + 15 z_{3,5} + 18 z_{3,6} + 21 z_{3,7} + 24 z_{3,8} + 27 z_{3,9} + 30 z_{3,10} + 20 z_{4,5} + 24 z_{4,6} + 28 z_{4,7} + 32 z_{4,8} + 36 z_{4,9} + 40 z_{4,10} + 30 z_{5,6} + 35 z_{5,7} + 40 z_{5,8} + 45 z_{5,9} + 50 z_{5,10} + 42 z_{6,7} + 48 z_{6,8} + 54 z_{6,9} + 60 z_{6,10} + 56 z_{7,8} + 63 z_{7,9} + 70 z_{7,10} + 72 z_{8,9} + 80 z_{8,10} + 90 z_{9,10} - w_{1,red} - 4 w_{2,red} - 4 w_{2,blue} - 9 w_{3,red} - 16 w_{4,red} - 16 w_{4,blue} - 25 w_{5,red} - 36 w_{6,red} - 36 w_{6,blue} - 49 w_{7,red} - 64 w_{8,red} - 64 w_{8,blue} - 81 w_{9,red} - 100 w_{10,red} - 100 w_{10,blue} \\leq 0$$" ], "text/plain": [ "2 z[1,2] + 3 z[1,3] + 4 z[1,4] + 5 z[1,5] + 6 z[1,6] + 7 z[1,7] + 8 z[1,8] + 9 z[1,9] + 10 z[1,10] + 6 z[2,3] + 8 z[2,4] + 10 z[2,5] + 12 z[2,6] + 14 z[2,7] + 16 z[2,8] + 18 z[2,9] + 20 z[2,10] + 12 z[3,4] + 15 z[3,5] + 18 z[3,6] + 21 z[3,7] + 24 z[3,8] + 27 z[3,9] + 30 z[3,10] + 20 z[4,5] + 24 z[4,6] + 28 z[4,7] + 32 z[4,8] + 36 z[4,9] + 40 z[4,10] + 30 z[5,6] + 35 z[5,7] + 40 z[5,8] + 45 z[5,9] + 50 z[5,10] + 42 z[6,7] + 48 z[6,8] + 54 z[6,9] + 60 z[6,10] + 56 z[7,8] + 63 z[7,9] + 70 z[7,10] + 72 z[8,9] + 80 z[8,10] + 90 z[9,10] - w[1,red] - 4 w[2,red] - 4 w[2,blue] - 9 w[3,red] - 16 w[4,red] - 16 w[4,blue] - 25 w[5,red] - 36 w[6,red] - 36 w[6,blue] - 49 w[7,red] - 64 w[8,red] - 64 w[8,blue] - 81 w[9,red] - 100 w[10,red] - 100 w[10,blue] ≤ 0" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@constraint(m, sum(i*j*z[i,j] for i in 1:10, j in (i+1):10) <=\n", " sum(i^2*w[i,c] for i in 1:10, c in [\"red\",\"blue\"] if iseven(i) || c == \"red\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's return to the network revenue management example and attempt to rewrite it in a generic way that scales to any number of airports. \n", "\n", "First, let's create some random data for our problem." ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "demand[(:BOS,:YYZ,:REG)] = 90\n" ] }, { "data": { "text/plain": [ "90" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Set the random seed to ensure we always\n", "# get the same stream of 'random' numbers\n", "srand(1988) \n", "\n", "# Lets create a vector of symbols, one for each airport\n", "airports = [:BOS, :MDW, :SFO, :YYZ]\n", "num_airport = length(airports)\n", "\n", "# We'll also create a vector of fare classes\n", "classes = [:REG, :DIS]\n", "\n", "# All the demand and price data for each triple of\n", "# (origin, destination, class) will be stored in\n", "# 'dictionaries', also known as 'maps'.\n", "demand = Dict()\n", "prices = Dict()\n", "\n", "# Generate a demand and price for each pair of airports\n", "# To keep the code simple we will generate info for\n", "# nonsense flights like BOS-BOS and SFO-SFO, but they\n", "# won't appear in our final model.\n", "for origin in airports, dest in airports\n", " # Generate demand:\n", " # - Regular demand is Uniform(50,90)\n", " # - Discount demand is Uniform(100,130)\n", " demand[(origin,dest,:REG)] = rand(50:90) \n", " demand[(origin,dest,:DIS)] = rand(100:130)\n", " # Generate prices:\n", " # - Regular price is Uniform(400,700)\n", " # - Discount price is Uniform(150,300)\n", " prices[(origin,dest,:REG)] = rand(400:700)\n", " prices[(origin,dest,:DIS)] = rand(150:300)\n", "end\n", "\n", "# Finally set all places to have the same capacity\n", "plane_cap = rand(150:200)\n", "\n", "# Lets look at a sample demand at random\n", "@show demand[(:BOS,:YYZ,:REG)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's build the model. We will have our decision variable x indexed by three things:\n", "\n", "1. Origin\n", "2. Destination\n", "3. Class\n", "\n", "The upper bound (the demand for each) will vary accordingly." ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "\\begin{alignat*}{1}\\min\\quad & 0\\\\\n", "\\text{Subject to} \\quad & 0 \\leq x_{BOS,MDW,REG} \\leq 55\\\\\n", " & 0 \\leq x_{BOS,MDW,DIS} \\leq 124\\\\\n", " & 0 \\leq x_{BOS,SFO,REG} \\leq 68\\\\\n", " & 0 \\leq x_{BOS,SFO,DIS} \\leq 126\\\\\n", " & 0 \\leq x_{BOS,YYZ,REG} \\leq 90\\\\\n", " & 0 \\leq x_{BOS,YYZ,DIS} \\leq 120\\\\\n", " & 0 \\leq x_{MDW,BOS,REG} \\leq 86\\\\\n", " & 0 \\leq x_{MDW,BOS,DIS} \\leq 117\\\\\n", " & 0 \\leq x_{MDW,SFO,REG} \\leq 75\\\\\n", " & 0 \\leq x_{MDW,SFO,DIS} \\leq 115\\\\\n", " & 0 \\leq x_{MDW,YYZ,REG} \\leq 63\\\\\n", " & 0 \\leq x_{MDW,YYZ,DIS} \\leq 100\\\\\n", " & 0 \\leq x_{SFO,BOS,REG} \\leq 67\\\\\n", " & 0 \\leq x_{SFO,BOS,DIS} \\leq 112\\\\\n", " & 0 \\leq x_{SFO,MDW,REG} \\leq 90\\\\\n", " & 0 \\leq x_{SFO,MDW,DIS} \\leq 120\\\\\n", " & 0 \\leq x_{SFO,YYZ,REG} \\leq 55\\\\\n", " & 0 \\leq x_{SFO,YYZ,DIS} \\leq 123\\\\\n", " & 0 \\leq x_{YYZ,BOS,REG} \\leq 53\\\\\n", " & 0 \\leq x_{YYZ,BOS,DIS} \\leq 115\\\\\n", " & 0 \\leq x_{YYZ,MDW,REG} \\leq 62\\\\\n", " & 0 \\leq x_{YYZ,MDW,DIS} \\leq 123\\\\\n", " & 0 \\leq x_{YYZ,SFO,REG} \\leq 61\\\\\n", " & 0 \\leq x_{YYZ,SFO,DIS} \\leq 118\\\\\n", "\\end{alignat*}\n", "" ], "text/plain": [ "Feasibility problem with:\n", " * 0 linear constraints\n", " * 24 variables\n", "Solver is Gurobi" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "nrm2 = Model(solver=GurobiSolver())\n", "\n", "@variable(nrm2, 0 <= x[o=airports,\n", " d=airports,\n", " c=classes; o!=d] <= demand[(o,d,c)])\n", "nrm2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The objective is to maximize the profit we make, summing over each ticket set:" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/latex": [ "\\begin{alignat*}{1}\\max\\quad & 551 x_{BOS,MDW,REG} + 243 x_{BOS,MDW,DIS} + 677 x_{BOS,SFO,REG} + 198 x_{BOS,SFO,DIS} + 697 x_{BOS,YYZ,REG} + 163 x_{BOS,YYZ,DIS} + 450 x_{MDW,BOS,REG} + 254 x_{MDW,BOS,DIS} + 534 x_{MDW,SFO,REG} + 237 x_{MDW,SFO,DIS} + 556 x_{MDW,YYZ,REG} + 158 x_{MDW,YYZ,DIS} + 587 x_{SFO,BOS,REG} + 256 x_{SFO,BOS,DIS} + 593 x_{SFO,MDW,REG} + 249 x_{SFO,MDW,DIS} + 674 x_{SFO,YYZ,REG} + 236 x_{SFO,YYZ,DIS} + 407 x_{YYZ,BOS,REG} + 157 x_{YYZ,BOS,DIS} + 647 x_{YYZ,MDW,REG} + 247 x_{YYZ,MDW,DIS} + 658 x_{YYZ,SFO,REG} + 285 x_{YYZ,SFO,DIS}\\\\\n", "\\text{Subject to} \\quad & 0 \\leq x_{BOS,MDW,REG} \\leq 55\\\\\n", " & 0 \\leq x_{BOS,MDW,DIS} \\leq 124\\\\\n", " & 0 \\leq x_{BOS,SFO,REG} \\leq 68\\\\\n", " & 0 \\leq x_{BOS,SFO,DIS} \\leq 126\\\\\n", " & 0 \\leq x_{BOS,YYZ,REG} \\leq 90\\\\\n", " & 0 \\leq x_{BOS,YYZ,DIS} \\leq 120\\\\\n", " & 0 \\leq x_{MDW,BOS,REG} \\leq 86\\\\\n", " & 0 \\leq x_{MDW,BOS,DIS} \\leq 117\\\\\n", " & 0 \\leq x_{MDW,SFO,REG} \\leq 75\\\\\n", " & 0 \\leq x_{MDW,SFO,DIS} \\leq 115\\\\\n", " & 0 \\leq x_{MDW,YYZ,REG} \\leq 63\\\\\n", " & 0 \\leq x_{MDW,YYZ,DIS} \\leq 100\\\\\n", " & 0 \\leq x_{SFO,BOS,REG} \\leq 67\\\\\n", " & 0 \\leq x_{SFO,BOS,DIS} \\leq 112\\\\\n", " & 0 \\leq x_{SFO,MDW,REG} \\leq 90\\\\\n", " & 0 \\leq x_{SFO,MDW,DIS} \\leq 120\\\\\n", " & 0 \\leq x_{SFO,YYZ,REG} \\leq 55\\\\\n", " & 0 \\leq x_{SFO,YYZ,DIS} \\leq 123\\\\\n", " & 0 \\leq x_{YYZ,BOS,REG} \\leq 53\\\\\n", " & 0 \\leq x_{YYZ,BOS,DIS} \\leq 115\\\\\n", " & 0 \\leq x_{YYZ,MDW,REG} \\leq 62\\\\\n", " & 0 \\leq x_{YYZ,MDW,DIS} \\leq 123\\\\\n", " & 0 \\leq x_{YYZ,SFO,REG} \\leq 61\\\\\n", " & 0 \\leq x_{YYZ,SFO,DIS} \\leq 118\\\\\n", "\\end{alignat*}\n", "" ], "text/plain": [ "Maximization problem with:\n", " * 0 linear constraints\n", " * 24 variables\n", "Solver is Gurobi" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "@objective(nrm2, Max, sum(prices[(o,d,c)]*x[o,d,c] for \n", " o in airports, d in airports, c in classes if o != d))\n", "nrm2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Our first set of constraints enforces that all the legs leaving the hub airport must not oversell the plane capacity:" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Adding constraint for hub (MDW) to BOS\n", "Adding constraint for hub (MDW) to SFO\n", "Adding constraint for hub (MDW) to YYZ\n" ] }, { "data": { "text/latex": [ "\\begin{alignat*}{1}\\max\\quad & 551 x_{BOS,MDW,REG} + 243 x_{BOS,MDW,DIS} + 677 x_{BOS,SFO,REG} + 198 x_{BOS,SFO,DIS} + 697 x_{BOS,YYZ,REG} + 163 x_{BOS,YYZ,DIS} + 450 x_{MDW,BOS,REG} + 254 x_{MDW,BOS,DIS} + 534 x_{MDW,SFO,REG} + 237 x_{MDW,SFO,DIS} + 556 x_{MDW,YYZ,REG} + 158 x_{MDW,YYZ,DIS} + 587 x_{SFO,BOS,REG} + 256 x_{SFO,BOS,DIS} + 593 x_{SFO,MDW,REG} + 249 x_{SFO,MDW,DIS} + 674 x_{SFO,YYZ,REG} + 236 x_{SFO,YYZ,DIS} + 407 x_{YYZ,BOS,REG} + 157 x_{YYZ,BOS,DIS} + 647 x_{YYZ,MDW,REG} + 247 x_{YYZ,MDW,DIS} + 658 x_{YYZ,SFO,REG} + 285 x_{YYZ,SFO,DIS}\\\\\n", "\\text{Subject to} \\quad & x_{MDW,BOS,REG} + x_{MDW,BOS,DIS} + x_{SFO,BOS,REG} + x_{SFO,BOS,DIS} + x_{YYZ,BOS,REG} + x_{YYZ,BOS,DIS} \\leq 182\\\\\n", " & x_{BOS,SFO,REG} + x_{BOS,SFO,DIS} + x_{MDW,SFO,REG} + x_{MDW,SFO,DIS} + x_{YYZ,SFO,REG} + x_{YYZ,SFO,DIS} \\leq 182\\\\\n", " & x_{BOS,YYZ,REG} + x_{BOS,YYZ,DIS} + x_{MDW,YYZ,REG} + x_{MDW,YYZ,DIS} + x_{SFO,YYZ,REG} + x_{SFO,YYZ,DIS} \\leq 182\\\\\n", " & 0 \\leq x_{BOS,MDW,REG} \\leq 55\\\\\n", " & 0 \\leq x_{BOS,MDW,DIS} \\leq 124\\\\\n", " & 0 \\leq x_{BOS,SFO,REG} \\leq 68\\\\\n", " & 0 \\leq x_{BOS,SFO,DIS} \\leq 126\\\\\n", " & 0 \\leq x_{BOS,YYZ,REG} \\leq 90\\\\\n", " & 0 \\leq x_{BOS,YYZ,DIS} \\leq 120\\\\\n", " & 0 \\leq x_{MDW,BOS,REG} \\leq 86\\\\\n", " & 0 \\leq x_{MDW,BOS,DIS} \\leq 117\\\\\n", " & 0 \\leq x_{MDW,SFO,REG} \\leq 75\\\\\n", " & 0 \\leq x_{MDW,SFO,DIS} \\leq 115\\\\\n", " & 0 \\leq x_{MDW,YYZ,REG} \\leq 63\\\\\n", " & 0 \\leq x_{MDW,YYZ,DIS} \\leq 100\\\\\n", " & 0 \\leq x_{SFO,BOS,REG} \\leq 67\\\\\n", " & 0 \\leq x_{SFO,BOS,DIS} \\leq 112\\\\\n", " & 0 \\leq x_{SFO,MDW,REG} \\leq 90\\\\\n", " & 0 \\leq x_{SFO,MDW,DIS} \\leq 120\\\\\n", " & 0 \\leq x_{SFO,YYZ,REG} \\leq 55\\\\\n", " & 0 \\leq x_{SFO,YYZ,DIS} \\leq 123\\\\\n", " & 0 \\leq x_{YYZ,BOS,REG} \\leq 53\\\\\n", " & 0 \\leq x_{YYZ,BOS,DIS} \\leq 115\\\\\n", " & 0 \\leq x_{YYZ,MDW,REG} \\leq 62\\\\\n", " & 0 \\leq x_{YYZ,MDW,DIS} \\leq 123\\\\\n", " & 0 \\leq x_{YYZ,SFO,REG} \\leq 61\\\\\n", " & 0 \\leq x_{YYZ,SFO,DIS} \\leq 118\\\\\n", "\\end{alignat*}\n", "" ], "text/plain": [ "Maximization problem with:\n", " * 3 linear constraints\n", " * 24 variables\n", "Solver is Gurobi" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "for d in airports\n", " if d != :MDW\n", " println(\"Adding constraint for hub (MDW) to $d\")\n", " @constraint(nrm2, \n", " sum(x[o,d,c] for o in airports, c in classes if o!=d) <= plane_cap)\n", " end\n", "end\n", "nrm2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, as an exercise, complete this model by adding constraints that each flight _to_ the hub is not oversold." ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Adding constraint for BOS to hub (MDW)\n", "Adding constraint for SFO to hub (MDW)\n", "Adding constraint for YYZ to hub (MDW)\n", "Optimize a model with 6 rows, 24 columns and 36 nonzeros\n", "Coefficient statistics:\n", " Matrix range [1e+00, 1e+00]\n", " Objective range [2e+02, 7e+02]\n", " Bounds range [5e+01, 1e+02]\n", " RHS range [2e+02, 2e+02]\n", "Presolve time: 0.00s\n", "Presolved: 6 rows, 24 columns, 36 nonzeros\n", "\n", "Iteration Objective Primal Inf. Dual Inf. Time\n", " 0 8.0150000e+05 2.254000e+03 0.000000e+00 0s\n", " 11 4.4785500e+05 0.000000e+00 0.000000e+00 0s\n", "\n", "Solved in 11 iterations and 0.00 seconds\n", "Optimal objective 4.478550000e+05\n", "getvalue(x) = x: 3 dimensions, 24 entries:\n", " [BOS,MDW,DIS] = 0.0\n", " [BOS,MDW,REG] = 55.0\n", " [BOS,SFO,DIS] = 0.0\n", " [BOS,SFO,REG] = 46.0\n", " [BOS,YYZ,DIS] = 0.0\n", " [BOS,YYZ,REG] = 81.0\n", " [MDW,BOS,DIS] = 42.0\n", " [MDW,BOS,REG] = 86.0\n", " [MDW,SFO,DIS] = 0.0\n", " [MDW,SFO,REG] = 75.0\n", " [MDW,YYZ,DIS] = 0.0\n", " [MDW,YYZ,REG] = 63.0\n", " [SFO,BOS,DIS] = 0.0\n", " [SFO,BOS,REG] = 54.0\n", " [SFO,MDW,DIS] = 0.0\n", " [SFO,MDW,REG] = 90.0\n", " [SFO,YYZ,DIS] = 0.0\n", " [SFO,YYZ,REG] = 38.0\n", " [YYZ,BOS,DIS] = 0.0\n", " [YYZ,BOS,REG] = 0.0\n", " [YYZ,MDW,DIS] = 59.0\n", " [YYZ,MDW,REG] = 62.0\n", " [YYZ,SFO,DIS] = 0.0\n", " [YYZ,SFO,REG] = 61.0\n", "getobjectivevalue(nrm2) = 447855.0\n" ] }, { "data": { "text/plain": [ "447855.0" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Constraints here!\n", "\n", "for o in airports\n", " if o != :MDW\n", " println(\"Adding constraint for$o to hub (MDW)\")\n", " @constraint(nrm2, \n", " sum(x[o,d,c] for d in airports, c in classes if o!=d) <= plane_cap)\n", " end\n", "end\n", " \n", "# @constraint(nrm2, constr[o=airports;o!=:MDW], \n", "# sum(x[o,d,c] for d in airports, c in classes if o!=d) <= plane_cap)\n", "nrm2\n", " \n", "# Now solve the model\n", "solve(nrm2)\n", "@show getvalue(x)\n", "@show getobjectivevalue(nrm2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Exercise: Chebyshev Center\n", "\n", "Take some polyhedron $P = \\{x \\in \\mathbb{R}^n : Ax \\leq b\\}$. Now take the largest ball\n", "$$\n", "B = \\{x + u : ||u||_2 \\leq r\\}\n", "$$\n", "given by center $x$ and radius $r$ such that $B \\subset P$. The center $x$ is the Chebyshev center of the polyhedron.\n", "\n", "We can find the Chebyshev center using linear optimization. We enforce that the ball $B$ lies inside the polyhedron $P$ by requiring that, for each constraint $i$ ($=1\\ldots,m$),\n", "$$\n", "A_i(x+u) \\leq b_i \\quad \\forall u : ||u||_2 \\leq r.\n", "$$\n", "Since $\\sup_{u : ||u||_2 \\leq r} A_i x = r||A_i||_2$, the constraint is equivalent to\n", "$$\n", "A_ix + r ||A_i||_2 \\leq b_i.\n", "$$\n", "Therefore, the problem of finding the Chebyshev center is equivalent to\n", "\n", "\\begin{align*}\n", "\\max_{x,r} \\quad& r \\\\\n", "\\text{s.t.} \\quad& A_i x + ||A_i||_2r \\leq b_i \\quad \\forall i = 1,\\ldots,m.\n", "\\end{align*}\n", "\n", "Your exercise is to code up this optimization problem and solve it, using the problem data I give you." ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimize a model with 4 rows, 3 columns and 12 nonzeros\n", "Coefficient statistics:\n", " Matrix range [1e+00, 2e+00]\n", " Objective range [1e+00, 1e+00]\n", " Bounds range [0e+00, 0e+00]\n", " RHS range [1e+00, 1e+00]\n", "Presolve removed 4 rows and 3 columns\n", "Presolve time: 0.00s\n", "Presolve: All rows and columns removed\n", "Iteration Objective Primal Inf. Dual Inf. Time\n", " 0 4.4721360e-01 0.000000e+00 0.000000e+00 0s\n", "\n", "Solved in 0 iterations and 0.00 seconds\n", "Optimal objective 4.472135955e-01\n", "getvalue(x) = [0.0,0.0]\n", "getvalue(r) = 0.4472135954999579\n" ] } ], "source": [ "n = 2\n", "m = 4\n", "# Store LHS as vector-of-vectors\n", "A = Vector{Float64}[\n", " [ 2, 1], [ 2,-1],\n", " [-1, 2], [-1,-2]]\n", "b = ones(m)\n", "\n", "# Build JuMP model\n", "model = Model(solver=GurobiSolver())\n", "\n", "# Your model goes here!\n", "# HINT: The dot(x,y) function will be useful for writing the constraints\n", "@variable(model, r)\n", "@variable(model, x[1:2])\n", "for i in 1:m\n", " @constraint(model, dot(A[i],x) + norm(A[i])*r <= b[i])\n", "end\n", "@objective(model, Max, r)\n", "\n", "# Now solve the model\n", "solve(model)\n", "@show getvalue(x)\n", "@show getvalue(r);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you solved the problem above correctly, you can visualize the solution using the code below." ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "collapsed": false }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "\u001b[1m\u001b[34mINFO: Nothing to be done\n", "\u001b[0m\u001b[1m\u001b[34mINFO: METADATA is out-of-date — you may not have the latest version of Gadfly\n", "\u001b[0m\u001b[1m\u001b[34mINFO: Use Pkg.update() to get the latest versions of your packages\n", "\u001b[0m" ] }, { "data": { "image/svg+xml": [ "\n", "\n" ], "text/html": [ "\n", "\n" ], "text/plain": [ "Plot(...)" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Use Gadfly to display the solution\n", "Pkg.add(\"Gadfly\")\n", "\n", "using Gadfly\n", "Gadfly.set_default_plot_size(8cm, 8cm)\n", "# Plot lines over [-1.5, 1.5]\n", "domain = linspace(-1, +1)\n", "# Plot circle across all angles\n", "θ = linspace(0,2π)\n", "plot(\n", "# a_1 x_1 + a_2 x_2 = b\n", "# --> x_2 = (b - a_1 x_1)/a_2\n", "[layer(x=domain,\n", " y=(b[i]-A[i]*domain)/A[i],\n", " Geom.line,\n", " Theme(line_width=2px,\n", " default_color=colorant\"blue\")) for i in 1:4]...,\n", "# x_1' = x_1 + rθ\n", "# x_2' = x_2 + rθ\n", "layer(x=getvalue(x) + getvalue(r)*cos(θ),\n", " y=getvalue(x) + getvalue(r)*sin(θ),\n", " Geom.path,\n", " Theme(line_width=5px,\n", " default_color=colorant\"red\")),\n", "Coord.Cartesian(ymin=-1,ymax=+1)\n", ")" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 0.5.0", "language": "julia", "name": "julia-0.5" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.5.0" } }, "nbformat": 4, "nbformat_minor": 2 }