# Haverly's Pooling Problem

## Objective and Prerequisites

One of the new features of Gurobi 9.0 is the addition of a bilinear solver, which enables finding the optimal solution of non-convex quadratic programming  problems (i.e. QPs, QCQPs, MIQPs, and MIQCQPs). This notebook will show you how to use this feature by tackling the Haverly's Pooling Problem.

To fully understand the content of this notebook, the reader should be familiar with the following:

- Python.
- Gurobi Python interface.
- Knowledge of building mathematical optimization models.

**Note:** You can download the repository containing this and other examples by clicking [here](https://github.com/Gurobi/modeling-examples/archive/master.zip). In order to run this Jupyter Notebook properly, you must have a Gurobi license. If you do not have one, you can request an [evaluation license](https://www.gurobi.com/downloads/request-an-evaluation-license/?utm_source=Github&utm_medium=website_JupyterME&utm_campaign=CommercialDataScience) as a *commercial user*, or download a [free license](https://www.gurobi.com/academia/academic-program-and-licenses/?utm_source=Github&utm_medium=website_JupyterME&utm_campaign=AcademicDataScience) as an *academic user*.

---
## Motivation

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 [3].

---
## Problem Description

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.

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) [2].

---
## Define Data

In this example, there are three sources of oil, each of which contains a different percentage of sulfur. Given the pooling layout, these materials are to be blended in such a way to create two separate oil products. These final products have requirements on the percentage of sulfur they can contain, as well as how much of the product can be produced.

Now, we define the data used in Haverly's 1978 paper:


```python
import numpy as np
import pandas as pd
from itertools import product
from gurobipy import *

attrs = {'sulfur'}

sources, cost, content = multidict({
    "A": [6, {'sulfur': 3}],
    "B": [16, {'sulfur': 1}],
    "C": [10, {'sulfur': 2}]
})

targets, price, demand, max_tol = multidict({
    "X": [9, 100, {'sulfur': 2.5}],
    "Y": [15, 200, {'sulfur': 1.5}]
})

pools = {"P"}


s2p = {("A", "P"),
       ("B", "P")}
# The function `product` deploys the Cartesian product of elements in sets A and B
p2t = set(product(pools, targets))
s2t = {("C", "X"),
       ("C", "Y")}
```



---
## Mathematical Formulation

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.

A mathematical optimization model has five components, namely:

- Sets and indices.
- Parameters.
- Decision variables.
- Objective function(s).
- Constraints.

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) [1], [4] and [5].

We now present a Bilinear Program for the Standard Pooling Problem:

#### Sets and Indices

$G=(V,E)$: Directed graph.

$i,j \in V$: Set of nodes.

$(i,j) \in E \subset  V \times V$: Set of edges.

$N(i)^+ = \{j \in V \mid (i,j) \in E \}$: Set of successor nodes receiving outflow from node $i$.

$N(j)^- = \{i \in V \mid (i,j) \in E \}$: Set of predecessor nodes sending inflow to node $i$.

$k \in \text{Attrs}$: Set of attributes.

$s \in \text{Sources} \subset V$: Set of source nodes, i.e. $N(s)^-= \emptyset$.

$t \in \text{Targets} \subset V$: Set of target nodes, i.e. $N(t)^+= \emptyset$.

$p \in \text{Pools} = V \setminus (\text{Sources} \cup \text{Targets})$: Set of pools.

#### Parameters

$\text{Cost}_s \in \mathbb{R}^+$: Cost of acquiring one unit of raw material at source node $s$.

$\text{Content}_{s,k} \in \mathbb{R}^+$: Content of attribute $k$ in raw material at source node $s$.

$\text{Price}_t \in \mathbb{R}^+$: Price for selling one unit of final blend at target node $t$.

$\text{Demand}_t \in \mathbb{R}^+$: Minimum number of units required of final blend at target node $t$.

$\text{Max_tol}_{t,k} \in \mathbb{R}^+$: Maximum tolerance for attribute $k$ in final blend at target node $t$.

#### Decision Variables

$\text{flow}_{i,j} \in [0, \text{UB}_{i,j}]$: Flow from node $i$ to node $j$.

$\text{quality}_{p,k} \in \mathbb{R}^+$: Concentration of attribute $k$ at pool $p$.

#### Objective Function

- **Profit**: Maximize total profits.

\begin{equation}
\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}}}
\tag{0}
\end{equation}

#### Constraints


- **Flow conservation**: Total inflow of pool $p$ must be equal to its toal outflow (nothing is stored in them).

\begin{equation}
\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}
\tag{1}
\end{equation}

- **Target demand**: Total inflow of target $t$ cannot exceed maximum demand.

\begin{equation}
\sum_{i \in N(t)^-}{\text{flow}_{i,t}} \leq \text{Demand}_t \quad \forall t \in \text{Targets}
\tag{2}
\end{equation}

- **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).

\begin{equation}
\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}
\tag{3}
\end{equation}

- **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).

\begin{equation}
\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}
\tag{4}
\end{equation}

---

## Python Implementation

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 spatial Branch and Bound (sB&B).


```python
haverly = Model("Pooling")

# Set global parameters
haverly.params.nonConvex = 2

# Declare decision variables

# flow
ik = haverly.addVars(s2t, name="Source2Target")
ij = haverly.addVars(s2p, name="Source2Pool")
jk = haverly.addVars(p2t, name="Pool2Target")
# quality
prop = haverly.addVars(pools, attrs, name="Proportion")

# Deploy constraint sets

# 1. Flow conservation
haverly.addConstrs((ij.sum('*',j) == jk.sum(j,'*') for j in pools),
                     name="Flow_conservation")
# 2. Target demand
haverly.addConstrs((ik.sum('*',k) + jk.sum('*',k) <= demand[k] for k in targets),
                     name="Target_demand")
# 3. Pool concentration
haverly.addConstrs((quicksum(content[i][attr]*ij[i,j]
                               for i in sources if (i,j) in s2p)
                      == prop[j,attr]*jk.sum(j,'*') for j in pools for attr in attrs),
                     name="Pool_concentration")
# 4. Target (max) tolerances
haverly.addConstrs((quicksum(content[i][attr]*ik[i,k]
                               for i in sources if (i,k) in s2t)
                      + quicksum(prop[j,attr]*jk[j,k]
                                 for j in pools if (j,k) in p2t)
                      <= max_tol[k][attr]*(ik.sum('*',k) + jk.sum('*',k))
                      for k in targets for attr in max_tol[k].keys()),
                     name="Target_max_tolerances")

# Deploy Objective Function

# 0. Total profit
obj = quicksum(price[k]*(ik.sum('*',k) + jk.sum('*',k))
               for k in targets) \
- quicksum(cost[i]*(ij.sum(i,'*') + ik.sum(i,'*'))
           for i in sources)
haverly.setObjective(obj, GRB.MAXIMIZE)

# Find the optimal solution
haverly.optimize()
```

    Using license file /Users/orojuan/gurobi.lic
    Set parameter TokenServer to value Juans-MacBook-Pro-3.local
    Changed value of parameter nonConvex to 2
       Prev: -1  Min: -1  Max: 2  Default: -1
    Gurobi Optimizer version 9.0.0 build v9.0.0rc2 (mac64)
    Optimize a model with 3 rows, 7 columns and 8 nonzeros
    Model fingerprint: 0x777f4d01
    Model has 3 quadratic constraints
    Coefficient statistics:
      Matrix range     [1e+00, 1e+00]
      QMatrix range    [1e+00, 1e+00]
      QLMatrix range   [5e-01, 3e+00]
      Objective range  [1e+00, 2e+01]
      Bounds range     [0e+00, 0e+00]
      RHS range        [1e+02, 2e+02]
    
    Continuous model is non-convex -- solving as a MIP.
    
    Found heuristic solution: objective -0.0000000
    Presolve time: 0.00s
    Presolved: 16 rows, 10 columns, 34 nonzeros
    Presolved model has 4 bilinear constraint(s)
    Variable types: 10 continuous, 0 integer (0 binary)
    
    Root relaxation: objective 2.100000e+03, 4 iterations, 0.00 seconds
    
        Nodes    |    Current Node    |     Objective Bounds      |     Work
     Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
    
         0     0 2100.00000    0    3   -0.00000 2100.00000      -     -    0s
         0     0 2100.00000    0    3   -0.00000 2100.00000      -     -    0s
         0     0 1937.83784    0    3   -0.00000 1937.83784      -     -    0s
         0     0 1921.35002    0    4   -0.00000 1921.35002      -     -    0s
         0     0 1921.23235    0    4   -0.00000 1921.23235      -     -    0s
         0     0 1921.17768    0    4   -0.00000 1921.17768      -     -    0s
         0     2 1921.17768    0    4   -0.00000 1921.17768      -     -    0s
    *    6     6               2     321.5204987 1882.20468   485%   3.0    0s
    *   12     4               3     346.2144981 1846.72927   433%   2.4    0s
    *   14     4               4     348.0276050 1846.72927   431%   3.1    0s
    *   22    10               5     355.0833492 1778.85855   401%   3.0    0s
    *   23     5               5     363.4227077 1778.85855   389%   2.9    0s
    *   24     5               6     363.6955948 1778.85855   389%   2.9    0s
    *   25     5               6     370.7666965 1778.85855   380%   2.8    0s
    *   26     5               7     372.0510296 1778.85855   378%   2.8    0s
    *   32     2               7     377.9877375 1746.46493   362%   2.5    0s
    *   34     2               8     380.1494347 1746.46493   359%   2.4    0s
    *   36     4               8     383.8453897 1746.46493   355%   2.5    0s
    *   40     5               9     387.9905769 1715.10036   342%   2.5    0s
    *   41     5               9     390.8479187 1715.10036   339%   2.4    0s
    *   45     5              10     395.5742065 1684.76584   326%   2.4    0s
    *   47     5              10     397.5729772 1684.76584   324%   2.3    0s
    *   52     7              12     399.9188266 1684.76584   321%   2.3    0s
    *   60     2              16     399.9999622 1655.46244   314%   2.2    0s
    *   70     2              18     400.0000031 1627.19132   307%   2.1    0s
    
    Cutting planes:
      RLT: 3
    
    Explored 105 nodes (257 simplex iterations) in 0.06 seconds
    Thread count was 8 (of 8 available processors)
    
    Solution count 10: 400 400 399.919 ... 377.988
    
    Optimal solution found (tolerance 1.00e-04)
    Best objective 4.000000031377e+02, best bound 4.000000031377e+02, gap 0.0000%


### Analysis

Let's see the optimal flows found:


```python
def _print_table(rows, columns, variables):
    table = pd.DataFrame(columns=columns, index=rows, data=0.0)
    for row, col in variables.keys():
        value = variables[row, col].x
        if abs(value) > 1e-6:
            table.loc[row, col] = np.round(value, 1)
    print(table)
    
def print_solution():
    print("Flows from Sources to Targets")
    _print_table(sources.copy(), targets.copy(), ik)
    print("Flows from Pools to Targets")
    _print_table(pools.copy(), targets.copy(), jk)
    print("Flows from Sources to Pools")
    _print_table(sources.copy(), pools.copy(), ij)
```


```python
print_solution()
```

    Flows from Sources to Targets
         X      Y
    A  0.0    0.0
    B  0.0    0.0
    C  0.0  100.0
    Flows from Pools to Targets
         X      Y
    P  0.0  100.0
    Flows from Sources to Pools
           P
    A    0.0
    B  100.0
    C    0.0




---
## Changing the Data 

We now consider how the optimal solution changes by modifying some parameters in the model.

### Increasing Demand

The maximum demand of product $\text{X}$ increased from 100 to 600:


```python
dem_x = haverly.getConstrByName("Target_demand[X]")
dem_x.setAttr("rhs", 600) 
haverly.optimize()
```

    Gurobi Optimizer version 9.0.0 build v9.0.0rc2 (mac64)
    Optimize a model with 3 rows, 7 columns and 8 nonzeros
    Model fingerprint: 0xf2ff1488
    Model has 3 quadratic constraints
    Variable types: 7 continuous, 0 integer (0 binary)
    Coefficient statistics:
      Matrix range     [1e+00, 1e+00]
      QMatrix range    [1e+00, 1e+00]
      QLMatrix range   [5e-01, 3e+00]
      Objective range  [1e+00, 2e+01]
      Bounds range     [0e+00, 0e+00]
      RHS range        [2e+02, 6e+02]
    
    MIP start from previous solve produced solution with objective 324 (0.01s)
    MIP start from previous solve produced solution with objective 526.694 (0.02s)
    MIP start from previous solve produced solution with objective 599.488 (0.02s)
    MIP start from previous solve produced solution with objective 599.976 (0.02s)
    MIP start from previous solve produced solution with objective 600 (0.02s)
    MIP start from previous solve produced solution with objective 600 (0.03s)
    Loaded MIP start from previous solve with objective 600
    
    Presolve time: 0.00s
    Presolved: 16 rows, 10 columns, 34 nonzeros
    Presolved model has 4 bilinear constraint(s)
    Variable types: 10 continuous, 0 integer (0 binary)
    
    Root relaxation: objective 3.600000e+03, 4 iterations, 0.00 seconds
    
        Nodes    |    Current Node    |     Objective Bounds      |     Work
     Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
    
         0     0 3600.00000    0    3  600.00000 3600.00000   500%     -    0s
         0     0 3600.00000    0    3  600.00000 3600.00000   500%     -    0s
         0     0 3381.81818    0    4  600.00000 3381.81818   464%     -    0s
         0     0 3230.76924    0    3  600.00000 3230.76924   438%     -    0s
         0     0 3217.40285    0    4  600.00000 3217.40285   436%     -    0s
         0     0 3207.95528    0    4  600.00000 3207.95528   435%     -    0s
         0     2 3207.95528    0    4  600.00000 3207.95528   435%     -    0s
    *  126     5              54     600.0000064  600.70747  0.12%   3.5    0s
    
    Cutting planes:
      RLT: 3
    
    Explored 141 nodes (472 simplex iterations) in 0.06 seconds
    Thread count was 8 (of 8 available processors)
    
    Solution count 6: 600 600 599.976 ... 324
    
    Optimal solution found (tolerance 1.00e-04)
    Best objective 6.000000063584e+02, best bound 6.000000063584e+02, gap 0.0000%



```python
print_solution()
```

    Flows from Sources to Targets
           X    Y
    A    0.0  0.0
    B    0.0  0.0
    C  300.0  0.0
    Flows from Pools to Targets
           X    Y
    P  300.0  0.0
    Flows from Sources to Pools
           P
    A  300.0
    B    0.0
    C    0.0




### Decreasing Cost

The price of extracting crude oil from source $\text{B}$ decreased from $\$16$ to $\$13$:


```python
dem_x.setAttr("rhs", 100) # reinstate the model to its initial state
ij["B", "P"].obj = -13 # the coefficient is negative since we're modifying a cost
haverly.optimize()
```

    Gurobi Optimizer version 9.0.0 build v9.0.0rc2 (mac64)
    Optimize a model with 3 rows, 7 columns and 8 nonzeros
    Model fingerprint: 0x94aace23
    Model has 3 quadratic constraints
    Variable types: 7 continuous, 0 integer (0 binary)
    Coefficient statistics:
      Matrix range     [1e+00, 1e+00]
      QMatrix range    [1e+00, 1e+00]
      QLMatrix range   [5e-01, 3e+00]
      Objective range  [1e+00, 2e+01]
      Bounds range     [0e+00, 0e+00]
      RHS range        [1e+02, 2e+02]
    
    MIP start from previous solve produced solution with objective -14.0625 (0.01s)
    MIP start from previous solve produced solution with objective -0 (0.02s)
    MIP start from previous solve produced solution with objective 1.88255 (0.03s)
    MIP start from previous solve produced solution with objective 721.662 (0.03s)
    MIP start from previous solve produced solution with objective 726.254 (0.03s)
    Loaded MIP start from previous solve with objective 726.254
    
    Presolve time: 0.00s
    Presolved: 16 rows, 10 columns, 34 nonzeros
    Presolved model has 4 bilinear constraint(s)
    Variable types: 10 continuous, 0 integer (0 binary)
    
    Root relaxation: objective 2.100000e+03, 4 iterations, 0.00 seconds
    
        Nodes    |    Current Node    |     Objective Bounds      |     Work
     Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time
    
         0     0 2100.00000    0    3  726.25404 2100.00000   189%     -    0s
         0     0 2100.00000    0    3  726.25404 2100.00000   189%     -    0s
         0     0 1935.80379    0    3  726.25404 1935.80379   167%     -    0s
         0     0 1919.29441    0    4  726.25404 1919.29441   164%     -    0s
         0     0 1918.94856    0    4  726.25404 1918.94856   164%     -    0s
         0     0 1918.78695    0    4  726.25404 1918.78695   164%     -    0s
         0     2 1918.78695    0    4  726.25404 1918.78695   164%     -    0s
    *    4     6               2     740.2106233 1900.53786   157%   2.8    0s
    *   10     4               3     750.0000000 1844.51703   146%   2.6    0s
    
    Cutting planes:
      RLT: 3
    
    Explored 49 nodes (143 simplex iterations) in 0.08 seconds
    Thread count was 8 (of 8 available processors)
    
    Solution count 7: 750 740.211 726.254 ... -14.0625
    
    Optimal solution found (tolerance 1.00e-04)
    Best objective 7.500000000000e+02, best bound 7.500000000000e+02, gap 0.0000%



```python
print_solution()
```

    Flows from Sources to Targets
         X    Y
    A  0.0  0.0
    B  0.0  0.0
    C  0.0  0.0
    Flows from Pools to Targets
         X      Y
    P  0.0  200.0
    Flows from Sources to Pools
           P
    A   50.0
    B  150.0
    C    0.0




---
## Conclusion

This notebook showed how easy it is to solve Bilinear Programs using Gurobi.

The Haverly's pooling problem is indeed a simple instance of the Standard Pooling Problem, as it considers only one attribute. It was modeled using what the literature calls the P-Formulation, where the number of bilinear terms is proportional to the number of attributes. For intances with more specifications, it may be worth considering using alternative formulations —such as the Q-Formulation— to improve the performance of the optimization process. The Jupyter Notebook `Standard Pooling Problem` presents and compares both formulations with a more challenging scenario.

---

## References

1. Dombrowski, J. (2015, June 07). McCormick envelopes. Retrieved from https://optimization.mccormick.northwestern.edu/index.php/McCormick_envelopes
2. 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.
3. Haverly, C. A. (1978). Studies of the behavior of recursion for the pooling problem. Acm sigmap bulletin, (25), 19-28.
4. Liberti, L. (2008). Introduction to global optimization. Ecole Polytechnique.
5. Zhuang E. (2015, June 06). Spatial branch and bound method. Retrieved from
https://optimization.mccormick.northwestern.edu/index.php/Spatial_branch_and_bound_method

Copyright © 2019 Gurobi Optimization, LLC

63ec521c2c160a2316e904a0ac8ff705b095bf9d
ipynb
Jupyter Notebook
haverly/haverly.ipynb
Zhong-HY/modeling-examples
66355938dde764a320f0a6a70a9e5a69c696a18f
Apache-2.0
haverly/haverly.ipynb
Zhong-HY/modeling-examples
66355938dde764a320f0a6a70a9e5a69c696a18f
Apache-2.0
haverly/haverly.ipynb
Zhong-HY/modeling-examples
66355938dde764a320f0a6a70a9e5a69c696a18f
Apache-2.0
Qwen/Qwen-72B
1. YES
2. YES