{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# JuMP" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## New: Udated to the latest JuMP (0.20 at time of writing) ##" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The problem\n", "\n", "Given a single product, p plants and m markets, define the best route between p and m as to minimize transport costs:\n", "\n", "$ min \\sum_{p}{\\sum_{m} {c_{p,m} * x_{p,m}}} $\n", "\n", "s.t.\n", "\n", "$\\sum_{m} x_{p,m} <= a_p $\n", "\n", "$\\sum_{p} x_{p,m} >= b_m $\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Importing the libraries\n", "\n", "You will need to import as a minima the JuMP module. If you wish to specify a solver engine rather than letting JuMP select a suitable one, you will need to import also the module relative to the solver, e.g. Ipopt or GLPK" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Import of the JuMP, GLPK and CSV modules (the latter one just to import the data from a header based table, as in the original trasnport example in GAMS \n", "using JuMP, GLPK, CSV, DataFrames" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Defining the \"sets\"\n", "\n", "JuMP doesn't really have a concept of sets, but it uses the native containers available in the core Julia language\\\\Variables, parameters and constraints can be indexed using these containers.\n", "While many works with position-based lists, I find more readable using dictionaries instead. So the “sets” are represented as lists, but then everything else is a dictionary with the elements of the list as keys.\n", "One note: it seems that Julia/JuMP don't like much the “-” symbol, so I replaced it to “_”." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3-element Array{String,1}:\n", " \"new_york\"\n", " \"chicago\" \n", " \"topeka\" " ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Define sets #\n", "# Sets\n", "# i canning plants / seattle, san-diego /\n", "# j markets / new-york, chicago, topeka / ;\n", "plants = [\"seattle\",\"san_diego\"] # canning plants\n", "markets = [\"new_york\",\"chicago\",\"topeka\"] # markets" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Definition of the \"parameters\"\n", "\n", "Capacity of plants and demand of markets are directly defined as dictionaries, while the distance is first read as a DataFrame from a white-space separated table and then it is converted in a “(plant, market) ⇒ value” dictionary. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Dict{String,Int64} with 3 entries:\n", " \"chicago\" => 300\n", " \"new_york\" => 325\n", " \"topeka\" => 275" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Define parameters #\n", "# Parameters\n", "# a(i) capacity of plant i in cases\n", "# / seattle 350\n", "# san-diego 600 /\n", "a = Dict( # capacity of plant i in cases\n", " \"seattle\" => 350,\n", " \"san_diego\" => 600,\n", ")\n", " \n", "# b(j) demand at market j in cases\n", "# / new-york 325\n", "# chicago 300\n", "# topeka 275 / ;\n", "b = Dict( # demand at market j in cases\n", " \"new_york\" => 325,\n", " \"chicago\" => 300,\n", " \"topeka\" => 275,\n", ")\n", " " ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "

2 rows × 4 columns

plantsnew_yorkchicagotopeka
StringFloat64Float64Float64
1seattle2.51.71.8
2san_diego2.51.81.4
" ], "text/latex": [ "\\begin{tabular}{r|cccc}\n", "\t& plants & new\\_york & chicago & topeka\\\\\n", "\t\\hline\n", "\t& String & Float64 & Float64 & Float64\\\\\n", "\t\\hline\n", "\t1 & seattle & 2.5 & 1.7 & 1.8 \\\\\n", "\t2 & san\\_diego & 2.5 & 1.8 & 1.4 \\\\\n", "\\end{tabular}\n" ], "text/plain": [ "2×4 DataFrame\n", "│ Row │ plants │ new_york │ chicago │ topeka │\n", "│ │ \u001b[90mString\u001b[39m │ \u001b[90mFloat64\u001b[39m │ \u001b[90mFloat64\u001b[39m │ \u001b[90mFloat64\u001b[39m │\n", "├─────┼───────────┼──────────┼─────────┼─────────┤\n", "│ 1 │ seattle │ 2.5 │ 1.7 │ 1.8 │\n", "│ 2 │ san_diego │ 2.5 │ 1.8 │ 1.4 │" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Table d(i,j) distance in thousands of miles\n", "# new-york chicago topeka\n", "# seattle 2.5 1.7 1.8\n", "# san-diego 2.5 1.8 1.4 ;\n", "d_table = CSV.read(IOBuffer(\"\"\"\n", "plants new_york chicago topeka\n", "seattle 2.5 1.7 1.8\n", "san_diego 2.5 1.8 1.4\n", "\"\"\"), delim=\" \", ignorerepeated=true, copycols=true)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "2×3 Array{Float64,2}:\n", " 0.225 0.153 0.162\n", " 0.225 0.162 0.126" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "d = Dict( (r[:plants],m) => r[Symbol(m)] for r in eachrow(d_table), m in markets)\n", "# Here we are converting the table in a \"(plant, market) => distance\" dictionary\n", "# r[:plants]: the first key, row field using a fixed header\n", "# m: the second key\n", "# r[Symbol(m)]: the value, the row field with a dynamic header\n", " \n", "# Scalar f freight in dollars per case per thousand miles /90/ ;\n", "f = 90 # freight in dollars per case per thousand miles \n", " \n", "# Parameter c(i,j) transport cost in thousands of dollars per case ;\n", "# c(i,j) = f * d(i,j) / 1000 ;\n", "# We first declare an empty dictionary and then we fill it with the values\n", "c = Dict() # transport cost in thousands of dollars per case ;\n", "[ c[p,m] = f * d[p,m] / 1000 for p in plants, m in markets] " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Declaration of the model\n", "\n", "Here we declare a JuML optimisation model and we give it a name. This name will be then passed as first argument to all the subsequent operations, like creation of variables, constraints and objective function.\n", "We can, if we wish, works with several models at the same time.\n", "If we do not specify a solver, we let JuML use a suitable solver for the type of problem. Aside to specify the solver, we can also pass it solver-level options, e.g.: mymodel = Model(solver=IpoptSolver(print_level=0)) " ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$$ \\begin{alignat*}{1}\\text{feasibility}\\\\\n", "\\text{Subject to} \\quad\\end{alignat*}\n", " $$" ], "text/plain": [ "A JuMP Model\n", "Feasibility problem with:\n", "Variables: 0\n", "Model mode: AUTOMATIC\n", "CachingOptimizer state: EMPTY_OPTIMIZER\n", "Solver name: GLPK" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Model declaration (transport model)\n", "trmodel = Model(with_optimizer(GLPK.Optimizer,msg_lev=GLPK.MSG_ON))\n", "# or \n", "# trmodel = Model(with_optimizer(Ipopt.Optimizer,print_level=0)) # would require \"using Ipopt\" instead of GLPK" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Declaration of the model variables\n", "\n", "Variables can have multiple-dimensions - that is, being indexed under several indexes -, and bounds are given at the same time as their declaration.\n", "Differently from GAMS, we don't need to define the variable that is on the left hand side of the objective function. " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "## Define variables ##\n", "# Variables\n", "# x(i,j) shipment quantities in cases\n", "# z total transportation costs in thousands of dollars ;\n", "# Positive Variable x ;\n", "@variables trmodel begin\n", " x[p in plants, m in markets] >= 0 # shipment quantities in cases\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Declaration of the model constraints\n", "\n", "As in GAMS, each constraint can actually be a “family” of constraints: " ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "## Define contrains ##\n", "# supply(i) observe supply limit at plant i\n", "# supply(i) .. sum (j, x(i,j)) =l= a(i)\n", "# demand(j) satisfy demand at market j ; \n", "# demand(j) .. sum(i, x(i,j)) =g= b(j);\n", "@constraints trmodel begin\n", " supply[p in plants], # observe supply limit at plant p\n", " sum(x[p,m] for m in markets) <= a[p]\n", " demand[m in markets], # satisfy demand at market m\n", " sum(x[p,m] for p in plants) >= b[m]\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Declaration of the model objective\n", "\n", "Contrary to constraints and variables, the objective is always a unique function. Note that it is at this point that we specify the direction of the optimisation. " ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$$ 0.225 x_{seattle,new_york} + 0.225 x_{san_diego,new_york} + 0.153 x_{seattle,chicago} + 0.162 x_{san_diego,chicago} + 0.162 x_{seattle,topeka} + 0.12599999999999997 x_{san_diego,topeka} $$" ], "text/plain": [ "0.225 x[seattle,new_york] + 0.225 x[san_diego,new_york] + 0.153 x[seattle,chicago] + 0.162 x[san_diego,chicago] + 0.162 x[seattle,topeka] + 0.12599999999999997 x[san_diego,topeka]" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Objective\n", "@objective trmodel Min begin\n", " sum(c[p,m]*x[p,m] for p in plants, m in markets)\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Human-readable visualisation of the model (optional)\n", "\n", "If we wish, we can get the optimisation model printed in a human-readable fashion, so we can check that all is correct." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Min 0.225 x[seattle,new_york] + 0.225 x[san_diego,new_york] + 0.153 x[seattle,chicago] + 0.162 x[san_diego,chicago] + 0.162 x[seattle,topeka] + 0.12599999999999997 x[san_diego,topeka]\n", "Subject to\n", " x[seattle,new_york] + x[san_diego,new_york] ≥ 325.0\n", " x[seattle,chicago] + x[san_diego,chicago] ≥ 300.0\n", " x[seattle,topeka] + x[san_diego,topeka] ≥ 275.0\n", " x[seattle,new_york] + x[seattle,chicago] + x[seattle,topeka] ≤ 350.0\n", " x[san_diego,new_york] + x[san_diego,chicago] + x[san_diego,topeka] ≤ 600.0\n", " x[seattle,new_york] ≥ 0.0\n", " x[seattle,chicago] ≥ 0.0\n", " x[seattle,topeka] ≥ 0.0\n", " x[san_diego,new_york] ≥ 0.0\n", " x[san_diego,chicago] ≥ 0.0\n", " x[san_diego,topeka] ≥ 0.0\n" ] } ], "source": [ "print(trmodel)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Resolution of the model\n", "\n", "It is at this point that the solver is called and the model is passed to the solver engine for its solution. The return value is the status of the optimisation (“:Optimal” if all went fine) " ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 0: obj = 0.000000000e+00 inf = 9.000e+02 (3)\n", " 4: obj = 1.561500000e+02 inf = 0.000e+00 (0)\n", "* 5: obj = 1.536750000e+02 inf = 0.000e+00 (0)\n" ] } ], "source": [ "optimize!(trmodel)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Visualisation of the results\n", "\n", "While you can do any fancy output you may wish after you retrieve the optimal value of the variables with getvalue(var_name), you can just println(getvalue(x)) to get a basic output.\n", "Notice that you can also easily retrieve the dual value associated to the constraint with getdual(constraint_name). " ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Objective value: 153.675\n", "2-dimensional DenseAxisArray{Float64,2,...} with index sets:\n", " Dimension 1, [\"seattle\", \"san_diego\"]\n", " Dimension 2, [\"new_york\", \"chicago\", \"topeka\"]\n", "And data, a 2×3 Array{Float64,2}:\n", " 50.0 300.0 0.0\n", " 275.0 0.0 275.0\n", "Shadow prices of supply:\n", "seattle = 0.0\n", "san_diego = 0.0\n", "Shadow prices of demand:\n", "new_york = 0.225\n", "chicago = 0.153\n", "topeka = 0.12599999999999997\n" ] }, { "data": { "text/plain": [ "3-element Array{Nothing,1}:\n", " nothing\n", " nothing\n", " nothing" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "status = termination_status(trmodel)\n", "if (status == MOI.OPTIMAL) || (status == MOI.LOCALLY_SOLVED)\n", " println(\"Objective value: \",objective_value(trmodel))\n", " println(value.(x))\n", " println(\"Shadow prices of supply:\")\n", " [println(\"$p = $(dual(supply[p]))\") for p in plants]\n", " println(\"Shadow prices of demand:\")\n", " [println(\"$m = $(dual(demand[m]))\") for m in markets]\n", "else\n", " println(\"Model didn't solved\")\n", " println(status)\n", "end\n", "\n", "# Expected result:\n", "# obj= 153.675\n", "#['seattle','new-york'] = 50\n", "#['seattle','chicago'] = 300\n", "#['seattle','topeka'] = 0\n", "#['san-diego','new-york'] = 275\n", "#['san-diego','chicago'] = 0\n", "#['san-diego','topeka'] = 275" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "@webio": { "lastCommId": null, "lastKernelId": null }, "kernelspec": { "display_name": "Julia 1.2.0", "language": "julia", "name": "julia-1.2" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "1.2.0" } }, "nbformat": 4, "nbformat_minor": 1 }