{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "**Description**: Shows how to implement sensitivity analysis in JuMP.\n", "\n", "**Author**: Jack Dunn\n", "\n", "**License**: \"Creative
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Sensitivity Analysis using JuMP\n", "\n", "In this notebook, we will see how to use JuMP to conduct sensitivity analysis after solving a linear program and replicate the output offered by other packages (e.g. the Sensitivity Report in Excel Solver).\n", "\n", "We will consder a production planning problem by J E Beasley, made available on [OR-Notes](http://people.brunel.ac.uk/~mastjjb/jeb/or/lpsens_solver.html). \n", "\n", "### Problem Statement and Model Formulation\n", "\n", "(This section is reproduced exactly from the above link)\n", "\n", "A company manufactures four variants of the same product and in the final part of the manufacturing process there are assembly, polishing and packing operations. For each variant the time required for these operations is shown below (in minutes) as is the profit per unit sold.\n", "\n", "| Variant | Assembly | Polish | Pack | Profit |\n", "| :-----: | :------: | :----: | :--: | :----: |\n", "| 1 | 2 | 3 | 2 | 1.50 |\n", "| 2 | 4 | 2 | 3 | 2.50 |\n", "| 3 | 3 | 3 | 2 | 3.00 |\n", "| 4 | 7 | 4 | 5 | 4.50 |\n", "\n", "Given the current state of the labour force the company estimate that, each year, they have 100000 minutes of assembly time, 50000 minutes of polishing time and 60000 minutes of packing time available. How many of each variant should the company make per year and what is the associated profit?\n", "\n", "#### Variables\n", "\n", "Let $x_i \\geq 0$ be the number of units of variant i ($i=1,2,3,4$) made per year\n", "\n", "#### Constraints\n", "\n", "The operation time limits give the following constraints:\n", "\n", "- $2 x_1 + 4 x_2 + 3 x_3 + 7 x_4 \\leq 100000~~~$ (assembly) \n", "- $3 x_1 + 2 x_2 + 3 x_3 + 4 x_4 \\leq 50000~~~$ (polish) \n", "- $2 x_1 + 3 x_2 + 2 x_3 + 5 x_4 \\leq 60000~~~$ (pack)\n", "\n", "#### Objective\n", "\n", "Presumably to maximise profit - hence we have\n", "\n", "- $\\max 1.5 x_1 + 2.5 x_2 + 3.0 x_3 + 4.5 x_4$\n", "\n", "### Complete Formulation\n", "\n", "$$\\begin{align*}\n", "\\max~~~ &1.5 x_1 + 2.5 x_2 + 3.0 x_3 + 4.5 x_4 \\\\\n", "\\textrm{s.t.}~~~ &2 x_1 + 4 x_2 + 3 x_3 + 7 x_4 \\leq 100000 \\\\\n", " &3 x_1 + 2 x_2 + 3 x_3 + 4 x_4 \\leq 50000 \\\\\n", " &2 x_1 + 3 x_2 + 2 x_3 + 5 x_4 \\leq 60000 \\\\\n", " &x_i \\geq 0,~~~i = 1, 2, 3, 4\n", "\\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solving the model in JuMP\n", "\n", "Now we can formulate and solve this model in JuMP:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Optimize a model with 3 rows, 4 columns and 12 nonzeros\n", "Coefficient statistics:\n", " Matrix range [2e+00, 7e+00]\n", " Objective range [2e+00, 4e+00]\n", " Bounds range [0e+00, 0e+00]\n", " RHS range [5e+04, 1e+05]\n", "Presolve time: 0.00s\n", "Presolved: 3 rows, 4 columns, 12 nonzeros\n", "\n", "Iteration Objective Primal Inf. Dual Inf. Time\n", " 0 1.4000000e+31 7.875000e+30 1.400000e+01 0s\n", " 2 5.8000000e+04 0.000000e+00 0.000000e+00 0s\n", "\n", "Solved in 2 iterations and 0.00 seconds\n", "Optimal objective 5.800000000e+04\n", "getvalue(x) = [0.0,16000.000000000002,5999.999999999999,0.0]\n", "getobjectivevalue(model) = 58000.0\n" ] } ], "source": [ "# Define the data\n", "m = 3\n", "n = 4\n", "c = [1.5; 2.5; 3.0; 4.5]\n", "A = [2 4 3 7;\n", " 3 2 3 4;\n", " 2 3 2 5]\n", "b = [100000; 50000; 60000]\n", "\n", "# Import necessary packages and define model\n", "using JuMP\n", "using Gurobi # We need Gurobi for Sensitivity Analysis later\n", "model = Model(solver=GurobiSolver())\n", "\n", "# Define the variables\n", "@variable(model, x[i=1:n] >= 0)\n", "\n", "# Add the objective\n", "@objective(model, Max, sum{c[i] * x[i], i=1:n})\n", "\n", "# Add the constraints row-by-row, naming them according to each resource\n", "# The names will allow us to refer to the constraints during sensitivity analysis\n", "@constraint(model, assembly, dot(vec(A[1,:]), x) <= b[1])\n", "@constraint(model, polish, dot(vec(A[2,:]), x) <= b[2])\n", "@constraint(model, pack, dot(vec(A[3,:]), x) <= b[3])\n", "\n", "# Solve the model and show the optimal solution and objective value\n", "solve(model)\n", "@show getvalue(x)\n", "@show getobjectivevalue(model);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that the optimal production plan is to make 16000 units of variant 1 and 6000 units of variant 2, with an optimal profit of \\$58000" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Sensitivity Analysis\n", "\n", "Once we have solved a model, it is often useful to analyze the sensitivity of the solution to the model parameters. Other modeling tools like Excel Solver can produce a Sensitivity Report, which summarizes all of the sensitivity information in one table. \n", "\n", "The Sensitivity Report produced for the production planning solution above is as follows (image from OR-Tools):\n", "\n", "![Sensitivity Report](https://i.imgur.com/6tYVVwH.gif)\n", "\n", "The table contains information on the shadow prices and reduced costs in the model, as well as the ranges on the cost coefficients and right-hand side values for which the current basis is optimal.\n", "\n", "We will now reproduce this table using our JuMP model\n", "\n", "#### Final Values\n", "\n", "We can get the final values of the variables with `getValue()` as before:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "4-element Array{Float64,1}:\n", " 0.0\n", " 16000.0\n", " 6000.0\n", " 0.0" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x_final = getvalue(x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can get the final values of the constraints by calculating $Ax$:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 82000.0\n", " 50000.0\n", " 60000.0" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "con_final = A * x_final" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Reduced Costs/Shadow Prices\n", "\n", "We can extract the reduced costs by calling `getDual()` on the variables:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "4-element Array{Float64,1}:\n", " -1.5\n", " 0.0\n", " 0.0\n", " -0.2" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "red_costs = getdual(x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Similarly, we can extract the shadow prices by using `getDual()` on our constraint references:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.0" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "getdual(assembly)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.8" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "getdual(polish)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "0.3" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "getdual(pack)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can put these together into a vector as well:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3-element Array{Float64,1}:\n", " 0.0\n", " 0.8\n", " 0.3" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "shadow_prices = getdual([assembly; polish; pack])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Current Values\n", "\n", "The \"Objective Coefficient\" and \"Constraint R. H. Side\" columns simply contain the values of $c$ and $b$, respectively:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "4-element Array{Float64,1}:\n", " 1.5\n", " 2.5\n", " 3.0\n", " 4.5" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "obj_coeff = c" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3-element Array{Int64,1}:\n", " 100000\n", " 50000\n", " 60000" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "con_rhs = b" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Allowable Increase/Decrease\n", "\n", "We now want to retrieve the range of parameter values for which the optimal basis does not change. There are multiple ways to do this.\n", "\n", "One approach is to get the optimal basis using the `MathProgBase.getbasis()` function, and this can then be used to calculate the allowable ranges of increase and decrease using standard linear programming theory of sensitivity analysis.\n", "\n", "Alternatively, some solvers (like Gurobi) provide this information directly without the need for us to compute it manually. In this case, we can access the data through the Gurobi API directly, which is what we will do in the rest of this example. " ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "Gurobi Model: \n", " type : LP\n", " sense : maximize\n", " number of variables = 4\n", " number of linear constraints = 3\n", " number of quadratic constraints = 0\n", " number of sos constraints = 0\n", " number of non-zero coeffs = 12\n", " number of non-zero qp objective terms = 0\n", " number of non-zero qp constraint terms = 0\n" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Import Gurobi so we can access the API\n", "import Gurobi\n", "\n", "# Get a reference to the internal Gurobi model so we can use the API\n", "g = getrawsolver(model)" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "We can now access the sensitivity information directly. To do this, we make use of the `get_dblattrarray()` function from [Gurobi.jl](https://github.com/JuliaOpt/Gurobi.jl/blob/master/src/grb_attrs.jl#L58), which allows us to retrieve [attributes from the Gurobi API](http://www.gurobi.com/documentation/6.5/refman/attributes.html) that are arrays of floating point values. \n", "\n", "When calling `get_dblattrarray()`, we have to specify the internal Gurobi model, the name of the attribute we want to retrieve, the index at which to starting reading from the array, and the number of values to read.\n", "\n", "In our case, we use the `SARHSLow` and `SARHSUp` attributes to get the lower and upper bounds on the RHS values, and in each case we start at the first value and read in a total of $m$ values. Similarly, we use `SAObjLow` and `SAObjUp` to get the lower and upper bounds for the objective coefficients, and in this case we read in all $n$ values." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "rhs_lb = [82000.0,40000.0,33333.33333333333]\n", "rhs_ub = [1.0e100,90000.0,75000.0]\n", "obj_lb = [-1.0e100,2.357142857142857,2.499999999999999,-1.0e100]\n", "obj_ub = [3.0000000000000004,4.5,3.75,4.7]\n" ] } ], "source": [ "# RHS value lower and upper bounds\n", "rhs_lb = Gurobi.get_dblattrarray(g, \"SARHSLow\", 1, m)\n", "rhs_ub = Gurobi.get_dblattrarray(g, \"SARHSUp\", 1, m)\n", "@show rhs_lb\n", "@show rhs_ub\n", "\n", "# Objective coefficient lower and upper bounds\n", "obj_lb = Gurobi.get_dblattrarray(g, \"SAObjLow\", 1, n)\n", "obj_ub = Gurobi.get_dblattrarray(g, \"SAObjUp\", 1, n)\n", "@show obj_lb\n", "@show obj_ub;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The order of values in these arrays is not necessarily obvious in larger problems, and generally we do not know the order in which the information is returned by Gurobi. We can use the `getLinearIndex()` function on our variables and constraints to find their position in these arrays:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "4-element Array{Int64,1}:\n", " 1\n", " 2\n", " 3\n", " 4" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x_order = map(linearindex, x)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3-element Array{Int64,1}:\n", " 1\n", " 2\n", " 3" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "con_order = map(linearindex, [assembly, polish, pack])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that the variables and constraints are already ordered for us, but this isn't true all the time, so it pays to always rearrange the arrays according to this ordering in case the order is different" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [], "source": [ "rhs_lb_sorted = rhs_lb[con_order];\n", "rhs_ub_sorted = rhs_ub[con_order];\n", "obj_lb_sorted = obj_lb[x_order];\n", "obj_ub_sorted = obj_ub[x_order];" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, we can use these lower and upper bounds to obtain the allowable increase and decrease on each objective coefficient and RHS value:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "rhs_dec = con_rhs - rhs_lb_sorted = [18000.0,10000.0,26666.66666666667]\n", "rhs_inc = rhs_ub_sorted - con_rhs = [1.0e100,40000.0,15000.0]\n", "obj_dec = obj_coeff - obj_lb_sorted = [1.0e100,0.1428571428571428,0.5000000000000009,1.0e100]\n", "obj_inc = obj_ub_sorted - obj_coeff = [1.5000000000000004,2.0,0.75,0.20000000000000018]\n" ] } ], "source": [ "@show rhs_dec = con_rhs - rhs_lb_sorted;\n", "@show rhs_inc = rhs_ub_sorted - con_rhs;\n", "\n", "@show obj_dec = obj_coeff - obj_lb_sorted;\n", "@show obj_inc = obj_ub_sorted - obj_coeff;" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Final Sensitivity Table\n", "\n", "We can put all of this together to form the final Sensitivity Report tables.\n", "\n", "First, the variables:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "4x5 Array{Float64,2}:\n", " 0.0 -1.5 1.5 1.5 1.0e100 \n", " 16000.0 0.0 2.5 2.0 0.142857\n", " 6000.0 0.0 3.0 0.75 0.5 \n", " 0.0 -0.2 4.5 0.2 1.0e100 " ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "var_sensitivity = [x_final red_costs obj_coeff obj_inc obj_dec]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And similarly for the constraints:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "3x5 Array{Float64,2}:\n", " 82000.0 0.0 100000.0 1.0e100 18000.0\n", " 50000.0 0.8 50000.0 40000.0 10000.0\n", " 60000.0 0.3 60000.0 15000.0 26666.7" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "con_sensitivity = [con_final shadow_prices con_rhs rhs_inc rhs_dec]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This leads to the following tables, which we see are identical to those produced by Excel Solver:\n", "\n", "**Variables:**\n", "\n", "| Name | Final Value | Red. Cost | Obj. Coeff | Allow. Inc. | Allow. Dec. |\n", "| :---: | ----------: | --------: | ---------: | ----------: | ----------: |\n", "| $x_1$ | 0 | -1.5 | 1.5 | 1.50 | 1E+100 |\n", "| $x_2$ | 16000 | 0.0 | 2.5 | 2.00 | 0.1429 |\n", "| $x_3$ | 6000 | 0.0 | 3.0 | 0.75 | 0.5 |\n", "| $x_4$ | 0 | -0.2 | 4.5 | 0.20 | 1E+100 |\n", "\n", "**Constraints:**\n", "\n", "| Name | Final Value | Shadow Price | RHS | Allow. Inc. | Allow. Dec. |\n", "| :------: | ----------: | -----------: | -----: | ----------: | ----------: |\n", "| Assembly | 82000 | 0.0 | 100000 | 1E+100 | 18000 |\n", "| Polish | 50000 | 0.8 | 50000 | 40000 | 10000 |\n", "| Pack | 60000 | 0.3 | 60000 | 15000 | 26667 |\n" ] } ], "metadata": { "kernelspec": { "display_name": "Julia 0.4.3", "language": "julia", "name": "julia-0.4" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.4.3" } }, "nbformat": 4, "nbformat_minor": 0 }