pooling/std_pooling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Standard Pooling Problem\n",
    "\n",
    "## Objective and Prerequisites\n",
    "\n",
    "Companies across numerous industries – including petrochemical refining, wastewater treatment, and mining – use mathematical optimization to solve the pooling problem. In this example, we’ll guide you through the process of building a mixed-integer quadratically-constrained programming (MIQCP) model of a pooling problem using the Gurobi Python API and show you how to generate an optimal solution to the problem with the Gurobi Optimizer.\n",
    "\n",
    "This modeling example is at the advanced level, where we assume that you know Python and the Gurobi Python API and that you have advanced knowledge of building mathematical optimization models. Typically, the objective function and/or constraints of these examples are complex or require advanced features of the Gurobi Python API.\n",
    "\n",
    "**Download the Repository** <br /> \n",
    "You can download the repository containing this and other examples by clicking [here](https://github.com/Gurobi/modeling-examples/archive/master.zip). \n",
    "\n",
    "---\n",
    "## Motivation\n",
    "\n",
    "The pooling problem is a challenging problem in the petrochemical refining, wastewater treatment and mining industries. This problem can be regarded as a generalization of the minimum-cost flow problem and the blending problem. It is indeed important because of the significant savings it can generate, so it comes at no surprise that it has been studied extensively since Haverly pointed out the non-linear structure of this problem in 1978 [5].\n",
    "\n",
    "---\n",
    "## Problem Description\n",
    "\n",
    "The Minimum-Cost Flow Problem (MCFP) seeks to find the cheapest way of sending a certain amount of flow from a set of source nodes to a set of target nodes, possibily via transshipment nodes⁠, in a directed capacitated network. The Blending Problem is a type of MCFP with only source and target nodes, where raw materials with different attribute qualities are blended together to create end products in such a way that their attribute qualities are within tolerances.\n",
    "\n",
    "The Pooling Problem combines features of both problems, as flow streams from different sources are mixed at intermediate pools and blended again at the target nodes. The non-linearity is in fact the direct result of considering pools, as the quality of a given attribute at a pool —defined as the weighted average of the qualities of the incoming streams— is an unknown quantity and thus needs to be captured by a decision variable. We refer to this problem as the Standard Pooling Problem when the network can be represented by a tripartite graph, i.e. three disjoint sets of nodes such that no nodes within the same set are adjacent. In a nutshell, it can be stated as follows: Given a list of source nodes with raw materials containing known attribute qualities, what is the cheapest way of mixing these materials at intermediate pools so as to meet the demand and tolerances at multiple target nodes? (Gupte et al., 2017) [4]. Several different formulations for the Standard Pooling Problem and its extensions exist in the literature, which can be classified into two main categories: one that consists of flow and quality variables, and the other that uses flow proportions instead of quality variables. Both categories will be considered in this notebook.\n",
    "\n",
    "---\n",
    "## Problem Instance\n",
    "\n",
    "As an illustrative example, we will solve the second Pooling Problem posed by Rehfeldt and Tisljar in 1997 and cited by Audet et al. in 2004:\n",
    "\n",
    "![Graph](rehfeldt_tisljar2_graph.png)\n",
    "\n",
    "![Table](rehfeldt_tisljar2_table.png)\n",
    "\n",
    "To that end, let's declare the required data structures to represent this problem instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install gurobipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from itertools import product\n",
    "\n",
    "import gurobipy as gp\n",
    "from gurobipy import GRB\n",
    "\n",
    "# tested with Gurobi v9.1.0 and Python 3.7.0\n",
    "\n",
    "attrs = {'den', 'bnz', 'roz', 'moz'}\n",
    "\n",
    "sources, cost, supply, content = gp.multidict({\n",
    "    \"s1\": [49.2, 6097.56, {'den': 0.82, 'bnz':3, 'roz':99.2,'moz':90.5}],\n",
    "    \"s2\": [62.0, 16129, {'den': 0.62, 'bnz':0, 'roz':87.9,'moz':83.5}],\n",
    "    \"s3\": [300.0, 500, {'den': 0.75, 'bnz':0, 'roz':114,'moz':98.7}]\n",
    "})\n",
    "\n",
    "targets, price, demand, min_tol, max_tol = gp.multidict({\n",
    "    \"t1\": [190, 500, {'den': 0.74, 'roz':95,'moz':85}, {'den': 0.79}],\n",
    "    \"t2\": [230, 500, {'den': 0.74, 'roz':96,'moz':88}, {'den': 0.79, 'bnz':0.9}],\n",
    "    \"t3\": [150, 500, {'den': 0.74, 'roz':91}, {'den': 0.79}]\n",
    "})\n",
    "\n",
    "pools, cap = gp.multidict({\n",
    "    \"p1\": 1250,\n",
    "    \"p2\": 1750\n",
    "})\n",
    "\n",
    "# The function `product` deploys the Cartesian product of elements in sets A and B\n",
    "s2p = set(product(sources, pools))\n",
    "p2t = set(product(pools, targets))\n",
    "s2t = {(\"s1\", \"t2\"),\n",
    "       (\"s2\", \"t1\"),\n",
    "       (\"s2\", \"t3\"),\n",
    "       (\"s3\", \"t1\")}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Solution Approach\n",
    "\n",
    "Mathematical programming is a declarative approach where the modeler formulates a mathematical optimization model that captures the key aspects of a complex decision problem. The Gurobi Optimizer solves such models using state-of-the-art mathematics and computer science.\n",
    "\n",
    "A mathematical optimization model has five components, namely:\n",
    "\n",
    "- Sets and indices.\n",
    "- Parameters.\n",
    "- Decision variables.\n",
    "- Objective function(s).\n",
    "- Constraints.\n",
    "\n",
    "A quadratic constraint that involves only products of disjoint pairs of variables is often called a bilinear constraint, and a model that contains bilinear constraints is often called a Bilinear Program. Bilinear constraints are a special case of non-convex quadratic constraints. This type of problems is typically solved using spatial Branch and Bound (sB&B). This algorithm explores the entire search space, so it provides a globally valid lower bound on the optimal objective value and —given enough time— it will find a globally optimal solution (subject to tolerances). The interested reader is referred to [references](#references) [3], [6] and [7].\n",
    "\n",
    "We now present two alternative Bilinear Programs for the Standard Pooling Problem:\n",
    "\n",
    "### P-formulation (Concentration)\n",
    "\n",
    "#### Sets and Indices\n",
    "\n",
    "$G=(V,E)$: Directed graph.\n",
    "\n",
    "$i,j \\in V$: Set of nodes.\n",
    "\n",
    "$(i,j) \\in E \\subset  V \\times V$: Set of edges.\n",
    "\n",
    "$N(i)^+ = \\{j \\in V \\mid (i,j) \\in E \\}$: Set of successor nodes receiving outflow from node $i$.\n",
    "\n",
    "$N(j)^- = \\{i \\in V \\mid (i,j) \\in E \\}$: Set of predecessor nodes sending inflow to node $i$.\n",
    "\n",
    "$k \\in \\text{Attrs}$: Set of attributes.\n",
    "\n",
    "$s \\in \\text{Sources} \\subset V$: Set of source nodes, i.e. $N(s)^-= \\emptyset$.\n",
    "\n",
    "$t \\in \\text{Targets} \\subset V$: Set of target nodes, i.e. $N(t)^+= \\emptyset$.\n",
    "\n",
    "$p \\in \\text{Pools} = V \\setminus (\\text{Sources} \\cup \\text{Targets})$: Set of pools.\n",
    "\n",
    "#### Parameters\n",
    "\n",
    "$\\text{Cost}_s \\in \\mathbb{R}^+$: Cost of acquiring one unit of raw material at source node $s$.\n",
    "\n",
    "$\\text{Supply}_s \\in \\mathbb{R}^+$: Maximum number of units of raw material available at source node $s$.\n",
    "\n",
    "$\\text{Content}_{s,k} \\in \\mathbb{R}^+$: Content of attribute $k$ in raw material at source node $s$.\n",
    "\n",
    "$\\text{Price}_t \\in \\mathbb{R}^+$: Price for selling one unit of final blend at target node $t$.\n",
    "\n",
    "$\\text{Demand}_t \\in \\mathbb{R}^+$: Minimum number of units required of final blend at target node $t$.\n",
    "\n",
    "$\\text{Min_tol}_{t,k} \\in \\mathbb{R}^+$: Minimum tolerance for attribute $k$ in final blend at target node $t$.\n",
    "\n",
    "$\\text{Max_tol}_{t,k} \\in \\mathbb{R}^+$: Maximum tolerance for attribute $k$ in final blend at target node $t$.\n",
    "\n",
    "$\\text{Cap}_p \\in \\mathbb{R}^+$: Maximum Capacity to store intermediate blend at pool $p$.\n",
    "\n",
    "$\\text{UB}_{i,j}\\in \\mathbb{R}^+$: Maximum flow from node $i$ to node $j$.\n",
    "\n",
    "#### Decision Variables\n",
    "\n",
    "$\\text{flow}_{i,j} \\in [0, \\text{UB}_{i,j}]$: Flow from node $i$ to node $j$.\n",
    "\n",
    "$\\text{quality}_{p,k} \\in \\mathbb{R}^+$: Concentration of attribute $k$ at pool $p$.\n",
    "\n",
    "#### Objective Function\n",
    "\n",
    "- **Profit**: Maximize total profits.\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{Max} \\quad Z = \\sum_{t \\in \\text{Targets}}{\\sum_{i \\in N(t)^-}{\\text{Price}_t \\cdot \\text{flow}_{i,t}}} - \\sum_{s \\in \\text{Sources}}{\\sum_{j \\in N(s)^+}{\\text{Cost}_s \\cdot \\text{flow}_{s,j}}}\n",
    "\\tag{0}\n",
    "\\end{equation}\n",
    "\n",
    "#### Constraints\n",
    "\n",
    "- **Flow conservation**: Total inflow of pool $p$ must be equal to its total outflow (nothing is stored in them).\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{t \\in N(p)^+}{\\text{flow}_{p,t}} - \\sum_{s \\in N(p)^-}{\\text{flow}_{s,p}} = 0 \\quad \\forall p \\in \\text{Pools}\n",
    "\\tag{1}\n",
    "\\end{equation}\n",
    "\n",
    "- **Source capacity**: Total outflow of source $s$ cannot exceed its capacity.\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{j \\in N(s)^+}{\\text{flow}_{s,j}} \\leq \\text{Supply}_s \\quad \\forall s \\in \\text{Sources}\n",
    "\\tag{2}\n",
    "\\end{equation}\n",
    "\n",
    "- **Pool capacity**: Total outflow of pool $p$ cannot exceed its capacity.\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{t \\in N(p)^+}{\\text{flow}_{p,t}} \\leq \\text{Cap}_p \\quad \\forall p \\in \\text{Pools}\n",
    "\\tag{3}\n",
    "\\end{equation}\n",
    "\n",
    "- **Target demand**: Total inflow of target $t$ must at least meet its minimum demand.\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{i \\in N(t)^-}{\\text{flow}_{i,t}} \\geq \\text{Demand}_t \\quad \\forall t \\in \\text{Targets}\n",
    "\\tag{4}\n",
    "\\end{equation}\n",
    "\n",
    "- **Pool concentration**: Concentration of attribute $k$ at pool $p$ is expressed as the weighted average (linear blending) of the concentrations associated to the incoming flows (notice the bilinear terms on the right-hand side).\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{s \\in N(p)^-}{\\text{Content}_{s,k} \\cdot \\text{flow}_{s,p}} = \\text{quality}_{p,k} \\cdot \\sum_{t \\in N(p)^+}{\\text{flow}_{p,t}} \\quad \\forall (p,k) \\in \\text{Pools} \\times \\text{Attrs}\n",
    "\\tag{5}\n",
    "\\end{equation}\n",
    "\n",
    "- **Target tolerances**: Concentration of attribute $k$ at target $t$ is also the result of linear blending, and must be within tolerances (notice the bilinear terms on the second expression of the left-hand side).\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{s \\in N(t)^- \\cap \\text{Sources}}{\\text{Content}_{s,k} \\cdot \\text{flow}_{s,t}}+ \\sum_{p \\in N(t)^- \\cap \\text{Pools}}{\\text{quality}_{p,k} \\cdot \\text{flow}_{p,t}} \\geq \\text{Min_tol}_{t,k} \\cdot \\sum_{i \\in N(t)^-}{\\text{flow}_{i,t}} \\quad \\forall (t,k) \\in \\text{Targets} \\times \\text{Attrs}\n",
    "\\tag{6.1}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{s \\in N(t)^- \\cap \\text{Sources}}{\\text{Content}_{s,k} \\cdot \\text{flow}_{s,t}}+ \\sum_{p \\in N(t)^- \\cap \\text{Pools}}{\\text{quality}_{p,k} \\cdot \\text{flow}_{p,t}} \\leq \\text{Max_tol}_{t,k} \\cdot \\sum_{i \\in N(t)^-}{\\text{flow}_{i,t}} \\quad \\forall (t,k) \\in \\text{Targets} \\times \\text{Attrs}\n",
    "\\tag{6.2}\n",
    "\\end{equation}\n",
    "\n",
    "The number of bilinear terms in this formulation is proportional to the number of attributes. An alternative formulation relies on decision variables that represent fractions of flow instead of concentrations, so that the bilinear terms are no longer associated to the number of attributes. Two types of decision variables may be used:\n",
    "\n",
    "- fraction of total inflow at pool $p$, coming from source $s$.\n",
    "- fraction of total outflow at pool $p$, going to terminal $t$.\n",
    "\n",
    "A model based on the first option is now presented:\n",
    "\n",
    "### Q-formulation (Proportion)\n",
    "\n",
    "#### Decision Variables\n",
    "\n",
    "$\\text{flow}_{i,j} \\in [0, \\text{UB}_{i,j}]$: Flow from node $i$ to node $j$.\n",
    "\n",
    "$\\text{prop}_{s,p} \\in \\mathbb{R}^+$: fraction of total inflow at pool $p$, coming from source $s$.\n",
    "\n",
    "**Note:** The $\\text{flow}$ variables from sources to pools are replaced by the $\\text{prop}$ variables.\n",
    "\n",
    "#### Objective Function\n",
    "\n",
    "- **Profit**: Maximize total profits (notice the bilinear terms on the second expression).\n",
    "\n",
    "\\begin{equation}\n",
    "\\text{Max} \\quad Z = \\sum_{t \\in \\text{Targets}}{\\sum_{i \\in N(t)^-}{\\text{Price}_t \\cdot \\text{flow}_{i,t}}} - \\sum_{s \\in \\text{Sources}}{\\text{cost}_s \\cdot \\left( \\sum_{t \\in N(s)^+ \\cap \\text{Targets}}{\\text{flow}_{s,t}} + \\sum_{p \\in N(s)^+ \\cap \\text{Pools}}{\\text{prop}_{s,p} \\cdot \\sum_{t \\in N(p)^+}{\\text{flow}_{p,t}}} \\right) }\n",
    "\\tag{0}\n",
    "\\end{equation}\n",
    "\n",
    "#### Constraints\n",
    "\n",
    "- **Source capacity**: Total outflow of source $s$ cannot exceed its capacity (notice the bilinear terms on the first expression of the left-hand side).\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{p \\in N(s)^+ \\cap \\text{Pools}}{\\text{prop}_{s,p} \\cdot \\sum_{t \\in N(p)^+}{\\text{flow}_{p,t}}} + \\sum_{t \\in N(s)^+ \\cap \\text{Targets}}{\\text{flow}_{s,t}} \\leq \\text{Supply}_s \\quad \\forall s \\in \\text{Sources}\n",
    "\\tag{1}\n",
    "\\end{equation}\n",
    "\n",
    "- **Pool capacity**: Total outflow of pool $p$ cannot exceed its capacity.\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{t \\in N(p)^+}{\\text{flow}_{p,t}} \\leq \\text{Cap}_p \\quad \\forall p \\in \\text{Pools}\n",
    "\\tag{2}\n",
    "\\end{equation}\n",
    "\n",
    "- **Target demand**: Total inflow of target $t$ must at least meet its minimum demand.\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{i \\in N(t)^-}{\\text{flow}_{i,t}} \\geq \\text{Demand}_t \\quad \\forall t \\in \\text{Targets}\n",
    "\\tag{3}\n",
    "\\end{equation}\n",
    "\n",
    "- **Pool inflow**: The sum of the contributions of all incoming sources to pool $p$ must equal one.\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{s \\in N(p)^-}{\\text{prop}_{s,p}} = 1 \\quad \\forall p \\in \\text{Pools}\n",
    "\\tag{4}\n",
    "\\end{equation}\n",
    "\n",
    "- **Target tolerances**: Concentration of attribute $k$ at target $t$ is also the result of linear blending, and must be within tolerances (notice the bilinear terms on the second expression of the left-hand side).\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{s \\in N(t)^- \\cap \\text{Sources}}{\\text{Content}_{s,k} \\cdot \\text{flow}_{s,t}} + \\sum_{p \\in N(t)^- \\cap \\text{Pools}}{\\text{flow}_{p,t} \\cdot \\sum_{s \\in N(p)^-}{\\text{content}_{s,k} \\cdot \\text{prop}_{s,p}}}\n",
    "\\geq \\text{Min_tol}_{t,k} \\cdot \\sum_{i \\in N(t)^-}{\\text{flow}_{i,t}} \\\\\n",
    "\\forall (t,k) \\in \\text{Targets} \\times \\text{Attrs}\n",
    "\\tag{5.1}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\sum_{s \\in N(t)^- \\cap \\text{Sources}}{\\text{Content}_{s,k} \\cdot \\text{flow}_{s,t}} + \\sum_{p \\in N(t)^- \\cap \\text{Pools}}{\\text{flow}_{p,t} \\cdot \\sum_{s \\in N(p)^-}{\\text{content}_{s,k} \\cdot \\text{prop}_{s,p}}}\n",
    "\\leq \\text{Max_tol}_{t,k} \\cdot \\sum_{i \\in N(t)^-}{\\text{flow}_{i,t}} \\\\\n",
    "\\forall (t,k) \\in \\text{Targets} \\times \\text{Attrs}\n",
    "\\tag{5.2}\n",
    "\\end{equation}\n",
    "\n",
    "One drawback is that if some of the source-to-pool edges have flow capacity, we need to define additional constraints instead of just specifying upper bounds on the associated decision variables. Such constraints can be defined with bilinear terms as follows:\n",
    "\n",
    "- **Flow limit**: Flow from source $s$ to pool $p$ cannot exceed the installed capacity (notice the bilinear terms on the left-hand side). In the P-formulation, declaring this is as easy as setting the upper bound of the associated $\\text{flow}$ variable. However, this variable no longer exists in the Q-formulation so we need to model the capacity as a constraint.  \n",
    "\n",
    "\\begin{equation}\n",
    "\\text{prop}_{s,p} \\cdot \\sum_{t \\in N(p)^+}{\\text{flow}_{p,t}} \\leq \\text{UB}_{s,p} \\quad \\forall (i,j) \\in E \\cap \\left( \\text{Sources} \\times \\text{Pools} \\right) \\mid \\text{UB}_{i,j} < \\infty\n",
    "\\tag{6}\n",
    "\\end{equation}\n",
    "\n",
    "---\n",
    "## Python Implementation\n",
    "\n",
    "Solving Bilinear Programs with Gurobi is as easy as configuring the global parameter `nonConvex`. When setting this parameter to a value of 2, non-convex quadratic problems are solved by means of translating them into bilinear form and applying sB&B. We will first deploy the P-formulation model, and then compare it with the Q-formulation model:\n",
    "\n",
    "### P-formulation (Concentration)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using license file c:\\gurobi\\gurobi.lic\n",
      "Set parameter TokenServer to value SANTOS-SURFACE-\n",
      "Changed value of parameter nonConvex to 2\n",
      "   Prev: -1  Min: -1  Max: 2  Default: -1\n",
      "Changed value of parameter timelimit to 300.0\n",
      "   Prev: inf  Min: 0.0  Max: inf  Default: inf\n",
      "Gurobi Optimizer version 9.1.0 build v9.1.0rc0 (win64)\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "Optimize a model with 10 rows, 24 columns and 38 nonzeros\n",
      "Model fingerprint: 0xa0238bae\n",
      "Model has 20 quadratic constraints\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 1e+00]\n",
      "  QMatrix range    [1e+00, 1e+00]\n",
      "  QLMatrix range   [1e-02, 1e+02]\n",
      "  Objective range  [5e+01, 3e+02]\n",
      "  Bounds range     [8e+02, 8e+02]\n",
      "  RHS range        [5e+02, 2e+04]\n",
      "\n",
      "Continuous model is non-convex -- solving as a MIP.\n",
      "\n",
      "Presolve removed 1 rows and 0 columns\n",
      "Presolve time: 0.00s\n",
      "Presolved: 125 rows, 49 columns, 311 nonzeros\n",
      "Presolved model has 24 bilinear constraint(s)\n",
      "Variable types: 49 continuous, 0 integer (0 binary)\n",
      "\n",
      "Root relaxation: objective 9.577646e+05, 56 iterations, 0.00 seconds\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0 957764.599    0   10          - 957764.599      -     -    0s\n",
      "     0     0 957764.599    0   13          - 957764.599      -     -    0s\n",
      "     0     0 878230.977    0   24          - 878230.977      -     -    0s\n",
      "     0     0 878230.977    0   23          - 878230.977      -     -    0s\n",
      "     0     0 878230.977    0   23          - 878230.977      -     -    0s\n",
      "     0     2 878230.977    0   23          - 878230.977      -     -    0s\n",
      "* 7132  6079             193    231816.40594 856384.539   269%   4.0    0s\n",
      "* 8015  6395              88    365634.75700 856384.539   134%   4.4    0s\n",
      "* 8017  6395              89    365637.89650 856384.539   134%   4.4    0s\n",
      "* 8020  6395              88    365784.24521 856384.539   134%   4.4    0s\n",
      "* 9522  6174              64    390877.94434 856384.539   119%   4.4    0s\n",
      "H10416  6354                    411530.70256 856384.539   108%   4.8    0s\n",
      " 73005 47549 infeasible   54      411530.703 836666.047   103%   6.4    5s\n",
      " 155608 107510 785256.517   55   20 411530.703 836482.933   103%   5.5   10s\n",
      " 226878 157501 498463.046  119   20 411530.703 836446.021   103%   5.5   15s\n",
      " 291816 201793 836370.825   70   18 411530.703 836425.538   103%   5.4   20s\n",
      " 362875 249419 430396.666  100   21 411530.703 836413.992   103%   5.5   25s\n",
      " 438817 298604 434094.581   60   16 411530.703 836406.472   103%   5.6   30s\n",
      " 512215 347480 827773.956   59   16 411530.703 836401.082   103%   5.6   35s\n",
      " 585983 396819     cutoff   75      411530.703 836396.959   103%   5.6   40s\n",
      " 662151 446378 653687.436   75   16 411530.703 836393.458   103%   5.7   45s\n",
      " 742681 502708 609442.923   93   16 411530.703 836391.049   103%   5.6   50s\n",
      " 817744 554087 419889.913   89   13 411530.703 836388.846   103%   5.6   55s\n",
      " 882873 596119     cutoff   91      411530.703 836387.401   103%   5.7   60s\n",
      " 959654 651891 471795.319  113   20 411530.703 836385.850   103%   5.7   65s\n",
      " 1036851 705020 836384.089   51   18 411530.703 836384.629   103%   5.6   70s\n",
      " 1104638 751885 423286.104  103   16 411530.703 836383.705   103%   5.6   75s\n",
      " 1175201 801794 836375.090   64   18 411530.703 836382.762   103%   5.6   80s\n",
      " 1250281 854347 836370.765   73   18 411530.703 836382.011   103%   5.6   85s\n",
      " 1319910 901548 836372.202   56   20 411530.703 836381.408   103%   5.6   90s\n",
      " 1394141 955225 836370.340   77   15 411530.703 836380.837   103%   5.6   95s\n",
      " 1465387 1006492 836377.056   55   15 411530.703 836380.302   103%   5.5  100s\n",
      " 1539934 1059497     cutoff   87      411530.703 836379.822   103%   5.5  105s\n",
      " 1617244 1113943 827758.589   62   16 411530.703 836379.402   103%   5.5  110s\n",
      " 1680790 1157475     cutoff  119      411530.703 836379.031   103%   5.5  115s\n",
      " 1747029 1205583 779311.019   67   20 411530.703 836378.685   103%   5.5  120s\n",
      " 1823570 1260783 564225.818  109   17 411530.703 836378.345   103%   5.5  125s\n",
      " 1898494 1310768     cutoff   79      411530.703 836378.038   103%   5.5  130s\n",
      " 1963288 1356267 456248.349  122   18 411530.703 836377.783   103%   5.5  135s\n",
      " 2029351 1403946 836374.293   60   15 411530.703 836377.541   103%   5.5  140s\n",
      " 2106054 1458608 501936.436  121   13 411530.703 836377.278   103%   5.5  145s\n",
      " 2172663 1503450 836370.391   66   24 411530.703 836377.082   103%   5.5  150s\n",
      " 2241601 1552251 608573.498   99   18 411530.703 836376.880   103%   5.5  155s\n",
      " 2317497 1602717 412706.294   91   20 411530.703 836376.658   103%   5.5  160s\n",
      " 2390945 1654509     cutoff   78      411530.703 836376.484   103%   5.5  165s\n",
      " 2467128 1706983 561483.281  134   10 411530.703 836376.297   103%   5.5  170s\n",
      " 2543540 1760478 499198.432  119   15 411530.703 836376.108   103%   5.5  175s\n",
      " 2617389 1813273 811013.608   63   20 411530.703 836375.946   103%   5.5  180s\n",
      " 2693117 1865482 infeasible  113      411530.703 836375.791   103%   5.5  185s\n",
      " 2772751 1924177 593299.156  102   16 411530.703 836375.644   103%   5.5  190s\n",
      " 2837297 1969289     cutoff   83      411530.703 836375.518   103%   5.5  195s\n",
      " 2905908 2015550 836370.104   77   14 411530.703 836375.398   103%   5.5  200s\n",
      " 2981394 2066668 577916.822  104   19 411530.703 836375.265   103%   5.5  205s\n",
      " 3055692 2117808     cutoff   83      411530.703 836375.132   103%   5.5  210s\n",
      " 3118158 2161500 823541.490   62   20 411530.703 836375.039   103%   5.5  215s\n",
      " 3177273 2205071 766248.485   65   24 411530.703 836374.948   103%   5.5  220s\n",
      " 3241326 2250085 591440.916  110   17 411530.703 836374.857   103%   5.5  225s\n",
      " 3297912 2288107     cutoff  132      411530.703 836374.775   103%   5.5  230s\n",
      " 3359422 2329124     cutoff   72      411530.703 836374.702   103%   5.5  235s\n",
      " 3435060 2381339     cutoff   71      411530.703 836374.602   103%   5.5  240s\n",
      " 3512498 2435257 827714.793   55   20 411530.703 836374.513   103%   5.5  245s\n",
      " 3585947 2486526 811569.039   63   20 411530.703 836374.432   103%   5.5  250s\n",
      " 3653077 2532411     cutoff  126      411530.703 836374.365   103%   5.5  255s\n",
      " 3714363 2573644 525568.001  110   12 411530.703 836374.303   103%   5.5  260s\n",
      " 3791329 2626887 493137.296   85   16 411530.703 836374.227   103%   5.5  265s\n",
      " 3859267 2674192 infeasible  121      411530.703 836374.164   103%   5.5  270s\n",
      " 3926788 2721494     cutoff  129      411530.703 836374.095   103%   5.5  275s\n",
      " 4000791 2775237 infeasible  141      411530.703 836374.020   103%   5.5  280s\n",
      " 4070200 2824363 836370.786   59   24 411530.703 836373.960   103%   5.5  285s\n",
      " 4134032 2866820 836370.415   78   20 411530.703 836373.911   103%   5.5  290s\n",
      " 4202036 2911142 673328.150  100   17 411530.703 836373.854   103%   5.5  295s\n",
      " 4268134 2954435 771842.062   67   24 411530.703 836373.802   103%   5.5  300s\n",
      "\n",
      "Explored 4268195 nodes (23664518 simplex iterations) in 300.00 seconds\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 6: 411531 390878 365784 ... 231816\n",
      "\n",
      "Time limit reached\n",
      "Best objective 4.115307025599e+05, best bound 8.363738023290e+05, gap 103.2348%\n"
     ]
    }
   ],
   "source": [
    "p_pooling = gp.Model(\"Pooling\")\n",
    "\n",
    "# Set global parameters\n",
    "p_pooling.params.nonConvex = 2\n",
    "p_pooling.params.timelimit = 5*60 # time limit of 5 minutes\n",
    "\n",
    "# Declare decision variables\n",
    "\n",
    "# flow\n",
    "ik = p_pooling.addVars(s2t, name=\"Source2Target\")\n",
    "ij = p_pooling.addVars(s2p, name=\"Source2Pool\")\n",
    "jk = p_pooling.addVars(p2t, name=\"Pool2Target\")\n",
    "ik[\"s1\",\"t2\"].ub = 750\n",
    "ik[\"s3\",\"t1\"].ub = 750\n",
    "# quality\n",
    "prop = p_pooling.addVars(pools, attrs, name=\"Proportion\")\n",
    "\n",
    "# Deploy constraint sets\n",
    "\n",
    "# 1. Flow conservation\n",
    "p_pooling.addConstrs((ij.sum('*',j) == jk.sum(j,'*') for j in pools),\n",
    "                     name=\"Flow_conservation\")\n",
    "# 2. Source capacity\n",
    "p_pooling.addConstrs((ij.sum(i,'*') + ik.sum(i,'*') <= supply[i] for i in sources),\n",
    "                     name=\"Source_capacity\")\n",
    "# 3. Pool capacity\n",
    "p_pooling.addConstrs((jk.sum(j,'*') <= cap[j] for j in pools),\n",
    "                     name=\"Pool_capacity\")\n",
    "# 4. Target demand\n",
    "p_pooling.addConstrs((ik.sum('*',k) + jk.sum('*',k) >= demand[k] for k in targets),\n",
    "                     name=\"Target_demand\")\n",
    "# 5. Pool concentration\n",
    "p_pooling.addConstrs((gp.quicksum(content[i][attr]*ij[i,j]\n",
    "                               for i in sources if (i,j) in s2p)\n",
    "                      == prop[j,attr]*jk.sum(j,'*') for j in pools for attr in attrs),\n",
    "                     name=\"Pool_concentration\")\n",
    "# 6.1 Target (min) tolerances\n",
    "p_pooling.addConstrs((gp.quicksum(content[i][attr]*ik[i,k]\n",
    "                               for i in sources if (i,k) in s2t)\n",
    "                      + gp.quicksum(prop[j,attr]*jk[j,k]\n",
    "                                 for j in pools if (j,k) in p2t)\n",
    "                      >= min_tol[k][attr]*(ik.sum('*',k) + jk.sum('*',k))\n",
    "                      for k in targets for attr in min_tol[k].keys()),\n",
    "                     name=\"Target_min_tolerances\")\n",
    "# 6.2 Target (max) tolerances\n",
    "p_pooling.addConstrs((gp.quicksum(content[i][attr]*ik[i,k]\n",
    "                               for i in sources if (i,k) in s2t)\n",
    "                      + gp.quicksum(prop[j,attr]*jk[j,k]\n",
    "                                 for j in pools if (j,k) in p2t)\n",
    "                      <= max_tol[k][attr]*(ik.sum('*',k) + jk.sum('*',k))\n",
    "                      for k in targets for attr in max_tol[k].keys()),\n",
    "                     name=\"Target_max_tolerances\")\n",
    "\n",
    "# Deploy Objective Function\n",
    "\n",
    "# 0. Total profit\n",
    "obj = gp.quicksum(price[k]*(ik.sum('*',k) + jk.sum('*',k))\n",
    "               for k in targets) \\\n",
    "- gp.quicksum(cost[i]*(ij.sum(i,'*') + ik.sum(i,'*'))\n",
    "           for i in sources)\n",
    "p_pooling.setObjective(obj, GRB.MAXIMIZE)\n",
    "\n",
    "# Find the optimal solution\n",
    "p_pooling.optimize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The P-formulation for this instance has:\n",
    "\n",
    "- 24 decision variables.\n",
    "- 10 linear constraints.\n",
    "- 20 bilinear constraints.\n",
    "- a linear objective function.\n",
    "\n",
    "As can be seen, we observe a gap of 103.23% after reaching the time limit of five minutes (at this point, the incumbent solution induces a total profit of 411,530.70 USD). In fact, even after 20 minutes the solver does not make much progress in closing the gap.\n",
    "\n",
    "### Q-Formulation (proportion)\n",
    "\n",
    "Let's now see how the q-formulation model performs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Changed value of parameter nonConvex to 2\n",
      "   Prev: -1  Min: -1  Max: 2  Default: -1\n",
      "Changed value of parameter timelimit to 300.0\n",
      "   Prev: inf  Min: 0.0  Max: inf  Default: inf\n",
      "Gurobi Optimizer version 9.1.0 build v9.1.0rc0 (win64)\n",
      "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n",
      "Optimize a model with 7 rows, 16 columns and 22 nonzeros\n",
      "Model fingerprint: 0x5737f095\n",
      "Model has 18 quadratic objective terms\n",
      "Model has 15 quadratic constraints\n",
      "Coefficient statistics:\n",
      "  Matrix range     [1e+00, 1e+00]\n",
      "  QMatrix range    [6e-01, 1e+02]\n",
      "  QLMatrix range   [1e-02, 1e+02]\n",
      "  Objective range  [9e+01, 2e+02]\n",
      "  QObjective range [1e+02, 6e+02]\n",
      "  Bounds range     [1e+00, 8e+02]\n",
      "  RHS range        [1e+00, 2e+03]\n",
      "  QRHS range       [5e+02, 2e+04]\n",
      "\n",
      "Continuous model is non-convex -- solving as a MIP.\n",
      "\n",
      "Presolve time: 0.00s\n",
      "Presolved: 95 rows, 35 columns, 315 nonzeros\n",
      "Presolved model has 18 bilinear constraint(s)\n",
      "Variable types: 35 continuous, 0 integer (0 binary)\n",
      "\n",
      "Root relaxation: objective 1.810743e+06, 41 iterations, 0.00 seconds\n",
      "\n",
      "    Nodes    |    Current Node    |     Objective Bounds      |     Work\n",
      " Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time\n",
      "\n",
      "     0     0 1810743.26    0   13          - 1810743.26      -     -    0s\n",
      "     0     0 865587.763    0   13          - 865587.763      -     -    0s\n",
      "     0     0 595135.123    0   13          - 595135.123      -     -    0s\n",
      "     0     0 564193.525    0   15          - 564193.525      -     -    0s\n",
      "     0     0 556780.563    0   15          - 556780.563      -     -    0s\n",
      "     0     0 549889.932    0   15          - 549889.932      -     -    0s\n",
      "     0     0 549889.932    0   18          - 549889.932      -     -    0s\n",
      "     0     0 549889.932    0   18          - 549889.932      -     -    0s\n",
      "     0     2 549889.932    0   18          - 549889.932      -     -    0s\n",
      "*   56    56               7    327332.37947 465046.855  42.1%   6.8    0s\n",
      "H   59    62                    433679.50783 465046.855  7.23%   6.5    0s\n",
      "H  107    86                    439182.58928 465046.855  5.89%   4.9    0s\n",
      "*  251     0              25    439182.58937 439182.589  0.00%   3.0    0s\n",
      "\n",
      "Cutting planes:\n",
      "  RLT: 18\n",
      "\n",
      "Explored 263 nodes (860 simplex iterations) in 0.11 seconds\n",
      "Thread count was 8 (of 8 available processors)\n",
      "\n",
      "Solution count 3: 439183 433680 327332 \n",
      "\n",
      "Optimal solution found (tolerance 1.00e-04)\n",
      "Warning: max constraint violation (1.3855e-06) exceeds tolerance\n",
      "Best objective 4.391825892869e+05, best bound 4.391825893721e+05, gap 0.0000%\n"
     ]
    }
   ],
   "source": [
    "q_pooling = gp.Model(\"Pooling\")\n",
    "\n",
    "# Set global parameters\n",
    "q_pooling.params.nonConvex = 2\n",
    "q_pooling.params.timelimit = 5*60\n",
    "\n",
    "# Declare decision variables\n",
    "\n",
    "# flow\n",
    "ik = q_pooling.addVars(s2t, name=\"Source2Target\")\n",
    "jk = q_pooling.addVars(p2t, name=\"Pool2Target\")\n",
    "ik[\"s1\",\"t2\"].ub = 750\n",
    "ik[\"s3\",\"t1\"].ub = 750\n",
    "# proportion\n",
    "p_ij = q_pooling.addVars(s2p, ub=1.0, name=\"Prop_Source2Pool\")\n",
    "\n",
    "# Deploy constraint sets\n",
    "\n",
    "# 1. Source capacity\n",
    "q_pooling.addConstrs((gp.quicksum(p_ij[i,j]*jk.sum(j,'*')\n",
    "                               for j in pools if (i,j) in s2p)\n",
    "                      + ik.sum(i,'*') <= supply[i] for i in sources),\n",
    "                     name=\"Source_capacity\")\n",
    "# 2. Pool capacity\n",
    "q_pooling.addConstrs((jk.sum(j,'*') <= cap[j] for j in pools),\n",
    "                     name=\"Pool_capacity\")\n",
    "# 3. Target demand\n",
    "q_pooling.addConstrs((ik.sum('*',k) + jk.sum('*',k) >= demand[k] for k in targets),\n",
    "                     name=\"Target_demand\")\n",
    "# 4. Pool inflow\n",
    "q_pooling.addConstrs((p_ij.sum('*',j) == 1 for j in pools),\n",
    "                     name=\"Pool_inflow\")\n",
    "# 5.1 Target (min) tolerances\n",
    "q_pooling.addConstrs((gp.quicksum(content[i][attr]*ik[i,k]\n",
    "                               for i in sources if (i,k) in s2t)\n",
    "                      + gp.quicksum(content[i][attr]*p_ij[i,j]*jk[j,k]\n",
    "                                 for (i,j) in s2p if (j,k) in p2t)\n",
    "                      >= min_tol[k][attr]*(ik.sum('*',k) + jk.sum('*',k))\n",
    "                      for k in targets for attr in min_tol[k].keys()),\n",
    "                     name=\"Target_min_tolerances\")\n",
    "# 5.2 Target (max) tolerances\n",
    "q_pooling.addConstrs((gp.quicksum(content[i][attr]*ik[i,k]\n",
    "                               for i in sources if (i,k) in s2t)\n",
    "                      + gp.quicksum(content[i][attr]*p_ij[i,j]*jk[j,k]\n",
    "                                 for (i,j) in s2p if (j,k) in p2t)\n",
    "                      <= max_tol[k][attr]*(ik.sum('*',k) + jk.sum('*',k))\n",
    "                      for k in targets for attr in max_tol[k].keys()),\n",
    "                     name=\"Target_max_tolerances\")\n",
    "\n",
    "# Deploy Objective Function\n",
    "\n",
    "# 0. Total profit\n",
    "obj = gp.quicksum(price[k]*(ik.sum('*',k) + jk.sum('*',k))\n",
    "               for k in targets) \\\n",
    "- gp.quicksum(cost[i]*(gp.quicksum(p_ij[i,j]*jk.sum(j,'*')\n",
    "                             for j in pools if (i,j) in s2p)\n",
    "                    + ik.sum(i,'*'))\n",
    "           for i in sources)\n",
    "q_pooling.setObjective(obj, GRB.MAXIMIZE)\n",
    "\n",
    "# Find the optimal solution\n",
    "q_pooling.optimize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Q-formulation for this instance has:\n",
    "\n",
    "- 16 decision variables.\n",
    "- 7 linear constraints.\n",
    "- 15 bilinear constraints.\n",
    "- a bilinear objective function.\n",
    "\n",
    "Notice it has fewer decision variables and also fewer bilinear constraints. Now Gurobi was able to find the optimal solution of 439,182.59 USD in less than one second.\n",
    "\n",
    "---\n",
    "## Analysis\n",
    "\n",
    "Let's see the optimal flows found:\n",
    "\n",
    "### Flow from Sources to Targets\n",
    "The following table determines the flows from source nodes to target nodes. For example, from source node s2 to target node t1 there is a flow of 966.7$\\times 10^{2}$ bbl."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>t1</th>\n",
       "      <th>t2</th>\n",
       "      <th>t3</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>s1</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>s2</th>\n",
       "      <td>966.7</td>\n",
       "      <td>0.0</td>\n",
       "      <td>200.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>s3</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       t1   t2     t3\n",
       "s1    0.0  0.0    0.0\n",
       "s2  966.7  0.0  200.0\n",
       "s3    0.0  0.0    0.0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rows = sources.copy()\n",
    "columns = targets.copy()\n",
    "s2t_plan = pd.DataFrame(columns=columns, index=rows, data=0.0)\n",
    "\n",
    "for source, target in ik.keys():\n",
    "    if (abs(ik[source, target].x) > 1e-6):\n",
    "        s2t_plan.loc[source, target] = np.round(ik[source, target].x, 1)\n",
    "s2t_plan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Flow from Pools to Targets\n",
    "The following table defines the flows from pool nodes to target nodes. For example, from pool node p1 to target node t1 there is a flow of 92.8$\\times 10^{2}$ bbl."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>t1</th>\n",
       "      <th>t2</th>\n",
       "      <th>t3</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>p1</th>\n",
       "      <td>92.8</td>\n",
       "      <td>990.6</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>p2</th>\n",
       "      <td>1450.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>300.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "        t1     t2     t3\n",
       "p1    92.8  990.6    0.0\n",
       "p2  1450.0    0.0  300.0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rows = pools.copy()\n",
    "columns = targets.copy()\n",
    "p2t_plan = pd.DataFrame(columns=columns, index=rows, data=0.0)\n",
    "\n",
    "for pool, target in jk.keys():\n",
    "    if (abs(jk[pool, target].x) > 1e-6):\n",
    "        p2t_plan.loc[pool, target] = np.round(jk[pool, target].x, 1)\n",
    "p2t_plan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Flow from Sources to Pools\n",
    "The following table shows the flows from source nodes to pool nodes. For example, from source node s2 to pool node p1  there is a flow of 258.3$\\times 10^{2}$ bbl."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>p1</th>\n",
       "      <th>p2</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>s1</th>\n",
       "      <td>325.0</td>\n",
       "      <td>1750.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>s2</th>\n",
       "      <td>258.3</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>s3</th>\n",
       "      <td>500.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       p1      p2\n",
       "s1  325.0  1750.0\n",
       "s2  258.3     0.0\n",
       "s3  500.0     0.0"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rows = sources.copy()\n",
    "columns = pools.copy()\n",
    "s2p_plan = pd.DataFrame(columns=columns, index=rows, data=0.0)\n",
    "\n",
    "for source, pool in p_ij.keys():\n",
    "    if (abs(p_ij[source, pool].x) > 1e-6):\n",
    "        flow = p_ij[source, pool].x * p2t_plan.loc[pool,:].sum()\n",
    "        s2p_plan.loc[source, pool] = np.round(flow, 1)\n",
    "s2p_plan"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "## Conclusions\n",
    "\n",
    "This notebook showed how easy it is to solve Bilinear Programs using Gurobi. It also highlighted the dramatic difference in performance of alternative formulations when solving challenging problems, such as the Standard Pooling Problem. It is thus of utmost importance to analyze carefully the context of the problem at hand, and to weigh the pros and cons of alternative models. \n",
    "\n",
    "---\n",
    "<a id='references'></a>\n",
    "## References\n",
    "\n",
    "1. Alfaki, M. (2012). Models and solution methods for the pooling problem.\n",
    "2. Audet, C., Brimberg, J., Hansen, P., Digabel, S. L., & Mladenović, N. (2004). Pooling problem: Alternate formulations and solution methods. Management science, 50(6), 761-776.\n",
    "3. Dombrowski, J. (2015, June 07). McCormick envelopes. Retrieved from https://optimization.mccormick.northwestern.edu/index.php/McCormick_envelopes\n",
    "4. Gupte, A., Ahmed, S., Dey, S. S., & Cheon, M. S. (2017). Relaxations and discretizations for the pooling problem. Journal of Global Optimization, 67(3), 631-669.\n",
    "5. Haverly, C. A. (1978). Studies of the behavior of recursion for the pooling problem. Acm sigmap bulletin, (25), 19-28.\n",
    "6. Liberti, L. (2008). Introduction to global optimization. Ecole Polytechnique.\n",
    "7. Zhuang E. (2015, June 06). Spatial branch and bound method. Retrieved from\n",
    "https://optimization.mccormick.northwestern.edu/index.php/Spatial_branch_and_bound_method\n",
    "\n",
    "Copyright © 2020 Gurobi Optimization, LLC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}