{ "cells": [ { "cell_type": "markdown", "metadata": { "hide_cell": true }, "source": [ "Make me look good. Click on the cell below and press Ctrl+Enter." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "hide_cell": true }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.core.display import HTML\n", "HTML(open('css/custom.css', 'r').read())" ] }, { "cell_type": "markdown", "metadata": { "hide_cell": true }, "source": [ "
SM286D · Introduction to Applied Mathematics with Python · Spring 2020 · Uhan
\n", "\n", "
Lesson 16.
\n", "\n", "

Linear programming with Python

" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## This lesson..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Advanced methods for formulating linear programs with Python using __Pyomo__:\n", " - Using set-based notation\n", " - Using xlwings to get data and store results for an optimization model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Set-based notation in Pyomo" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- We learned the basics of using Pyomo to model and solve linear programs in Project 3. \n", "\n", "- Many of you have also used Pyomo in SA305, your course on linear programming.\n", "\n", "- In SA305, you may have learned about the benefits of using __symbolic input parameters__ or __set-based notation__ when formulating linear programs.\n", " - With set-based notation, we can formulate one model for all problems with the same structure.\n", " - As a result, we can easily reuse models and code to solve similar problems.\n", "\n", "- In this lesson, we'll go over an example of how to use __set-based notation__ in Pyomo." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- To illustrate the use of set-based notation in Pyomo, we'll use the following example taken from Chapter 3 of [*Pyomo Optimization Modeling in Python, 2nd Edition*](https://drive.google.com/open?id=1wKrm7UZ6an7g73Et6455gXgc1YCvGx2i), referred to as POMP below." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Example. (POMP Chapter 3)__ \n", "Consider the following warehouse location problem. \n", "Let $N$ be a set of candidate warehouse locations, and let $M$ be a set of customer locations. For\n", "each warehouse $n \\in N$, the cost of delivering product to customer $m \\in M$ is given by $d_{n,m}$.\n", "We want to determine the warehouse locations that will minimize the total\n", "cost of product delivery. \n", "\n", "Below is a formulation of this problem.\n", "The binary variables $y_n$ are used to define whether or not a\n", "warehouse should be built: $y_n$ is 1 if warehouse $n$ is built and 0 otherwise.\n", "The variables $x_{n,m}$ indicate the fraction of demand for customer $m$ that is served by\n", "warehouse $n$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\\begin{alignat}{3}\n", "\\min_{x,y} \\quad & \\sum_{n\\in N} \\sum_{m \\in M} d_{n,m} x_{n,m} &\\quad&&\\quad& (\\rm{WL.1})\\\\\n", "\\mbox{s.t.} \\quad & \\sum_{n\\in N} x_{n,m} = 1, &\\quad& \\forall m \\in M &\\quad& (\\rm{WL.2}) \\\\\n", " & x_{n,m} \\le y_n, &\\quad&\\forall n\\in N, m\\in M &\\quad& (\\rm{WL.3}) \\\\\n", " & \\sum_{n \\in N} y_n \\le P &\\quad&&\\quad& (\\rm{WL.4}) \\\\\n", " & 0 \\le x_{n,m} \\le 1 &\\quad&\\forall n \\in N, m \\in M&\\quad& (\\rm{WL.5}) \\\\\n", " & y_n \\in \\{0,1\\} &\\quad&\\forall n \\in N&\\quad& (\\rm{WL.6})\n", "\\end{alignat}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Quick check.__ Is this a linear program?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "_Write your answer here. Double-click to edit._\n", "\n", "- No! Although the objective function is linear, and constraints (WL.2) - (W.5) are also linear, the variables $y_n$ are required to be binary. This requirement makes the above formulation an _integer_ linear program. You'll learn more about these in SA405.\n", "\n", "- Nevertheless, we can still use this model to learn about how to use set-based notation in Pyomo." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- In the code cell below, we begin by importing `pyomo.environ` as `pyo`, defining all the required data for the example, and defining a concrete model.\n", "\n", " - Note the use of __tuples__ as keys for the dictionary `d`. For example, in the code below, we have\n", " \n", " ```python\n", " d[('Harlingen', 'NYC')] = 1956\n", " ```\n", " The idea is that `d[(n, m)]` represents $d_{n,m}$ in the formulation above.\n", " \n", " - In Python, a tuple is a collection of items in round parentheses `()`, separated by commas.\n", " \n", " - For our purposes here, you can think of a tuple as a list-like object that you can use as a dictionary key. \n", " - Try using a list as a dictionary key. It won't work... 🤔\n", " - For more details on lists vs. tuples, take a look at [this article from Real Python](https://realpython.com/python-lists-tuples/)." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# First import Pyomo\n", "import pyomo.environ as pyo\n", "\n", "# Define data\n", "N = ['Harlingen', 'Memphis', 'Ashland']\n", "M = ['NYC', 'LA', 'Chicago', 'Houston']\n", "d = {\n", " ('Harlingen', 'NYC'): 1956,\n", " ('Harlingen', 'LA'): 1606,\n", " ('Harlingen', 'Chicago'): 1410,\n", " ('Harlingen', 'Houston'): 330,\n", " ('Memphis', 'NYC'): 1096,\n", " ('Memphis', 'LA'): 1792,\n", " ('Memphis', 'Chicago'): 531,\n", " ('Memphis', 'Houston'): 567,\n", " ('Ashland', 'NYC'): 485,\n", " ('Ashland', 'LA'): 2322,\n", " ('Ashland', 'Chicago'): 324,\n", " ('Ashland', 'Houston'): 1236 \n", "}\n", "P = 2\n", "\n", "# Define the model as a Concrete Model in Pyomo\n", "model = pyo.ConcreteModel()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Next, we define the sets $N$ and $M$ in Pyomo. Note that they are initialized with the lists `N` and `M` defined above." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# Define sets \n", "model.N = pyo.Set(initialize=N, doc='Set of candidate warehouse locations')\n", "model.M = pyo.Set(initialize=M, doc='Set of customer locations')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Below, we use the sets $N$ and $M$ to define the parameters $d_{n,m}$ for $n \\in N$ and $m \\in M$. Note that it is initialized using the dictionary `d` defined in the code above." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# Define parameters\n", "model.d = pyo.Param(model.N, model.M, initialize=d, doc='cost of delivering product to customer m from warehouse n')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Now, we define the decision variables $x_{n,m}$ for $n \\in N$ and $m \\in M$, and $y_n$ for $n \\in N$. Note that \n", " - we use the keyword argument `bounds` to restrict $x_{n,m}$ to values between 0 and 1, and\n", " \n", " - we use the keyword argument `within` to specify the domain of the decision variable $y_n$ to be `pyo.Binary`." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# Define decision variables\n", "model.x = pyo.Var(model.N, model.M, bounds=(0,1), doc='fraction of demand for customer m served by warehouse n')\n", "model.y = pyo.Var(model.N, within=pyo.Binary, doc='1 if warehouse n is open, 0 otherwise')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Next we use _construction rules_ to define the objective function and constraints. \n", " \n", "- We can add the objective function (WL.1)\n", "\n", " $$\n", " \\min_{x,y} \\quad \\sum_{n\\in N} \\sum_{m \\in M} d_{n,m} x_{n,m} \\qquad (\\rm{WL.1})\n", " $$\n", " \n", " like this:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# Define objective (WL.1)\n", "def obj_rule(model):\n", " return sum(model.d[n,m] * model.x[n,m] for n in model.N for m in model.M)\n", "model.obj = pyo.Objective(rule=obj_rule, sense=pyo.minimize)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Next, we can add the constraint (WL.2) that ensures each customer's demand is exactly satisfied\n", "\n", " $$\n", " \\sum_{n\\in N} x_{n,m} = 1 \\quad \\forall m \\in M \\qquad (\\rm{WL.2})\n", " $$\n", " \n", " like this:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# Define \"100% of each customer's demand satisfied\" constraint (WL.2)\n", "def one_per_cust_rule(model, m):\n", " return sum(model.x[n,m] for n in model.N) == 1 \n", "model.one_per_cust = pyo.Constraint(model.M, rule=one_per_cust_rule, doc=\"Each customer's demand is fully met\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- We can add constraints (WL.3) and (WL.4) in a similar way:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# Define \"warehouse is active\" constraint (WL.3)\n", "def warehouse_active_rule(model, n, m):\n", " return model.x[n,m] <= model.y[n] \n", "model.warehouse_active = pyo.Constraint(model.N, model.M, rule=warehouse_active_rule, doc=\"Warehouse can only deliver if it is built\")\n", "\n", "# Define warehouse capacity constraint (WL.4)\n", "def num_warehouses_rule(model):\n", " return sum(model.y[n] for n in model.N) <= P \n", "model.num_warehouses = pyo.Constraint(rule=num_warehouses_rule, doc=\"Can open at most P warehouses\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Now that we've finished writing the model in Pyomo, we need to solve it. \n", "\n", " - The code below sets the solver to `\"glpk\"`, which corresponds to GLPK, the open source GNU Linear Programming Kit solver. \n", "\n", " - Then the code solves `model` and saves the results object in `results`.\n", "\n", " - Note that in this example, we set `tee=False`, which means we will *not* see all the output from the solver." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "# Set the solver to be GLPK\n", "opt = pyo.SolverFactory(\"glpk\")\n", "\n", "# Call the solver\n", "# Note that we set tee = False, \n", "# so we don't see all the output from the solver\n", "results = opt.solve(model, tee=False)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Finally, we produce output to show the results from the solver.\n", " - First, we check to see if the solver termination condition is optimal. \n", " - If it is, we print out the value of the objective function and other important decisions made by the model." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Results\n", "\n", "Minimum total cost is 2745.00.\n", "\n", "Open Warehouses\n", "---------------\n", "Harlingen \n", "Ashland \n", "\n", "\n", "Customer Warehouse Fraction of Demand Satistifed\n", "---------- ---------- -----------------------------\n", "LA Harlingen 1.00\n", "Houston Harlingen 1.00\n", "NYC Ashland 1.00\n", "Chicago Ashland 1.00\n" ] } ], "source": [ "# Check to see if the solver termination condition was optimal\n", "# If it's optimal, print the optimal objective value and solution\n", "if str(results.solver.termination_condition) == \"optimal\":\n", " print(\"Results\\n\")\n", " \n", " # Print the objective function value\n", " print(f\"Minimum total cost is {pyo.value(model.obj):0.2f}.\\n\")\n", " print(\"Open Warehouses\")\n", " print(\"---------------\")\n", " \n", " # Below we loop over all the elements of the set N\n", " for n in model.N:\n", " \n", " # Here we check to see if the decision variable value is \n", " # not equal to 0 for this value of n\n", " if model.y[n].value != 0:\n", " # If the decision variable is not equal to 0, we\n", " # print the value of n. This will tell us which\n", " # warehouses the model decided to open.\n", " print(f\"{n:15s}\")\n", " \n", " print(\"\\n\")\n", " print(\"Customer Warehouse Fraction of Demand Satistifed\")\n", " print(\"---------- ---------- -----------------------------\")\n", " \n", " # Below we loop over all the elements of the sets N and M\n", " for n in model.N:\n", " for m in model.M:\n", " \n", " # Here we check to see if the decision variable value is\n", " # not equal to 0 for this combination of n and m.\n", " if model.x[n,m].value != 0:\n", " # If the decision variable is not equal to 0, we\n", " # print the values m, n, and the decision variable.\n", " # This will tell us how much demand was satisfied for \n", " # each customer from each servicing warehouse.\n", " print(f\"{m:10s} {n:10s} {model.x[n,m].value:<.2f}\")\n", "\n", "# Otherwise, print a warning message\n", "else:\n", " print(\"Solver termination condition was not optimal.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Notes and additional reading\n", "\n", "- Read pages 29-38 from POMP for an even more detailed treatment of this example.\n", "\n", " - If you read nothing else from POMP, read Section 3.3.3 of POMP, where the authors discuss the use of construction rules.\n", "\n", "- There are a few minor differences in the code above over what is presented in the book. \n", " - The changes are primarily based on the decision in the code below to use `pyo` as an alias when importing `pyomo.environ`. This choice makes the model slightly more complicated, but it prevents pollution of the name space and is therefore recommended as good Python practice. \n", " - In addition, the authors choose to work with $N$ and $M$ as Python lists and $d$ as a Python dictionary. Above, we define $N$ and $M$ as Pyomo sets and $d$ as a Pyomo parameter. These changes are not required, but can have advantages when you want the model to be more \"self-contained\". Both the authors and the code above work with $P$ as a Python variable. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Classwork" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Based on the example above, you should be able to solve Problems 1 and 2 below. \n", "\n", "- In Problem 2, you will grab the necessary data from an Excel workbook. You'll need to use of `xlwings`, so you may want to refer back to Lesson 10 on working with spreadsheets.\n", "\n", "- Before proceeding with Problems 1 and 2, __do Problem 0 first__." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Problem 0.__ Read the sections on \"Line Continuation\", \"Implicit Line Continuation\", and \"Explicit Line Continuation\" in [this article on Real Python](https://realpython.com/python-program-structure/).\n", "\n", "\"Line continuation\" might sound scary, but it's just about syntax. In particular, what if you have a line of Python code that is really, really long - so long that it can't all be displayed at once? This article outlines some techniques you can write your code to avoid this problem and make your code more readable.\n", "\n", "You'll find that the solutions to this lesson in particular, uses some of these techniques. You'll also find that we've already been using some of these techniques stealthily in previous lessons." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Problem 1. (Winston, Operations Research: Applications and Algorithms, Fourth Edition, Exercise 8 on page 93)__ Highland's TV-Radio Store must determine how many TVs and radios to keep in stock. A TV requires 10 sq. ft. of floorspace, whereas a radio requires 4 sq. ft. 200 sq. ft. of floorspace is available. A TV will earn Highland \\\\$60 in profits and a radio will earn \\\\$20. The store stocks only TVs and radios. Marketing requirements dictate that at least 60% of all appliances in stock be radios. TVs tie up \\\\$200 in capital and radios tie up \\\\$50. Highland wants to have at most \\\\$3000 worth of capital tied up at any time. \n", "\n", "Below is a linear program without set-based notation that can be used to maximize Highland's profit. \n", "\n", "__Linear program without set-based notation.__\n", "\n", "_Decision variables._ Let $t$ equal the number of TVs to keep in stock and let $r$ equal the number of radios to keep in stock.\n", "\n", "_Formulation._\n", "\\begin{alignat}{3}\n", "\\max \\quad & 60t + 20r \\\\\n", "\\mbox{s.t.} \\quad & 10t + 4r \\le 200 &\\qquad& (\\rm{floorspace}) &\\qquad& \\\\\n", " & r \\ge 0.60(t+r) & & (\\rm{marketing}) \\\\\n", " & 200t + 50r \\le 3000 && (\\rm{capital})\\\\\n", " & t,r \\ge 0 && \n", "\\end{alignat}\n", "\n", "Note that we can rewrite the marketing constraint as\n", "\n", "\\begin{alignat}{3}\n", " & 0.60 t - 0.40 r \\le 0 &\\qquad & (\\rm{ marketing}) \n", "\\end{alignat}\n", "\n", "Here is the same linear program, but with set-based notation. Make sure you understand how these two formulations are related.\n", "\n", "__Linear program with set-based notation.__\n", "\n", "_Indices and sets._\n", "\\begin{array}{ll}\n", " \\mbox{$p \\in P$} & \\mbox{products}, P=\\{\\mbox{tv, radio} \\} \\\\\n", "\\end{array}\n", "\n", "_Data._\n", "\\begin{array}{ll}\n", " \\mbox{profit}_p & \\mbox{profit earned for product $p$} \\\\\n", " \\mbox{floorspace}_p & \\mbox{sq. ft. of floorspace required for product $p$} \\\\\n", " \\mbox{marketing}_p & \\mbox{marketing percentage for product $p$} \\\\\n", " \\mbox{capital}_p & \\mbox{capital dollars tied up for product $p$} \\\\\n", " \\mbox{floorspace_avail} & \\mbox{total sq of floorspace available} \\\\\n", " \\mbox{capital_avail} & \\mbox{total dollars of capital available} \\\\\n", "\\end{array}\n", "\n", "_Decision variables._\n", "\\begin{array}{ll}\n", " \\mbox{PRODUCTS}_p & \\mbox{number of items of product $p$ to put on the floor} \\\\\n", "\\end{array}\n", "\n", "_Formulation._\n", "\\begin{alignat}{3}\n", "\\max \\quad & \\sum_{p\\in P} \\mbox{profit}_p * \\mbox{PRODUCTS}_p \\\\\n", "\\mbox{s.t.} \\quad & \\sum_{p\\in P} \\mbox{floorspace}_p * \\mbox{PRODUCTS}_p \\le \\mbox{floorspace_avail} &\\qquad& &\\qquad& \\\\\n", " & \\sum_{p \\in P} \\mbox{marketing}_p * \\mbox{PRODUCTS}_p \\le 0 && \\\\\n", " & \\sum_{p \\in P} \\mbox{capital}_p * \\mbox{PRODUCTS}_p \\le \\mbox{capital_avail} && \\\\\n", " & \\mbox{PRODUCTS}_p \\ge 0 && \\forall p\\in P\n", "\\end{alignat}\n", "\n", "
\n", "\n", "- Implement the linear program __with set-based notation__ in Pyomo, and solve it.\n", " \n", " - Note that the decision variables are _not_ constrained to be integer, so feasible solutions can have fractional variables. Nevertheless, the model is still useful: the same formulation, but with integer decision variables is solved on a smaller feasible set. So, the linear program gives an upper bound on the maximum possible profit.\n", "\n", "- Verify that both $(t,r) = (6,35)$ and $(t,r) = (7,32)$ are feasible solutions. Which of these gives the higher profit?\n", "\n", "- Find how close the feasible integer solutions above are to the upper bound on the maximum possible profit given by the linear program's optimal value. You should report your answer as a percentage, e.g. \n", "\n", " ```\n", " The answer with t = ... and r = ... is within 92.6% of the upper bound on the maximum profit.\n", " ```" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Results\n", "\n", "Maximum Total Profit is 1066.67 dollars.\n", "\n", "Product Number Proft Made \n", "------------- ----------- -----------\n", "tv 6.67 400.00\n", "radio 33.33 666.67\n" ] } ], "source": [ "# Import Pyomo\n", "import pyomo.environ as pyo\n", "\n", "# Define data\n", "productslist = ['tv', 'radio']\n", "profitdict = {'tv': 60, 'radio': 20}\n", "floorspacedict = {'tv': 10, 'radio': 4}\n", "marketingdict = {'tv': 0.6, 'radio': -0.4}\n", "capitaldict = {'tv': 200, 'radio': 50}\n", "floorspace_avail = 200\n", "capital_avail = 3000\n", "\n", "# Define the model as a Concrete Model in Pyomo\n", "model = pyo.ConcreteModel()\n", "\n", "# Define sets \n", "model.P = pyo.Set(initialize=productslist, doc='Set of Products')\n", "\n", "# Define parameters\n", "model.profit = pyo.Param(model.P, initialize=profitdict, \n", " doc='profit for item p')\n", "model.floorspace = pyo.Param(model.P, initialize=floorspacedict, \n", " doc='floorspace for item p')\n", "model.marketing = pyo.Param(model.P, initialize=marketingdict, \n", " doc='marketing for item p')\n", "model.capital = pyo.Param(model.P, initialize=capitaldict, \n", " doc='capital for item p')\n", "\n", "# Define decision variables\n", "model.PRODUCTS = pyo.Var(model.P, within=pyo.NonNegativeReals, \n", " doc='number of items to put on floor')\n", "\n", "# Define objective\n", "# Note that the default objective sense is minimize in Pyomo\n", "def profit_rule(model):\n", " return sum(model.profit[p] * model.PRODUCTS[p] for p in model.P)\n", "model.cost = pyo.Objective(rule=profit_rule, sense=pyo.maximize)\n", "\n", "# Define constraints\n", "def floorspace_rule(model):\n", " return (\n", " sum(model.PRODUCTS[p] * model.floorspace[p] for p in model.P) \n", " <= floorspace_avail \n", " )\n", "model.avail_floorspace_rule = pyo.Constraint(rule=floorspace_rule, \n", " doc='Floor space constraint')\n", "\n", "def marketing_rule(model):\n", " return (\n", " sum(model.PRODUCTS[p] * model.marketing[p] for p in model.P)\n", " <= 0 \n", " )\n", "model.avail_marketing_rule = pyo.Constraint(rule=marketing_rule, \n", " doc='Marketing constraint')\n", "\n", "def capital_rule(model):\n", " return (\n", " sum(model.PRODUCTS[p] * model.capital[p] for p in model.P) \n", " <= capital_avail \n", " )\n", "model.avail_capital_rule = pyo.Constraint(rule=capital_rule, \n", " doc='Capital constraint')\n", "\n", "# Set the solver to be GLPK\n", "opt = pyo.SolverFactory(\"glpk\")\n", "\n", "# Call the solver\n", "results = opt.solve(model, tee=False)\n", "\n", "# If it's optimal, print the optimal objective value and solution\n", "if str(results.solver.termination_condition) == \"optimal\":\n", " print(\"Results\\n\")\n", " print(f\"Maximum Total Profit is {pyo.value(model.cost):0.2f} dollars.\\n\")\n", " print(\"Product Number Proft Made \")\n", " print(\"------------- ----------- -----------\")\n", " for p in model.P:\n", " if model.PRODUCTS[p].value != 0:\n", " # Note that multiple strings (plain or f-string) in a print statement \n", " # (with nothing in between) will be concatenated\n", " print(f\"{p:13s} {model.PRODUCTS[p].value:11.2f} \"\n", " f\"{model.PRODUCTS[p].value*model.profit[p]:11.2f}\")" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Stocking 7 televisions and 32 radios is feasible and gives profit equal to $1060.\n", "Stocking 6 televisions and 35 radios is feasible and gives profit equal to $1060.\n", "Both stocking plans are within just $6.67 of the LP relaxation.\n", "Either integer solution is at least 99.37% optimal.\n" ] } ], "source": [ "# Note that multiple strings (plain or f-string) in a print statement \n", "# (with nothing in between) will be concatenated\n", "print(f'Stocking 7 televisions and 32 radios is feasible and gives '\n", " f'profit equal to ${str(7 * 60 + 32 * 20)}.')\n", "print(f'Stocking 6 televisions and 35 radios is feasible and gives '\n", " f'profit equal to ${str(6 * 60 + 35 * 20)}.')\n", "print(f'Both stocking plans are within just $6.67 of the LP relaxation.')\n", "print(f'Either integer solution is at least {100 * 1060 / 1066.67:.2f}% optimal.')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If we wanted to ensure that the number of each product is integer, we can change the type of decision variable on line 30 below.\n", "\n", "While the change in code is simple, how the resulting model is solved is _very_ different from the one above. You'll learn more about this in SA405." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Results\n", "\n", "Maximum Total Profit is 1060.00 dollars.\n", "\n", "Product Number Proft Made \n", "------------- ----------- -----------\n", "tv 6.00 360.00\n", "radio 35.00 700.00\n" ] } ], "source": [ "# Import Pyomo\n", "import pyomo.environ as pyo\n", "\n", "# Define data\n", "productslist = ['tv', 'radio']\n", "profitdict = {'tv': 60, 'radio': 20}\n", "floorspacedict = {'tv': 10, 'radio': 4}\n", "marketingdict = {'tv': 0.6, 'radio': -0.4}\n", "capitaldict = {'tv': 200, 'radio': 50}\n", "floorspace_avail = 200\n", "capital_avail = 3000\n", "\n", "# Define the model as a Concrete Model in Pyomo\n", "model = pyo.ConcreteModel()\n", "\n", "# Define sets \n", "model.P = pyo.Set(initialize=productslist, doc='Set of Products')\n", "\n", "# Define parameters\n", "model.profit = pyo.Param(model.P, initialize=profitdict, \n", " doc='profit for item p')\n", "model.floorspace = pyo.Param(model.P, initialize=floorspacedict, \n", " doc='floorspace for item p')\n", "model.marketing = pyo.Param(model.P, initialize=marketingdict, \n", " doc='marketing for item p')\n", "model.capital = pyo.Param(model.P, initialize=capitaldict, \n", " doc='capital for item p')\n", "\n", "# Define decision variables\n", "model.PRODUCTS = pyo.Var(model.P, within=pyo.NonNegativeIntegers, \n", " doc='number of items to put on floor')\n", "\n", "# Define objective\n", "# Note that the default objective sense is minimize in Pyomo\n", "def profit_rule(model):\n", " return sum(model.profit[p] * model.PRODUCTS[p] for p in model.P)\n", "model.cost = pyo.Objective(rule=profit_rule, sense=pyo.maximize)\n", "\n", "# Define constraints\n", "def floorspace_rule(model):\n", " return (\n", " sum(model.PRODUCTS[p] * model.floorspace[p] for p in model.P) \n", " <= floorspace_avail \n", " )\n", "model.avail_floorspace_rule = pyo.Constraint(rule=floorspace_rule, \n", " doc='Floor space constraint')\n", "\n", "def marketing_rule(model):\n", " return (\n", " sum(model.PRODUCTS[p] * model.marketing[p] for p in model.P)\n", " <= 0 \n", " )\n", "model.avail_marketing_rule = pyo.Constraint(rule=marketing_rule, \n", " doc='Marketing constraint')\n", "\n", "def capital_rule(model):\n", " return (\n", " sum(model.PRODUCTS[p] * model.capital[p] for p in model.P) \n", " <= capital_avail \n", " )\n", "model.avail_capital_rule = pyo.Constraint(rule=capital_rule, \n", " doc='Capital constraint')\n", "\n", "# Set the solver to be GLPK\n", "opt = pyo.SolverFactory(\"glpk\")\n", "\n", "# Call the solver\n", "results = opt.solve(model, tee=False)\n", "\n", "# If it's optimal, print the optimal objective value and solution\n", "if str(results.solver.termination_condition) == \"optimal\":\n", " print(\"Results\\n\")\n", " print(f\"Maximum Total Profit is {pyo.value(model.cost):0.2f} dollars.\\n\")\n", " print(\"Product Number Proft Made \")\n", " print(\"------------- ----------- -----------\")\n", " for p in model.P:\n", " if model.PRODUCTS[p].value != 0:\n", " # Note that multiple strings (plain or f-string) in a print statement \n", " # (with nothing in between) will be concatenated\n", " print(f\"{p:13s} {model.PRODUCTS[p].value:11.2f} \"\n", " f\"{model.PRODUCTS[p].value*model.profit[p]:11.2f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Problem 2.__\n", "A small engineering consulting firm has 3 senior designers available to work on the firm's 4 current projects over the next 2 weeks. Each designer has 80 hours to split among the projects, and the following table shows the manager's scoring (0 = nil to 100 = perfect) of the capability of each designer to contribute to each project, along with his estimate of the hours that each project will require. \n", "\n", "\\begin{array}{l|rrrr}\n", "& \\mbox{Project} \\\\ \\hline\n", " \\mbox{Designer} & \\mbox{1} & \\mbox{2} & \\mbox{3} & \\mbox{4} \\\\ \\hline\n", "\\mbox{1} & 90 & 80 & 10 & 50 \\\\\n", "\\mbox{2} & 60 & 70 & 50 & 65 \\\\\n", "\\mbox{3} & 70 & 40 & 80 & 85 \\\\ \\hline\n", "\\mbox{Required} & 70 & 50 & 85 & 35\n", "\\end{array}\n", "\n", "Below is a linear program using set-based notation that maximizes the number of hours spent by the most qualified designers for a given project, while ensuring that all projects are completed according to the given requirements.\n", "\n", "- Implement the linear program below with Pyomo. \n", "- To populate the data, look at the Excel workbook `engineering_consulting.xlsx`, located in the same folder as this notebook. Use xlwings to read the data from the Excel Workbook on the `Input` sheet. \n", "- Solve the linear program.\n", "- Use `xlwings` to write the results back to the Excel spreadsheet on the `Output` worksheet: \n", " - Write the value of the objective function to the cell `A2`.\n", " - Write the values of the decision variable $H_{d,p}$ in the cells `B7:E9`.\n", "\n", "_Indices and sets._\n", "\\begin{array}{ll}\n", " \\mbox{$d \\in D$} & \\mbox{designer}, D=\\{\\mbox{1, 2, 3} \\} \\\\\n", " \\mbox{$p \\in P$} & \\mbox{project}, P=\\{\\mbox{1, 2, 3, 4} \\} \\\\\n", "\\end{array}\n", "\n", "_Data._\n", "\\begin{array}{ll}\n", " \\mbox{score}_{d,p} & \\mbox{manager's score for designer $d$ working on project $p$} \\\\\n", " \\mbox{required}_p & \\mbox{manager's estimate of hours required to complete project $p$} \\\\\n", " \\mbox{hours}_d & \\mbox{number of hours designer $d$ has to split among projects} \\\\\n", "\\end{array}\n", "\n", "_Decision variables._\n", "\\begin{array}{ll}\n", " \\mbox{$H_{d,p}$} & \\mbox{number of hours worked by designer $d$ on project $p$} \\\\\n", "\\end{array}\n", "\n", "_Formulation._\n", "\\begin{alignat}{3}\n", "\\max \\quad & \\sum_{d \\in D}\\sum_{p\\in P} \\mbox{score}_{d,p} * H_{d,p} \\\\\n", "\\mbox{s.t.} \\quad & \\sum_{p\\in P} H_{d,p} \\le \\mbox{hours}_d && \\forall d \\in D \\\\\n", " & \\sum_{d \\in D} H_{d,p} \\ge req_p && \\forall p \\in P\\\\\n", " & H_{d,p} \\ge 0 && \\forall d \\in D, p\\in P\n", "\\end{alignat}" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Results\n", "\n", "Maximum Objective is 18825.00.\n", "\n", " Project\n", "Designer 1 2 3 4\n", "-------- ----- ----- ----- -----\n", " 1 70.0 10.0 0.0 0.0\n", " 2 0.0 40.0 5.0 35.0\n", " 3 0.0 0.0 80.0 0.0\n" ] } ], "source": [ "# Import xlwings\n", "import xlwings as xw\n", "\n", "# Define the model as a Concrete Model in Pyomo\n", "model = pyo.ConcreteModel()\n", "\n", "# Define sets \n", "model.D = pyo.Set(initialize=[1,2,3], doc='Set of designers')\n", "model.P = pyo.Set(initialize=[1,2,3,4], doc='Set of projects')\n", "\n", "# Define data and parameters\n", "#\n", "# Load data from Excel\n", "#\n", "# Load data workbook\n", "wb = xw.Book(r'engineering_consulting.xlsx')\n", "sinput = wb.sheets['Input']\n", "soutput = wb.sheets['Output']\n", "\n", "# Read ALL the entries from Excel into memory for score, \n", "# then assign them to the scoredict dictionary in a nested loop\n", "scoredict = {}\n", "rng1 = sinput.range('C3').expand('table')\n", "\n", "for d in model.D:\n", " for p in model.P:\n", " scoredict[d, p] = rng1[d - 1, p - 1].value\n", "\n", "model.score = pyo.Param(\n", " model.D, model.P, initialize=scoredict, \n", " doc='score if designer d is working on project p'\n", ")\n", "\n", "# Read ALL the entries from Excel into memory for required,\n", "# then assign them to the requireddict dictionary in a loop\n", "requireddict = {}\n", "rng2 = sinput.range('C7').expand('right')\n", "\n", "for p in model.P:\n", " requireddict[p] = rng2[p - 1].value\n", " \n", "model.required = pyo.Param(\n", " model.P, initialize=requireddict, \n", " doc='manager estimate for time required to complete project p'\n", ")\n", "\n", "# Below we assign a value of 80 for EACH designer.\n", "# This could change based on input from the client, \n", "# so having it be \"easy\" to change as a parameter \n", "# is potentially helpful for later.\n", "model.hours = pyo.Param(\n", " model.D, initialize=80,\n", " doc='number of hours designer d has to split among projects'\n", ")\n", "\n", "# Define decision variables\n", "model.H = pyo.Var(\n", " model.D, model.P, within=pyo.NonNegativeReals, \n", " doc='number of hours worked by designer d on project p'\n", ")\n", "\n", "# Define objective\n", "# Recall that default sense is minimize in Pyomo\n", "def qualified_rule(model):\n", " return sum(model.score[d, p] * model.H[d, p] \n", " for d in model.D for p in model.P)\n", "model.cost = pyo.Objective(rule=qualified_rule, sense=pyo.maximize)\n", "\n", "# Define constraints\n", "def hours_rule(model, d):\n", " return sum(model.H[d, p] for p in model.P) <= model.hours[d] \n", "model.avail_hours_rule = pyo.Constraint(model.D, rule=hours_rule, \n", " doc='Designer hours bound')\n", "\n", "def project_rule(model, p):\n", " return sum(model.H[d, p] for d in model.D) >= model.required[p] \n", "model.req_project_rule = pyo.Constraint(model.P, rule=project_rule, \n", " doc='Project hours bound')\n", "\n", "# Set the solver to be GLPK\n", "opt = pyo.SolverFactory(\"glpk\")\n", "\n", "# Call the solver\n", "results = opt.solve(model, tee=False)\n", "\n", "# If it's optimal, print the optimal objective value and solution\n", "if str(results.solver.termination_condition)==\"optimal\":\n", " print(\"Results\\n\")\n", " print(f\"Maximum Objective is {pyo.value(model.cost):0.2f}.\\n\")\n", " print(\" Project\")\n", " print(\"Designer 1 2 3 4\")\n", " print(\"-------- ----- ----- ----- -----\")\n", " for d in model.D:\n", " print(f\"{d:8d} {model.H[d, 1].value:5.1f} \"\n", " f\"{model.H[d, 2].value:5.1f} {model.H[d, 3].value:5.1f} \"\n", " f\"{model.H[d, 4].value:5.1f}\")\n", " \n", " # Now use xlwings to write output back to spreadsheet\n", " # Clear current output in worksheet\n", " wb.sheets['Output'].range('A2').clear()\n", " wb.sheets['Output'].range('B7:E9').clear()\n", " \n", " # Output objective value to Excel\n", " soutput.range('A2').value = pyo.value(model.cost)\n", " \n", " # Output solution to Excel\n", " counter = 0\n", " for d in model.D:\n", " soutput.range(f'B{counter + 7}').value = model.H[d,1].value\n", " soutput.range(f'C{counter + 7}').value = model.H[d,2].value\n", " soutput.range(f'D{counter + 7}').value = model.H[d,3].value\n", " soutput.range(f'E{counter + 7}').value = model.H[d,4].value\n", " counter = counter + 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Bonus classwork" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The problem below is _supplemental material_. You should only work on it after you have completed the two probelms above." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Problem 3. (Winston, Operations Research: Applications and Algorithms, Fourth Edition, Example 7 on page 72)__\n", "\n", "A post office requires different numbers of full-time employees on different days of the week. The number of full-time employees required on each day is given in the Excel workbook `post_office.xlsx`, located in the same folder as this notebook. Union rules state that each full-time employee must work five consecutive days and then receive two days off. For example, an employee who works Monday to Friday must be off on Saturday and Sunday. The post office wants to meet its daily requirements using only full-time employees. \n", "\n", " - Formulate and solve an LP that the post office can use to minimize the number of full-time employees who must be hired. \n", " - Use the data from `Sheet1`. \n", " - When you have solved the model, write your solution and optimal objective function value back into the spreadsheet in the cells provided.\n", " - *Hint.* let $X_i$ be the number of employees __starting their work week__ on day $i$." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "X : number of workers starting on day d\n", " Size=7, Index=D\n", " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", " 0 : 0 : 2.0 : None : False : False : NonNegativeReals\n", " 1 : 0 : 7.0 : None : False : False : NonNegativeReals\n", " 2 : 0 : 0.0 : None : False : False : NonNegativeReals\n", " 3 : 0 : 9.0 : None : False : False : NonNegativeReals\n", " 4 : 0 : 6.0 : None : False : False : NonNegativeReals\n", " 5 : 0 : 0.0 : None : False : False : NonNegativeReals\n", " 6 : 0 : 0.0 : None : False : False : NonNegativeReals\n" ] } ], "source": [ "# Let X[0] = number of workers starting on Monday, ..., \n", "# X[6] = number of workers starting on Sunday. \n", "\n", "# Then want to minimize X[0] + X[1] + X[2] + X[3] + X[4] + X[5] + X[6]\n", "# subject to: \n", "# X[0] + X[3] + X[4] + X[5] + X[6] >= MON workers required\n", "# X[0] + X[1] + X[4] + X[5] + X[6] >= TUE workers required\n", "# X[0] + X[1] + X[2] + X[5] + X[6] >= WED workers required\n", "# X[0] + X[1] + X[2] + X[3] + X[6] >= THU workers required\n", "# X[0] + X[1] + X[2] + X[3] + X[4] >= FRI workers required\n", "# X[1] + X[2] + X[3] + X[4] + X[5] >= SAT workers required\n", "# X[2] + X[3] + X[4] + X[5] + X[6] >= SUN workers required\n", "# X[0], ..., X[6] >= 0\n", "\n", "# Define the model as a Concrete Model in Pyomo\n", "model = pyo.ConcreteModel()\n", "\n", "# Define sets \n", "model.D = pyo.Set(initialize=[0, 1, 2, 3, 4, 5, 6], \n", " doc='Set of days of the week -> 0 is Monday')\n", "\n", "# Define data and parameters\n", "# Read data from Excel\n", "wb = xw.Book('post_office.xlsx')\n", "sht = wb.sheets('Sheet1')\n", "reqs = sht.range('B3').expand('down').value\n", "\n", "# Read ALL the entries from Excel for number of \n", "# required workers on each day of the week and \n", "# then assign them to the reqsdict dictionary\n", "reqsdict = {}\n", "rng1 = sht.range('B3').expand('down')\n", "\n", "for d in model.D:\n", " reqsdict[d] = rng1[d].value\n", "\n", "model.reqs = pyo.Param(model.D, initialize=reqsdict, \n", " doc='number of full time employees required on day d')\n", "\n", "# Define decision variables\n", "model.X = pyo.Var(model.D, within=pyo.NonNegativeReals, \n", " doc='number of workers starting on day d')\n", "\n", "# Define objective\n", "# Recall that default sense is minimize in Pyomo\n", "def total_workers_rule(model):\n", " return sum(model.X[d] for d in model.D)\n", "model.cost = pyo.Objective(rule=total_workers_rule, sense=pyo.minimize)\n", "\n", "# Define constraints\n", "def monday_rule(model):\n", " return model.X[0] + model.X[3] + model.X[4] + model.X[5] + model.X[6] \\\n", " >= model.reqs[0] \n", "model.avail_monday_rule = pyo.Constraint(rule=monday_rule, \n", " doc='Obey Monday workers constraint')\n", "\n", "def tuesday_rule(model):\n", " return model.X[0] + model.X[1] + model.X[4] + model.X[5] + model.X[6] \\\n", " >= model.reqs[1] \n", "model.avail_tuesday_rule = pyo.Constraint(rule=tuesday_rule, \n", " doc='Obey Tuesday workers constraint')\n", "\n", "def wednesday_rule(model):\n", " return model.X[0] + model.X[1] + model.X[2] + model.X[5] + model.X[6] \\\n", " >= model.reqs[2] \n", "model.avail_wednesday_rule = pyo.Constraint(rule=monday_rule, \n", " doc='Obey Wednesday workers constraint')\n", "\n", "def thursday_rule(model):\n", " return model.X[0] + model.X[1] + model.X[2] + model.X[3] + model.X[6] \\\n", " >= model.reqs[3] \n", "model.avail_thursday_rule = pyo.Constraint(rule=thursday_rule, \n", " doc='Obey Thursday workers constraint')\n", "\n", "def friday_rule(model):\n", " return model.X[0] + model.X[1] + model.X[2] + model.X[3] + model.X[4] \\\n", " >= model.reqs[4] \n", "model.avail_friday_rule = pyo.Constraint(rule=friday_rule, \n", " doc='Obey Friday workers constraint')\n", "\n", "def saturday_rule(model):\n", " return model.X[1] + model.X[2] + model.X[3] + model.X[4] + model.X[5] \\\n", " >= model.reqs[5] \n", "model.avail_saturday_rule = pyo.Constraint(rule=saturday_rule, \n", " doc='Obey Saturday workers constraint')\n", "\n", "def sunday_rule(model):\n", " return model.X[2] + model.X[3] + model.X[4] + model.X[5] + model.X[6] \\\n", " >= model.reqs[6] \n", "model.avail_sunday_rule = pyo.Constraint(rule=sunday_rule, \n", " doc='Obey Sunday workers constraint')\n", "\n", "# Set the solver to be GLPK\n", "opt = pyo.SolverFactory(\"glpk\")\n", "\n", "# Call the solver\n", "results = opt.solve(model, tee=False)\n", "\n", "# If it's optimal, print the optimal objective value and solution\n", "if str(results.solver.termination_condition) == \"optimal\":\n", "\n", " # display output\n", " model.X.display()\n", "\n", " # Now use xlwings to write output back to spreadsheet\n", " # Clear current output in worksheet\n", " sht.range('C3:C9').clear()\n", " \n", " # Output objective value to Excel\n", " sht.range('B11').value = pyo.value(model.cost)\n", " \n", " # Output solution to Excel\n", " counter = 0\n", " for d in model.D:\n", " sht.range(f'C{counter + 3}').value = model.X[d].value\n", " counter = counter + 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now solve the same problem but use the data on `Sheet2`. What difficulties do you anticipate in using the optimal solution to staff the post office? How can you modify the Pyomo model to provide a solution that the post office can actually use for this second set of data? " ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "X : number of workers starting on day d\n", " Size=7, Index=D\n", " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", " 0 : 0 : 5.66666666666667 : None : False : False : NonNegativeReals\n", " 1 : 0 : 4.66666666666667 : None : False : False : NonNegativeReals\n", " 2 : 0 : 0.0 : None : False : False : NonNegativeReals\n", " 3 : 0 : 8.66666666666667 : None : False : False : NonNegativeReals\n", " 4 : 0 : 2.66666666666667 : None : False : False : NonNegativeReals\n", " 5 : 0 : 0.0 : None : False : False : NonNegativeReals\n", " 6 : 0 : 0.0 : None : False : False : NonNegativeReals\n" ] } ], "source": [ "# Let X[0] = number of workers starting on Monday, ..., \n", "# X[6] = number of workers starting on Sunday. \n", "\n", "# Then want to minimize X[0] + X[1] + X[2] + X[3] + X[4] + X[5] + X[6]\n", "# subject to: \n", "# X[0] + X[3] + X[4] + X[5] + X[6] >= MON workers required\n", "# X[0] + X[1] + X[4] + X[5] + X[6] >= TUE workers required\n", "# X[0] + X[1] + X[2] + X[5] + X[6] >= WED workers required\n", "# X[0] + X[1] + X[2] + X[3] + X[6] >= THU workers required\n", "# X[0] + X[1] + X[2] + X[3] + X[4] >= FRI workers required\n", "# X[1] + X[2] + X[3] + X[4] + X[5] >= SAT workers required\n", "# X[2] + X[3] + X[4] + X[5] + X[6] >= SUN workers required\n", "# X[0], ..., X[6] >= 0\n", "\n", "# Define the model as a Concrete Model in Pyomo\n", "model = pyo.ConcreteModel()\n", "\n", "# Define sets \n", "model.D = pyo.Set(initialize=[0, 1, 2, 3, 4, 5, 6], \n", " doc='Set of days of the week -> 0 is Monday')\n", "\n", "# Define data and parameters\n", "# Read data from Excel\n", "wb = xw.Book('post_office.xlsx')\n", "sht = wb.sheets('Sheet2')\n", "reqs = sht.range('B3').expand('down').value\n", "\n", "# Read ALL the entries from Excel for number of \n", "# required workers on each day of the week and \n", "# then assign them to the reqsdict dictionary\n", "reqsdict = {}\n", "rng1 = sht.range('B3').expand('down')\n", "\n", "for d in model.D:\n", " reqsdict[d] = rng1[d].value\n", "\n", "model.reqs = pyo.Param(model.D, initialize=reqsdict, \n", " doc='number of full time employees required on day d')\n", "\n", "# Define decision variables\n", "model.X = pyo.Var(model.D, within=pyo.NonNegativeReals, \n", " doc='number of workers starting on day d')\n", "\n", "# Define objective\n", "# Recall that default sense is minimize in Pyomo\n", "def total_workers_rule(model):\n", " return sum(model.X[d] for d in model.D)\n", "model.cost = pyo.Objective(rule=total_workers_rule, sense=pyo.minimize)\n", "\n", "# Define constraints\n", "def monday_rule(model):\n", " return model.X[0] + model.X[3] + model.X[4] + model.X[5] + model.X[6] \\\n", " >= model.reqs[0] \n", "model.avail_monday_rule = pyo.Constraint(rule=monday_rule, \n", " doc='Obey Monday workers constraint')\n", "\n", "def tuesday_rule(model):\n", " return model.X[0] + model.X[1] + model.X[4] + model.X[5] + model.X[6] \\\n", " >= model.reqs[1] \n", "model.avail_tuesday_rule = pyo.Constraint(rule=tuesday_rule, \n", " doc='Obey Tuesday workers constraint')\n", "\n", "def wednesday_rule(model):\n", " return model.X[0] + model.X[1] + model.X[2] + model.X[5] + model.X[6] \\\n", " >= model.reqs[2] \n", "model.avail_wednesday_rule = pyo.Constraint(rule=monday_rule, \n", " doc='Obey Wednesday workers constraint')\n", "\n", "def thursday_rule(model):\n", " return model.X[0] + model.X[1] + model.X[2] + model.X[3] + model.X[6] \\\n", " >= model.reqs[3] \n", "model.avail_thursday_rule = pyo.Constraint(rule=thursday_rule, \n", " doc='Obey Thursday workers constraint')\n", "\n", "def friday_rule(model):\n", " return model.X[0] + model.X[1] + model.X[2] + model.X[3] + model.X[4] \\\n", " >= model.reqs[4] \n", "model.avail_friday_rule = pyo.Constraint(rule=friday_rule, \n", " doc='Obey Friday workers constraint')\n", "\n", "def saturday_rule(model):\n", " return model.X[1] + model.X[2] + model.X[3] + model.X[4] + model.X[5] \\\n", " >= model.reqs[5] \n", "model.avail_saturday_rule = pyo.Constraint(rule=saturday_rule, \n", " doc='Obey Saturday workers constraint')\n", "\n", "def sunday_rule(model):\n", " return model.X[2] + model.X[3] + model.X[4] + model.X[5] + model.X[6] \\\n", " >= model.reqs[6] \n", "model.avail_sunday_rule = pyo.Constraint(rule=sunday_rule, \n", " doc='Obey Sunday workers constraint')\n", "\n", "# Set the solver to be GLPK\n", "opt = pyo.SolverFactory(\"glpk\")\n", "\n", "# Call the solver\n", "results = opt.solve(model, tee=False)\n", "\n", "# If it's optimal, print the optimal objective value and solution\n", "if str(results.solver.termination_condition) == \"optimal\":\n", "\n", " # display output\n", " model.X.display()\n", "\n", " # Now use xlwings to write output back to spreadsheet\n", " # Clear current output in worksheet\n", " sht.range('C3:C9').clear()\n", " \n", " # Output objective value to Excel\n", " sht.range('B11').value = pyo.value(model.cost)\n", " \n", " # Output solution to Excel\n", " counter = 0\n", " for d in model.D:\n", " sht.range(f'C{counter + 3}').value = model.X[d].value\n", " counter = counter + 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now fix the domain of the decision variables to be `NonNegativeIntegers` so we can get a solution the post office can actually use." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "X : number of workers starting on day d\n", " Size=7, Index=D\n", " Key : Lower : Value : Upper : Fixed : Stale : Domain\n", " 0 : 0 : 6.0 : None : False : False : NonNegativeIntegers\n", " 1 : 0 : 4.0 : None : False : False : NonNegativeIntegers\n", " 2 : 0 : 0.0 : None : False : False : NonNegativeIntegers\n", " 3 : 0 : 9.0 : None : False : False : NonNegativeIntegers\n", " 4 : 0 : 3.0 : None : False : False : NonNegativeIntegers\n", " 5 : 0 : 0.0 : None : False : False : NonNegativeIntegers\n", " 6 : 0 : 0.0 : None : False : False : NonNegativeIntegers\n" ] } ], "source": [ "# Let X[0] = number of workers starting on Monday, ..., \n", "# X[6] = number of workers starting on Sunday. \n", "\n", "# Then want to minimize X[0] + X[1] + X[2] + X[3] + X[4] + X[5] + X[6]\n", "# subject to: \n", "# X[0] + X[3] + X[4] + X[5] + X[6] >= MON workers required\n", "# X[0] + X[1] + X[4] + X[5] + X[6] >= TUE workers required\n", "# X[0] + X[1] + X[2] + X[5] + X[6] >= WED workers required\n", "# X[0] + X[1] + X[2] + X[3] + X[6] >= THU workers required\n", "# X[0] + X[1] + X[2] + X[3] + X[4] >= FRI workers required\n", "# X[1] + X[2] + X[3] + X[4] + X[5] >= SAT workers required\n", "# X[2] + X[3] + X[4] + X[5] + X[6] >= SUN workers required\n", "# X[0], ..., X[6] >= 0\n", "\n", "# Define the model as a Concrete Model in Pyomo\n", "model = pyo.ConcreteModel()\n", "\n", "# Define sets \n", "model.D = pyo.Set(initialize=[0, 1, 2, 3, 4, 5, 6], \n", " doc='Set of days of the week -> 0 is Monday')\n", "\n", "# Define data and parameters\n", "# Read data from Excel\n", "wb = xw.Book('post_office.xlsx')\n", "sht = wb.sheets('Sheet2')\n", "reqs = sht.range('B3').expand('down').value\n", "\n", "# Read ALL the entries from Excel for number of \n", "# required workers on each day of the week and \n", "# then assign them to the reqsdict dictionary\n", "reqsdict = {}\n", "rng1 = sht.range('B3').expand('down')\n", "\n", "for d in model.D:\n", " reqsdict[d] = rng1[d].value\n", "\n", "model.reqs = pyo.Param(model.D, initialize=reqsdict, \n", " doc='number of full time employees required on day d')\n", "\n", "# Define decision variables\n", "model.X = pyo.Var(model.D, within=pyo.NonNegativeIntegers, \n", " doc='number of workers starting on day d')\n", "\n", "# Define objective\n", "# Recall that default sense is minimize in Pyomo\n", "def total_workers_rule(model):\n", " return sum(model.X[d] for d in model.D)\n", "model.cost = pyo.Objective(rule=total_workers_rule, sense=pyo.minimize)\n", "\n", "# Define constraints\n", "def monday_rule(model):\n", " return model.X[0] + model.X[3] + model.X[4] + model.X[5] + model.X[6] \\\n", " >= model.reqs[0] \n", "model.avail_monday_rule = pyo.Constraint(rule=monday_rule, \n", " doc='Obey Monday workers constraint')\n", "\n", "def tuesday_rule(model):\n", " return model.X[0] + model.X[1] + model.X[4] + model.X[5] + model.X[6] \\\n", " >= model.reqs[1] \n", "model.avail_tuesday_rule = pyo.Constraint(rule=tuesday_rule, \n", " doc='Obey Tuesday workers constraint')\n", "\n", "def wednesday_rule(model):\n", " return model.X[0] + model.X[1] + model.X[2] + model.X[5] + model.X[6] \\\n", " >= model.reqs[2] \n", "model.avail_wednesday_rule = pyo.Constraint(rule=monday_rule, \n", " doc='Obey Wednesday workers constraint')\n", "\n", "def thursday_rule(model):\n", " return model.X[0] + model.X[1] + model.X[2] + model.X[3] + model.X[6] \\\n", " >= model.reqs[3] \n", "model.avail_thursday_rule = pyo.Constraint(rule=thursday_rule, \n", " doc='Obey Thursday workers constraint')\n", "\n", "def friday_rule(model):\n", " return model.X[0] + model.X[1] + model.X[2] + model.X[3] + model.X[4] \\\n", " >= model.reqs[4] \n", "model.avail_friday_rule = pyo.Constraint(rule=friday_rule, \n", " doc='Obey Friday workers constraint')\n", "\n", "def saturday_rule(model):\n", " return model.X[1] + model.X[2] + model.X[3] + model.X[4] + model.X[5] \\\n", " >= model.reqs[5] \n", "model.avail_saturday_rule = pyo.Constraint(rule=saturday_rule, \n", " doc='Obey Saturday workers constraint')\n", "\n", "def sunday_rule(model):\n", " return model.X[2] + model.X[3] + model.X[4] + model.X[5] + model.X[6] \\\n", " >= model.reqs[6] \n", "model.avail_sunday_rule = pyo.Constraint(rule=sunday_rule, \n", " doc='Obey Sunday workers constraint')\n", "\n", "# Set the solver to be GLPK\n", "opt = pyo.SolverFactory(\"glpk\")\n", "\n", "# Call the solver\n", "results = opt.solve(model, tee=False)\n", "\n", "# If it's optimal, print the optimal objective value and solution\n", "if str(results.solver.termination_condition) == \"optimal\":\n", "\n", " # display output\n", " model.X.display()\n", "\n", " # Now use xlwings to write output back to spreadsheet\n", " # Clear current output in worksheet\n", " sht.range('C3:C9').clear()\n", " \n", " # Output objective value to Excel\n", " sht.range('B11').value = pyo.value(model.cost)\n", " \n", " # Output solution to Excel\n", " counter = 0\n", " for d in model.D:\n", " sht.range(f'C{counter + 3}').value = model.X[d].value\n", " counter = counter + 1" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.6" } }, "nbformat": 4, "nbformat_minor": 2 }