{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# The pooling problem\n",
"\n",
"In the pooling problem, source materials with are mixed in intermediate pools to create a number of output products. The output products are constrained by the quality of a certain attribute. In this notebook, we consider the small pooling example of Haverly:\n",
"\n",
"C. A. Haverly, *Studies of the behavior of the recursion for the pooling problem*, SIGMAP Bulletin, 25 (1978), pp. 19–28. [pdf](https://www.eq.uc.pt/~nuno/eps/Edicoes_anteriores/2012-13/aula08/Haverly%20(1978).pdf)\n",
"\n",
"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.\n",
"\n",
"The problem's difficulty stems from the fact that it is a non-convex problem, as some of the constraints contain bilinear terms. Gurobi introduced support for non-convex quadratically constrained models with release 9.0. This notebook demonstrates how the example pooling problem can be defined, modeled, and solved with Gurobi/Python.\n",
"\n",
"Most nonlinear solvers look for a local solution to this problem. In contrast, Gurobi finds the global optimum.\n",
"\n",
"We begin by importing the `gurobipy` package."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from gurobipy import *"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Define data\n",
"\n",
"Now, we define the data used in Haverly's 1978 paper."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Define the crude oil sources, cost per unit, and sulfur percentage (i.e., quality)\n",
"sources, costs, sulfur = multidict({ \n",
" \"A\" : [ 6, 3 ],\n",
" \"B\" : [ 16, 1 ],\n",
" \"C\" : [ 10, 2 ] \n",
"})\n",
"\n",
"# Define the crude oil sinks, market price, maximum sulfur percentage (quality), and maximum satisfiable demand\n",
"sinks, price, max_sulfur, demand = multidict({\n",
" \"X\" : [ 9, 2.5, 100 ],\n",
" \"Y\" : [ 15, 1.5, 200 ]\n",
"})\n",
"\n",
"# Define intermediate pools\n",
"pools = [ \"P\" ]\n",
"\n",
"# Node connections (i.e., graph edges)\n",
"connections = tuplelist([\n",
" (\"A\", \"P\"), (\"B\", \"P\"), # Connections to pool\n",
" (\"P\", \"X\"), (\"P\", \"Y\"), # Connections from pool to sink\n",
" (\"C\", \"X\"), (\"C\", \"Y\") # Direct connections\n",
"])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Visualizing the problem\n",
"\n",
"We can draw the directed graph with the `graphviz` package (see [here](https://graphviz.readthedocs.io/en/stable/manual.html))."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from graphviz import Digraph\n",
"dot = Digraph(graph_attr={'rankdir': 'LR'}, node_attr={'fontcolor':'blue', 'shape' : 'box'})\n",
"for pool in pools:\n",
" dot.node(pool, shape='doublecircle')\n",
"for sink in sinks:\n",
" dot.node(sink, \"{:s} @${:d}\\n\\u2264 {:0.1f}%\\n \\u2264 {:d}\".format(sink, price[sink], max_sulfur[sink], demand[sink]), shape='invhouse')\n",
"for source in sources:\n",
" dot.node(source, \"{:s}\\n {:d}% @${:d}\".format(source, sulfur[source], costs[source]))\n",
"for connection in connections:\n",
" dot.edge(*connection)\n",
"dot"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Mathematical formulation\n",
"\n",
"$N = \\{ A, B, C, P, X, Y \\}$ be the set of nodes in the graph, and let $\\mathcal{A}$ be the set of arcs. We create two sets of variables:\n",
"\n",
"* *flow* variables $x_{ij}$ for the arcs $(i,j) \\in \\mathcal{A}$, and\n",
"* *quality* (sulfur percentage) variables $y_i$ for the nodes $i \\in N$.\n",
"\n",
"For source nodes $i \\in S := \\{ A, B, C \\}$, let $f_i$ be the fixed sulfur percentage and let $c_i$ be the unit cost for obtaining oil from that source. For sink nodes $i \\in T := \\{ X, Y \\}$, let $u_i$ be the maximum sulfur percentage, let $p_i$ be the unit selling price, and let $d_i$ be the maximum demand. We can model this problem as follows:\n",
"\n",
" \n",
"\n",
"$$ \n",
"\\begin{alignat*}{4}\n",
"\\max_{x, y}\\ && \\sum_{i \\in T} \\sum_{(j,i) \\in \\mathcal{A}} p_i x_{ji} &- \\sum_{i \\in S} \\sum_{(i,j) \\in \\mathcal{A}} c_i x_{ij} \\\\\n",
"\\textrm{s.t. } && \\sum_{(j,i) \\in \\mathcal{A}} x_{ji} &= \\sum_{(i,j) \\in \\mathcal{A}} x_{ij} && \\forall i \\in N \\setminus (S \\cup T) \\quad && \\textrm{(flow balance)} \\\\\n",
"&& y_i \\sum_{(j,i) \\in \\mathcal{A}} x_{ji} &= \\sum_{(j,i) \\in \\mathcal{A}} y_j x_{ji} && \\forall i \\in N && \\textrm{(define quality)} \\\\\n",
"&& y_i &= f_i && \\forall i \\in S \\quad && \\textrm{(fix source quality)} \\\\\n",
"&& y_i &\\leq u_i && \\forall i \\in T && \\textrm{(bound sink quality)} \\\\\n",
"&& \\sum_{(j,i) \\in \\mathcal{A}} x_{ji} &\\leq d_i && \\forall i \\in T && \\textrm{(maximum demand)} \\\\\n",
"&& x, y &\\geq 0\n",
"\\end{alignat*}\n",
"$$\n",
"\n",
"\n",
"We can model this problem using `gurobipy`."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Variables and bounds"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def define_variables(m):\n",
" \n",
" # Total flow per connection\n",
" flow = m.addVars(connections, name=\"flow\")\n",
"\n",
" # Quality (percentage of sulfur after blending all inflows)\n",
" quality = m.addVars(sources + pools + sinks, name=\"quality\")\n",
"\n",
" # Fix quality for all sources\n",
" for source in sources:\n",
" quality[source].LB = sulfur[source]\n",
" quality[source].UB = sulfur[source]\n",
"\n",
" # Bound minimum quality for all sinks\n",
" for sink in sinks:\n",
" quality[sink].UB = max_sulfur[sink]\n",
" \n",
" return flow, quality"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Objective function (profit)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def define_objective(m, flow, quality):\n",
" \n",
" # Build objective function\n",
" total_costs = sum(costs[i] * flow[i,j] for s in sources for i,j in connections.select(s, '*'))\n",
" total_revenue = sum(price[j] * flow[i,j] for t in sinks for i,j in connections.select('*', t))\n",
" m.setObjective(total_revenue - total_costs, GRB.MAXIMIZE)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Constraints"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def define_constraints(m, flow, quality):\n",
" \n",
" # Add balance constraints\n",
" m.addConstrs((flow.sum('*', pool) == flow.sum(pool, '*') for pool in pools), name=\"balance\")\n",
"\n",
" # Add demand constraints\n",
" m.addConstrs((flow.sum('*', j) <= demand[j] for j in sinks), name=\"demand\")\n",
"\n",
" # Add quality constraints\n",
" m.addConstrs((quality[node] * flow.sum('*', node)\n",
" == sum(quality[i] * flow[i, j] for i,j in connections.select('*', node))\n",
" for node in pools + sinks), name=\"quality\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, we let Gurobi build and solve the model."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"m = Model()\n",
"\n",
"flow, quality = define_variables(m)\n",
"define_objective(m, flow, quality)\n",
"define_constraints(m, flow, quality)\n",
"\n",
"m.params.NonConvex=2\n",
"m.optimize()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Visualize flows"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def draw_solution(m, flow, quality):\n",
" \n",
" dot = Digraph(graph_attr={'rankdir': 'LR'}, node_attr={'fontcolor':'blue', 'shape' : 'box'})\n",
" for pool in pools: dot.node(pool, shape='doublecircle')\n",
" for sink in sinks: dot.node(sink, label=\"<%s
(%g)>\" % (sink, sum([f.X for f in flow.select('*', sink)])),\n",
" shape='invhouse')\n",
" for connection in connections:\n",
" dot.edge(*connection, label=(\"%g @ %g%%\" % (flow[connection].X, quality[connection[0]].X)) if flow[connection].X > 0 else '', \n",
" color='red' if flow[connection].X > 0 else 'lightgray',\n",
" fontname='courier', fontsize='10')\n",
" return dot\n",
"\n",
"draw_solution(m, flow, quality)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Changing the data \n",
"\n",
"### Increasing demand\n",
"\n",
"We increase the maximum demand of product X from 100 to 600 and solve the model again."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"demand['X'] = 600\n",
"m2 = Model()\n",
"\n",
"flow2, quality2 = define_variables(m2)\n",
"define_objective(m2, flow2, quality2)\n",
"define_constraints(m2, flow2, quality2)\n",
"\n",
"m2.params.NonConvex=2\n",
"m2.optimize()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"# Draw solution\n",
"draw_solution(m2, flow2, quality2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Changing source price\n",
"\n",
"We lower the price of extracting crude oil from source B from $\\$16$ to $\\$13$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"costs['B'] = 13\n",
"\n",
"m3 = Model()\n",
"\n",
"flow3, quality3 = define_variables(m3)\n",
"define_objective(m3, flow3, quality3)\n",
"define_constraints(m3, flow3, quality3)\n",
"\n",
"m3.params.NonConvex=2\n",
"m3.optimize()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true
},
"outputs": [],
"source": [
"# Draw solution\n",
"draw_solution(m3, flow3, quality3)"
]
}
],
"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.6.8"
}
},
"nbformat": 4,
"nbformat_minor": 2
}