\section{Introduction}

This manuscript presents a library to formulate pooling problems using Pyomo~\citep{Hart2017,Watson2011}, a modeling library for Python, without losing network structure information.
Traditionally, mixed-integer non-linear optimization (MINLP) solvers see ``flat'' optimization problems as a series of objectives, constraints, and variables, with little to no information about what the optimization problem represents.
We propose retaining the network structure: we use Python to model the network structure with nodes representing network components, and edges representing possible network flow.

We focus on the pooling problem because of its many industrial applications, including \citep{misener2009advances}: crude-oil scheduling \citep{lee1996mixed,li2007improving,li2012continuous,li2012scheduling}, water networks \citep{galan1998optimal,castro2007efficient}, natural gas production \citep{selot2008short,li2011stochastic}, fixed-charge transportation with product blending \citep{papageorgiou2012fixed}, hybrid energy systems \citep{baliban2012global}, multi-period blend scheduling \citep{kolodziej2013discretization}, and mining \citep{boland2015special}. 
Solving the pooling problem is NP-hard \citep{alfaki2013strong,Baltean-Lugojan2018,letsios2020approximation}, so deterministic global optimization solvers use algorithms such as branch \& bound to solve the problem. Our goal with explicitly using special structure information is practically solving larger problem sizes.
But embedding structural information within \galini{} could be similarly implemented for other optimization problems over networks, e.g.\ optimal power flow.

Our pooling network library automatically generates a Pyomo block with the problem formulated as a mixed-integer quadratically constrained program (MIQCQP). The library also provides functions to compute a feasible (but not necessarily optimal) problem solution \citep{Dey2015}, and functions to generate valid linear cuts \citep{luedtke2020strong}.
We integrate the pooling network library with \galini{} \citep{ceccon2020galini, ceccon2021galini}, an extensible MIQCQP solver. The pooling network library provides an initial feasible solution search strategy and a cut generator: users enable these components by changing their \galini{} configuration file.
Both these \galini{} extensions take advantage of the pooling problem embedded inside the Pyomo model. \galini{} itself does not require any change, so everything presented in this paper is an external extension to \galini{}. The advantage of this approach is that, although \galini{}  itself is unaware of the pooling problem, \galini{} can, through this library, behave as if it has perfect knowledge of its structure.

State-of-the-art MIP solvers detect constraints that commonly appear in network problems and deduce the network structure \citep{brown1984automatic,bixby1988finding,gulpinar2004extracting,achterberg2010mcf,salvagnin2016detecting}.
Out prior work \citep{Ceccon2016} deduced MINLP network structure: we detected pooling problems in 6\% of MINLPLib. But detecting special structure in a ``flat'' optimization problem is a heuristic approach that may miss special structure developed by the user.
Our pooling network library is different because the library keeps the original pooling problem structure together with the optimization problem.
This paper is published alongside the source code \citep{ceccon2021pooling}, licensed under the Apache License, version 2.

\section{Background}

\subsection{The pooling problem}

The pooling problem has a feed-forward network $T = (N, A)$, where $N$ and $A$ are the sets of nodes and edges, respectively. Nodes $N$ are partitioned in three layers: 1) inputs $I$ representing material availability, 2) pools $L$ representing material mixing, and 3) outputs $J$ representing product demand. By definition, $N = I \cup L \cup J$. The network allows connections from inputs to pools, inputs to outputs, and pools to outputs, i.e.\ $A = (I \times L) \cup (I \times J) \cup (L \times J)$. The notation $I_l$ represents the input nodes connected to pool $l$, i.e.\ $I_l = \{i : i \in I, (i, l) \in A\}$. Similarly, $I_j$ represents the set of inputs connected to output $j$.
The optimization objective is to minimize cost (maximize profit) in the network, while tracking material qualities $K$. Tracking these material qualities K across the network introduces bilinear terms to the problem.

The first optimization formulation to solve the pooling problem was the P-formulation \citep{haverly1978studies}.
   
   
   
   
   
The mathematically-equivalent Q-formulation \citep{ben1994global} introduces fractional flow variables $y_{il} = q_{il} \sum_{j \in J} y_{lj}$.
   
 



   

   
   
The Equation~\eqref{eq:pooling.pq.formulation} PQ-formulation \citep{quesada1995global,tawarmalani2013convexification} is equivalent to the Q-formulation with the addition of Equation~\eqref{eq:pq.cut}. Tables~\ref{table:pooling.problem.notation.sets} and \ref{table:pooling.problem.notation.variables} in \ref{app:notation} summarize the notation.
\allowdisplaybreaks{
\begin{subequations}
  \label{eq:pooling.pq.formulation}
  \begin{align}
    & & & \underset{y,v,q}{\text{min}} \sum_{i \in I, l \in L} c_i v_{ilj} - \sum_{l \in L, j \in J} d_j y_{lj} - \sum_{i \in I, j \in J} (d_j - c_i) y_{ij} \\
    & &\begin{array}{r} \text{Path} \\ \text{Definition} \end{array} &\left[ v_{ilj} = q_{il} y_{lj} \quad \forall i \in I, l \in L, j \in J
                                                                       \vphantom{\sum_{i \in I}} \right. \label{eq:pooling.pq.pathdef}\\
    & & \begin{array}{r} \text{Simplex} \end{array} & \left[ \sum_{i \in I} q_{il} = 1 \quad \forall l \in L \right. \label{eq:pooling.pq.simplex}\\
   
    & & \begin{array}{r} \text{Product} \\ \text{Quality} \end{array}&
                                                                       \left[
                                                                       \begin{array}{lr}
                                                                       \displaystyle \sum_{i \in I, l \in L} C_{ik} v_{ilj} + \displaystyle \sum_{i \in I} C_{ik} y_{ij}  & \left\{
                                                                       \begin{array}{l}
                                                                         \leq P_{jk}^U \left( \displaystyle \sum_{l \in L} y_{lj} + \sum_{i \in I} y_{ij} \right)\\
                                                                         \geq P_{jk}^L \left( \displaystyle  \sum_{l \in L} y_{lj} + \sum_{i \in I} y_{ij} \right)
                                                                       \end{array}
    \right.\\
    & \quad \forall k \in K, j \in J
    \end{array}
     \right. \label{eq:pooling.pq.product.quality}\\
 
    & & \begin{array}{r} \text{Input} \\ \text{Capacity} \end{array}&\left[ A_i^L \leq \sum_{l \in L, j \in J} v_{ilj} + \sum_{j \in J} y_{ij} \leq A_i^U \quad \forall i \in I \right.\label{eq:pooling.pq.input.capacity}\\
   
    & & \begin{array}{r} \text{Pool} \\ \text{Capacity} \end{array}&\left[ \sum_{l \in L} y_{lj} \leq S_l \quad \forall l \in L\right. \label{eq:pooling.pq.pool.capacity}\\
   
    & & \begin{array}{r} \text{Output} \\ \text{Capacity} \end{array}&\left[ D_j^L \leq \sum_{l \in L} y_{lj} + \sum_{i \in I} y_{ij} \leq D_j^U \quad \forall j \in J \right. \label{eq:pooling.pq.output.capacity}\\
   
    & & \begin{array}{r} \text{Reduction} \end{array}&\left[\sum_{i \in I} v_{ilj} = y_{lj} \quad \forall l \in L, j \in J \right. \label{eq:pooling.pq.reduction.1}\\
    & & \begin{array}{r} \text{Reduction} \end{array}&\left[\sum_{i \in I, l \in L} v_{ilj} \leq  c_l q_{il} \quad \forall i \in I, l \in L \right. \label{eq:pooling.pq.reduction.2}\\
   
    & & \begin{array}{r} \text{Other} \\ \text{Bounds} \end{array}&\left[
                                                                    \begin{array}{l}
                                                                      \sum_{j \in J} v_{ilj} \leq c_{il} \quad \forall i \in I, l \in L\\
                                                                      y_{lj} \leq c_{lj} \quad \forall l \in L, j \in J\\
                                                                      y_{ij} \leq c_{ij} \quad \forall i \in I, j \in J\\
                                                                      y_{ij}, y_{lj}, v_{ilj}, q_{il} \geq 0 \quad \forall i \in I, l \in L, j \in J\\
                                                                      \sum_{j \in J} v_{ilj} \leq c_l q_{il} \quad \forall i \in I, l \in L
                                                                    \end{array}
    \right.
    \label{eq:pooling.pq.bounds}
  \end{align}
\end{subequations}
}

\begin{equation}
  \label{eq:pq.cut}
  \sum_{i \in I} q_{il} y_{lj} = y_{lj} \quad \forall l \in L, j \in J
\end{equation}

Other pooling formulations incorporate other algorithmic trade-offs \citep{audet2004pooling,alfaki2013multi,boland2016new}.

\subsection{Pyomo}

Pyomo \citep{Hart2017, Friedman2013} is a Python-based algebraic modeling language to formulate MINLP. Other algebraic modeling languages include GAMS \citep{brooke1996gams}, AMPL \citep{fourer2003ampl}, AIMMS, YALMIP \citep{lofberg2004yalmip}, and JuMP \citep{dunning2017jump}.

Large scale optimization models describe complex physical systems comprised of many sub-systems.
Pyomo provides a \texttt{Block} component to structure models. Blocks encapsulate individual model components in reusable pieces that are combined together to form the final model. For example, users can develop blocks representing their system at each time period and combine multiple instances of the block to produce a scheduling problem.

Pyomo blocks can be nested into other blocks, creating a tree-like structure of blocks. Figure~\ref{fig:pyomo.block.pooling} shows how to organize a large model to optimize refinery production subject to US Environmental Protection Agency (EPA) regulation \citep{Misener2010d}. The root \texttt{ConcreteModel} has two blocks, one containing the refinery production model, the other containing EPA constraints. The EPA constraints contain complex non-linear expressions and logical disjunctions, so we organize them further in separate blocks. We add one block for each emission model to make the model easier to update and debug. Blocks may use variables defined in other blocks, both from their children blocks and their siblings.

\begin{figure}
  \centering
  \includegraphics[width=1.0\textwidth]{pyomo-block-pooling-problem.pdf}
  \caption{Pyomo model for the pooling problem with EPA constraints \citep{Misener2010d}. The model contains one reusable model with the pooling problem PQ-formulation and one block with the additional EPA constraints. Each block contains definitions for its variables and constraints, the EPA block uses variables from the PQ-formulation block.}
  \label{fig:pyomo.block.pooling}
\end{figure}

Developers can also create \textit{smart} blocks, i.e.\ new types of Pyomo blocks, which can create and update their internal constraints. These custom blocks can automatically generate relaxations of a particular expression or build an optimization model from an higher-level definition. \ref{app:pyomo_smart_block} shows how to define a new Pyomo block that relaxes a nonconvex bilinear term.


\subsection{\galini{}}

\galini{} is an open-source MIQCQP solver written in Python that can be extended at runtime.
\galini{} implements a generic branch \& cut algorithm: at each node of the branch \& bound algorithm \galini{} iteratively adds cuts and solves the linear relaxation of the original optimization problem.
Algorithm developers extend \galini{} by registering their classes using Python \textit{entry points}, a system-wide registry of Python classes.
This manuscript extends \galini{} with a new cut generator and an initial primal search heuristic. \galini{} uses cut generators inside the cut loop to augment the linear relaxation of the optimization problem. \galini{} calls the initial primal search heuristic before visiting the branch \& bound root node to find an initial feasible solution.

The cut generator interface requires to implement six methods:

\begin{itemize}
\item \texttt{before\_start\_at\_root(problem, relaxed\_problem)}: called before entering the cut loop at the root node,
\item \texttt{after\_end\_at\_root(problem, relaxed\_problem, solution)}: called after the cut loop at the root node,
\item \texttt{before\_start\_at\_node(problem, relaxed\_problem)}: called before entering the cut loop  at non-root nodes,
\item \texttt{after\_end\_at\_node(problem, relaxed\_problem, solution)}: called after the cut loop at non-root nodes,
\item \texttt{has\_converged(state)}: returns a true value if the cut generator won't generate any more cuts at this node,
\item \texttt{generate(problem, relaxed\_problem, solution, tree, node)}:
  return a list of cuts to be added to \texttt{relaxed\_problem}, can use information from the previous iteration \texttt{solution}.
\end{itemize}

The \texttt{InitialPrimalSearchStrategy} requires algorithm developers to implement the following method:

\begin{enumerate}
\item \texttt{\textbf{solve(model, tree, node)}}: find a feasible solution for the model if any, or \texttt{None} if not possible.
\end{enumerate}

\section{Pooling network library}
\label{chapt:network.structure}

\subsection{Network structure}

The library's \texttt{Network} class describes the network structure of the problem. This manuscript assumes that a user explicitly stores the network structure. But we could alternatively detect pooling network structure \citep{Ceccon2016} and then use pooling network library to develop algorithmic approaches, e.g. cutting planes and primal heuristic.
\ref{app:network_structure} shows how the pooling library builds and stores a problem. As an example, Figure~\ref{fig:pooling.problem.network.library} shows how the library represents the \texttt{adhya4} test problem \citep{adhya1999lagrangian}. 
\begin{figure}%
  \centering
  \includegraphics[width=0.90\textwidth]{network-library.pdf}
  \caption{Graphical representation of the pooling problem \texttt{adhya4} as implemented with the pooling network library. Nodes are grouped by layer, in this example there are 3 layers (one for inputs, one for pools, and one for outputs). Nodes have the attributes: \texttt{name}, \texttt{layer}, \texttt{capacity}, and \texttt{cost}, and the (optional) \texttt{attr} dictionary to add application-specific attributes. Edges have the attributes: \texttt{source}, \texttt{destination}, \texttt{cost} (per unit flow cost), \texttt{fixed\_cost}, \texttt{capacity}, and the (optional) \texttt{attr} dictionary to add application-specific attributes The figure shows the predecessors and successors of node \texttt{o2} highlighted in purple and green, respectively.}%
  \label{fig:pooling.problem.network.library}%
\end{figure}


\subsection{Automatic PQ-formulation generation}

The library provides a custom \texttt{PoolingPQFormulationBlock} Pyomo block to generate the pooling problem PQ-formulation in Optimization Problem~\eqref{eq:pooling.pq.formulation} from the network structure. The library does not check that the network passed to it is a valid pooling problem, but it expects a network with a pooling problem-like structure. If the network contains edges that are not part of a standard pooling problem, the library ignores these edges. The network is expected to have a three layer structure. If the network has more than three layers of nodes, then all the layers after the third are ignored.
The pooling PQ-formulation block can be used in an indexed Pyomo block to have multiple pooling problems in the same optimization problem. The block adds the variables and constraints to the block, so, after changing anything in the optimization problem, the PQ-formulation needs to be rebuilt by calling the \texttt{rebuild()} method on the block. Table~\ref{table:block.formulation} summarizes the block's variables and how they relate to the PQ-formulation. 
Figure~\ref{fig:pooling.problem.network.block} shows the structure of a \texttt{PoolingPQFormulationBlock} Pyomo block, with the pooling problem constraints stored together with the network structure.

\begin{table}
  \begin{center}
    \begin{tabularx}{\textwidth}{@{} l l X @{}}
      \toprule
      Type & Name & Description \\
      \midrule
      Variables & \texttt{q[i, l]} & Fractional flow: input $i$ to pool $l$, as fraction of total flow through $l$ \\
      & \texttt{v[i, l, j]} & Flow: input $i$ to output $j$ through pool $l$ \\
      & \texttt{y[l, j]} & Pool $l$ to output $j$ flow \\
      & \texttt{z[i, j]} & Input $i$ to output $j$ flow \\
      \midrule
      Constraints & \texttt{path\_definition[i, l, j]} & Equation~\eqref{eq:pooling.pq.pathdef}\\
      & \texttt{simplex[l]} & Equation~\eqref{eq:pooling.pq.simplex}\\
      & \texttt{product\_quality\_lower\_bound[j, k]} & Equation~\eqref{eq:pooling.pq.product.quality} \\
      & \texttt{product\_quality\_upper\_bound[j, k]} & Equation~\eqref{eq:pooling.pq.product.quality} \\
      & \texttt{input\_capacity[i]} & Equation~\eqref{eq:pooling.pq.input.capacity} \\
      & \texttt{pool\_capacity[l]} & Equation~\eqref{eq:pooling.pq.pool.capacity} \\
      & \texttt{output\_capacity[j]} & Equation~\eqref{eq:pooling.pq.output.capacity} \\
      & \texttt{reduction\_1[l, j]} & Equation~\eqref{eq:pooling.pq.reduction.1}\\
      & \texttt{reduction\_2[i, l]} & Equation~\eqref{eq:pooling.pq.reduction.2}\\
      \bottomrule
    \end{tabularx}
  \end{center}
  \caption{Top-level definitions of \texttt{PoolingPQFormulationBlock}. Users can use the block variables in constraints outside of the block. Users can also \texttt{deactivate()} the constraints.}
  \label{table:block.formulation}
\end{table}

The library also provides convenience functions to iterate over the indexes used in a pooling problem. For example, the function \texttt{index\_set\_ilj} will return an iterator over the indexes \texttt{(i, l, j)} of variables \texttt{v[i, l, j]}, equivalent to the variables $v_{ilj}$ in the pooling problem PQ-formulation.

\begin{figure}%
  \centering
  \includegraphics[width=0.40\textwidth]{pq-formulation-block.pdf}
  \caption{\texttt{PoolingPQFormulation} Pyomo block. The block automatically generates the variables and constraints for the pooling PQ-formulation. The block keeps a reference to the pooling problem network structure. Software that is aware of this structure can use it to have a precise knowledge of the pooling problem.}%
  \label{fig:pooling.problem.network.block}%
\end{figure}

\subsection{Convex relaxation and cut generator}

The pooling network library implements a convex relaxation of the pooling problem \citep{luedtke2020strong}.
The relaxation introduces a variable $z_{lj}$ that represents the total flow to output $j$ that does not pass through pool $l$:
\begin{equation*}
  \label{eq:pooling.cut.z.def}
  z_{lj} := \sum_{i \in I} y_{ij} + \sum_{l' \in L, l' \neq l} y_{l'j} \quad \forall l \in L, j \in J
\end{equation*}
\noindent
The contribution of flow through pool $l$ to the quality $k$ at output $j$ is defined:
\begin{equation*}
  \label{eq:poolin.cut.u.def}
  u_{ljk} := \sum_{i \in I} (P_{jk}^U - C_{ik})v_{ilj} \quad \forall l \in L, j \in J, k \in K
\end{equation*}
\noindent
The total flow not going through pool $l$ to the quality $k$ at output $j$ is:

\begin{equation*}
  \label{eq:pooling.cut.t.def}
  t_{ljk} := \sum_{i \in I} (P_{jk}^U - C_{ik}) y_{il} + \sum_{i \in I} \sum_{l' \in L, l' \neq l} (P_{jk}^U - C_{ik}) v_{il'j}
  \quad \forall l \in L, j \in J, k \in K
\end{equation*}
\noindent
The quality of attribute $k$ at pool $l$ is:
\begin{equation*}
  \label{eq:pooling.cut.p.def}
  p_{ljk} := \sum_{i \in I} (P_{jk}^U - C_{ik}) q_{il} \quad \forall l \in L, j \in J, k \in K
\end{equation*}
\noindent
Scaled flow variable $s_{lj}$ scales the flow from pool $l$ to output $j$ to be in $[0, 1]$:
\begin{equation*}
  \label{eq:pooling.cut.s.def}
  s_{lj} := \frac{1}{c_j}\sum_{i \in I} v_{ilj} \quad \forall l \in L, j \in J
\end{equation*}
\noindent
We add the McCormick envelope of $r_{ljk} = s_{lj} p_{ljk}$, where $p_{ljk}^L$ and $p_{ljk}^U$ are the bounds on $p_{ljk}$:
\begin{align}
  \label{eq:pooling.cut.uxp}
  r_{ljk} - p_{ljk}^L s_{lj} & \geq 0 \quad \forall l \in L, j \in J\\
  p_{ljk}^U s_{lj} - r_{ljk} & \geq 0 \quad \forall l \in L, j \in J\\
  r_{ljk} - p_{ljk}^L x & \leq p_{ljk} - p_{ljk}^L \quad \forall l \in L, j \in J\\
  p_{ljk}^U s_{lj}  - r_{ljk} & \leq p_{ljk}^U - p_{ljk} \quad \forall l \in L, j \in J
\end{align}
Parameters $\underline{\eta}_{ljk}$ and $\overline{\eta}_{ljk}$ track the lower and upper bounds on the excess of quality $k$ at output $j$ over the inputs $i$ that are connected to pool $l$.
\begin{align}
  \underline{\eta}_{ljk} & := \min\{(P_{jk}^U - C_{ik}) : i \in I_l\}  \quad \forall l \in L, j \in J\\
  \overline{\eta}_{ljk} & := \max\{(P_{jk}^U - C_{ik}) : i \in I_l\} \quad \forall l \in L, j \in J
\end{align}
\noindent
Similarly, the parameters $\underline{\beta}_{ljk}$ and $\overline{\beta}_{ljk}$ track the bounds on the excess of quality $k$ at output $j$ over the inputs $i$ that are \textit{not} connected to pool $l$.
\begin{align}
  \underline{\beta}_{ljk} & := \min\{(P_{jk}^U - C_{ik}) : i \in I_j \cup \bigcup_{l' \in L \setminus \{l\}} I_{l'}\} \quad \forall l \in L, j \in J\\
  \overline{\beta}_{ljk} & := \max\{(P_{jk}^U - C_{ik}) : i \in I_j \cup \bigcup_{l' \in L \setminus \{l\}} I_{l'}\} \quad \forall l \in L, j \in J
\end{align}

\noindent
Figure~\ref{fig:pooling.problem.network.library.cuts} shows the variables and sets that, for each pool, output, quality triplet, generate the convex relaxation described in this section.

\begin{figure}%
  \centering
  \includegraphics[width=0.6\textwidth]{pooling-problem-cuts.pdf}
  \caption{Variables and sets generating cuts for each pool, output, quality triplet. This example considers pool \texttt{o1} and output \texttt{p2}. Variable \texttt{s[l, j]} is the scaled flow from pool $l$ to output $j$. Variable \texttt{z[l, j]} is the flow to output $j$ not passing through pool $l$, i.e.\ the flow to output \texttt{p2} passing through pool \texttt{o2}. The inputs circled in blue are the inputs connected to the output $j$ through pool $l$ that compute $\underline{\eta}_{ljk}$ and $\overline{\eta}_{ljk}$. The inputs circled in orange are the inputs connected to output $j$ through the by-pass flow, and compute $\underline{\beta}_{ljk}$ and $\overline{\beta}_{ljk}$.}%
  \label{fig:pooling.problem.network.library.cuts}%
\end{figure}



Linear inequality \eqref{eq:pooling.cut.eq22} is valid for all $(l, j, k) \in L \times J \times K$ where $\overline{\beta}_{ljk} > 0$:

\begin{equation}
  \label{eq:pooling.cut.eq22}
  (\overline{\eta}_{ljk} - \underline{\eta}_{ljk})t_{ljk} + \underline{\eta}_{ljk}(\overline{\eta}_{ljk} s_{lj} - u_{ljk}) + \overline{\beta}_{ljk}(u_{ljk} - \underline{\eta}_{ljk} s_{lj}) \leq \overline{\beta}_{ljk}(p_{ljk}  - \underline{\eta}_{ljk})
\end{equation}

\noindent
Linear inequality \eqref{eq:pooling.cut.eq28} is valid for all $(l, j, k) \in L \times J \times K$ where $\underline{\beta}_{ljk} < 0$:

\begin{equation}
  \label{eq:pooling.cut.eq28}
  (\underline{\eta}_{ljk} - \underline{\beta}_{ljk})(\overline{\eta}_{ljk} s_{lj} - u_{ljk}) \leq -\underline{\beta}_{ljk} (\overline{\eta}_{ljk} - p_{ljk})
\end{equation}

\noindent
Convex non-linear inequality \eqref{eq:pooling.cut.eq15} is valid  if $\underline{\beta}_{ljk} < 0$:

\begin{equation}
  \label{eq:pooling.cut.eq15}
  (u_{ljk} - \underline{\beta}_{ljk} s_{lj})(u_{ljk} - \underline{\eta}_{ljk} s_{lj}) \leq -\underline{\beta}_{ljk} s_{lj} (p_{ljk} - \underline{\eta}_{ljk})
\end{equation}

\noindent
Convex non-linear inequality \eqref{eq:pooling.cut.eq18} is valid if $\overline{\beta}_{ljk} > 0$ and $\underline{\eta}_{ljk} < 0$ when $t_{ljk} > 0$:

\begin{equation}
  \label{eq:pooling.cut.eq18}
  (\overline{\eta}_{ljk} - \underline{\eta}_{ljk}) t_{ljk} + \overline{\beta}_{ljk} (\overline{\eta}_{ljk} s_{lj} - u_{ljk}) +
  \frac{\underline{\eta}_{ljk} t_{ljk} (u_{ljk} - \underline{\eta}_{ljk} s_{lj})}{t_{ljk} + u_{ljk} - \underline{\eta}_{ljk} s_{lj}} \leq
  \overline{\beta}_{ljk} (\overline{\eta}_{ljk} - p_{ljk})
\end{equation}

Library function \texttt{add\_all\_pooling\_inequalities(block, parent, network)}, which should only be called once, adds the pooling inequalities variables and equations to \texttt{block}, using flow variables from the \texttt{PoolingPQFormulation} block \texttt{parent} and the pooling problem \texttt{network}. This function also adds the Equation~\eqref{eq:pooling.cut.eq22} \& \eqref{eq:pooling.cut.eq28} inequalities. The function adds the following variables to the Pyomo block:

\begin{enumerate}
\item \texttt{z[l, j]}: the total flow $z_{lj}$ to output $j$ that does not go through pool $l$,
\item \texttt{s[l, j]}: scaled flow $s_{lj}$ from pool $l$ to output $j$,
\item \texttt{u[j, k, l]}: the contribution $u_{ljk}$ of flow through pool $l$ to the quality $k$ at output $j$,
\item \texttt{y[j, k, l]}: the contribution $t_{ljk}$ of the total flow not going through pool $l$ to the quality $k$ at output $j$,
\item \texttt{t[j, k, l]}: the quality $p_{ljk}$ of attribute $k$ going through pool $l$.
\end{enumerate}

The function uses Coramin to add the McCormick envelope of $r_{ljk} = s_{lj} p_{ljk}$ described in Equation~\eqref{eq:pooling.cut.uxp} to the Pyomo block. Finally, the function creates an \texttt{inequalities} \texttt{ConstraintList} that contains all the inequalities from Equation~\eqref{eq:pooling.cut.eq22} \& \eqref{eq:pooling.cut.eq28}. The function iterates over all the combinations of pool, output, quality in the network and checks the conditions of $\underline{\beta}_{ljk}$ and $\overline{\beta}_{ljk}$ to decide whether or not to add the inequality to the list.
This function adds all possible inequalities. Alternatively, we could add the inequalities as cutting planes, e.g.\ if a relaxation violates the inequalities by more than a parameter $\varepsilon$. Our computational tests on the large scale instances of \citet{Dey2015} and \citet{luedtke2020strong} indicate that the change in performance is not significant.

Equations~\eqref{eq:pooling.cut.eq15} \& \eqref{eq:pooling.cut.eq18} are quadratic and cannot be added directly to the linear relaxation of the pooling problem. The relaxation adds a gradient inequality based on Equation~\eqref{eq:pooling.cut.eq15} for every $(l, j, k) \in L \times J \times K$ where Equation~\eqref{eq:pooling.cut.eq15} is violated more than a parameter $\varepsilon$.
The relaxation adds a gradient inequality based on Equation~\eqref{eq:pooling.cut.eq18} for every $(l, j, k) \in L \times J \times K$ where Equation~\eqref{eq:pooling.cut.eq18} is violated more than a parameter $\varepsilon$ and $t_{ljk} > 0$.

The library function \texttt{add\_valid\_cuts(block, parent, network)} 
checks if any of the non-linear equations, i.e.\ Equations~\eqref{eq:pooling.cut.eq15} \& \eqref{eq:pooling.cut.eq18}, is violated and, if violated, generates a gradient inequality. This function can be called multiple times, e.g.\ in each iteration of a cut loop. 
This function iterates over all the pool, output, quality triplets and adds the cut to the \texttt{cuts} \texttt{ConstraintList} on the Pyomo block if any of the expressions in Equations~\eqref{eq:pooling.cut.eq15} \& \eqref{eq:pooling.cut.eq18} are violated by a parameter $\varepsilon$. Users wishing to generate cuts but not wanting to add the cuts to the model can use the function \texttt{generate\_valid\_cuts(block, parent, network)} that returns an iterator of cuts to be added.

The \texttt{PoolingQPFormulation} block also contains two convenience methods to make the pooling network library more user-friendly:

\begin{enumerate}
\item \texttt{add\_inequalities()}: calls the \texttt{add\_all\_pooling\_inequalities} function with the correct arguments to add the relaxation variables for each pool, output, quality triplet. This method also adds Equation~\eqref{eq:pooling.cut.eq22} \& \eqref{eq:pooling.cut.eq28} if $\overline{\beta}_{ljk} > 0$ and $\underline{\beta}_{ljk} < 0$ respectively,
\item \texttt{add\_cuts()}: calls the \texttt{add\_valid\_cuts} function to add gradient inequalities based on Equation~\eqref{eq:pooling.cut.eq15} \& \eqref{eq:pooling.cut.eq18}.
\end{enumerate}

Listing~\ref{lst:pooling.network.library.add.cuts} in \ref{app:code} uses these methods to add valid inequalities to the linear relaxation of the pooling problem. This example also shows how to use the \texttt{add\_cuts} method to generate cuts at every iteration of a simple cut loop.


\subsection{Mixed-integer programming restriction}
\label{sect:network.primal.heuristic}

The library implements a MIP restriction of the pooling problem \citep{Dey2015}.
The heuristic ``splits'' each problem pool $l$ in $\tau$ copies, each copy receives a predefined fraction of the total flow received by pool $l$. Each copy of the pool can send flow only to one output node. Equation~\eqref{eq:pooling.mip.restriction} describes the MIP restriction, which introduces auxiliary flow variables $w_{il_tj}$ representing the flow from input $i$ to output $j$ through pool $l$ and auxiliary pool $t$. Each auxiliary pool $t$ receives a fraction $\gamma_{lt}$ of the total flow into pool $l$. The restriction also introduces binary variables $\zeta_{l_tj}$ that equal $1$ if the auxiliary pool $t$ of pool $l$ sends flow to output $j$, since the heuristic requires that each auxiliary pool sends flow to at most one output, we require that $\sum_{j \in J} \zeta_{l_tj} = 1, \forall l \in L, t \in \{1,\dots,\tau\}$.

\begin{subequations}
  \label{eq:pooling.mip.restriction}
  \begin{align}
    & & & \underset{y,v,q}{\text{max}} \sum_{i \in I, j \in J} f_{ij} y_{ij} + \sum_{i \in I, l \in L, j \in J} (f_{il} + f_{lj}) v_{ilj}\\
   
    & & \begin{array}{r} \text{Flow} \\ \text{Balance} \end{array} &\left[ v_{ilj} = \sum_{t = 1}^{\tau} w_{il_tj} \quad \forall i \in I, l \in L, j \in J
                                                                     \vphantom{\sum_{i \in I}}  \right. \label{eq:pooling.mip.restriction.flow.balance}\\
   
    & & \begin{array}{r} \text{Fractional Flow} \\ \text{Definition} \end{array} &\left[ \sum_{j \in J} w_{il_tj} = \gamma_{lt} \sum_{j \in J} v_{ilj} \quad \forall i \in I, l \in L, t \in \{1,\dots,\tau\} \right. \label{eq:pooling.mip.restriction.fract.flow.def}\\
   
    & & \begin{array}{r} \text{Bounds} \end{array}&\left[
                                                                    \begin{array}{l}
                                                                      0 \leq w_{il_tj} \leq c_{lj} \zeta_{l_tj} \quad \forall i \in I, l \in L, t \in \{1,\dots,\tau\}, j \in J\\
                                                                      \sum_{j \in J} \zeta_{l_tj} = 1 \quad \forall l \in L, t \in \{1,\dots,\tau\}\\
                                                                      \zeta_{l_tj} \in \{0, 1\} \quad \forall l \in L, t \in \{1,\dots,\tau\}, j \in J
                                                                    \end{array}
    \right.\label{eq:pooling.mip.restriction.extra}\\
    & & \begin{array}{r} \text{PQ-formulation} \end{array}&\left[\text{Equation~\eqref{eq:pooling.pq.simplex} -- Equation~\eqref{eq:pooling.pq.bounds}} \right. \notag
  \end{align}
\end{subequations}

Figure~\ref{fig:pooling.problem.heuristic} diagrams the pooling problem and its MIP restriction.
The heuristic is implemented with two functions: one adds auxiliary variables and constraints to the \texttt{PoolingPQFormulation} block and a second restores the block to its original state. As a convenience to users, the library also implements a Python \textit{context manager} to automatically add and remove the MIP restriction. Listing~\ref{lst:network.add.heuristic} in \ref{app:code} shows how users can build the MIP heuristic.

The context manager \texttt{mip\_heuristic(block, network, tau, weights)} expects the following parameters:

\begin{itemize}
\item \texttt{block}: the \texttt{PoolingPQFormulation} block,
\item \texttt{network}: the pooling problem network,
\item \texttt{tau}: the number $\tau$ of auxiliary pools each pool is ``split'' into,
\item \texttt{weights}: a function to generate the weights $\gamma_{lt}$ for each pool $l$ and auxiliary pool $t$. By default, the library uses uniform weights $1/\tau$.
\end{itemize}

The \texttt{mip\_heuristic} context manager modifies the input block by deactivating the block's \texttt{path\_definition} constraints and adding a \texttt{\_pooling\_mip\_heuristic} sub-block with the following variables:

\begin{itemize}
\item \texttt{w[i,l,t,j]}: auxiliary flow variables $w_{il_tj}$,
\item \texttt{zeta[l,t,j]}: binary variables $\zeta_{l_tj}$.
\end{itemize}

The sub-block contains the following constraints:

\begin{itemize}
\item \texttt{flow\_balance[i,l,j]}: contains Equation~\eqref{eq:pooling.mip.restriction.flow.balance},
\item \texttt{flow\_balance\_2[i,l,t]}: contains Equation~\eqref{eq:pooling.mip.restriction.fract.flow.def},
\item \texttt{flow\_choice\_limit[i,l,t,j]}: contains the Equation~\eqref{eq:pooling.mip.restriction.extra} $w_{il_tj}$ bounds,
\item \texttt{flow\_choice[l,t]}: limits each auxiliary pool \texttt{t} to send flow to 1 output.
\end{itemize}

Solving the MIP restriction does not produce variable values $q_{il}$, so the library function \texttt{derive\_fractional\_flow\_variables(block)} derives the $q_{il}$ values from the other network flows, i.e.\ $q_{il} = \frac{v_{ilj}}{y_{lj}} = \frac{v_{ilj}}{\sum_{i \in I}v_{ilj}}$.

\begin{figure}%
  \centering
  \subfloat[\centering Example pooling problem with 3 inputs, 2 pools, and 2 outputs.]{{\includegraphics[width=0.5\textwidth]{pooling-problem.pdf} }}%
 
  \subfloat[\centering MIP restriction of the pooling problem in (a).]{{\includegraphics[width=0.5\textwidth]{pooling-problem-heuristic.pdf}}}%
  \caption{Pooling problem network structure of \texttt{adhya4} (a) together with its MIP restriction (b). The MIP restriction ``splits'' each pool into $\tau$ ($\tau = 2$ in this example) auxiliary pools, each auxiliary pool receives a predetermined fraction of the total flow into its pool and can send flow to only one output. There is no cost associated to the flow between a pool and its auxiliary pools. The blue arrows represent edges with a positive flow in the problem optimal solution and the MIP restriction solution, we can see that in this example the two optimal network configurations are different.}%
  \label{fig:pooling.problem.heuristic}%
\end{figure}

\section{\galini{} integration}
\subsection{Initial primal heuristic}
\label{galini.primal.heuristic}

The pooling network library integrates the Section~\ref{sect:network.primal.heuristic} MIP heuristic by providing a \texttt{PoolingPrimalSearchStrategy} class implementing the \galini{} \texttt{InitialPrimalSearchStrategy} interface and registering it with \galini{} to be available at runtime. \galini{} uses this interface to find an initial feasible solution at the root node only.
The primal search strategy detects the Pyomo block that contains the pooling problem and builds the MIP restriction using $\tau = 1$ 
\citep{Dey2015}. 
The search strategy then solves the MIP restriction using the same MIP solver used by the \galini{} branch \& cut algorithm. To limit the time spent finding the feasible solution, the search strategy sets the MIP solver time limit option to $60$ seconds and the optimality gap to $1\%$.

If the MIP heuristic is successful, the search strategy derives the values for the fractional flow variables $q_{il}$ and then it returns the feasible solution to the original user model to \galini{}.

The primal heuristic is manually activated by setting \galini{} configuration \texttt{intial\_search\_strategy} to \texttt{'pooling'}. \texttt{PoolingPrimalSearchStrategy} does not perform additional checks on the model's additional constraints. If the additional constraints are linear, the primal search strategy will succeed. If the model contains additional non-linear constraints, the primal search strategy will fail since the MIP solver is not able to handle them and Pyomo will throw an exception. We decided not to disable the non-linear constraints since by doing that the new model would not be a restriction of the original model, creating a subtle bug that is difficult to debug for users.


\subsection{Cuts generator}

The pooling network library implements a \galini{} \texttt{CutsGenerator} to generate valid cuts for the pooling problem at each node of the branch \& bound algorithm. The cut generator implements the following interface methods:

\begin{itemize}
\item \texttt{before\_start\_at\_root(problem, linear\_problem)}: this callback is called before the cut loop at the root node. It checks whether \texttt{linear\_problem} contains a \texttt{PoolingPQFormulation} block, if not sets a variable to disable the cut generator. This callback also adds valid inequalities to \texttt{linear\_problem}.
\item \texttt{generate(problem, linear\_problem, lp\_solution, tree, node)}: this method is called at each iteration of the cut loop to generate new cuts to be added to \texttt{linear\_problem}. The function calls the \texttt{generate\_valid\_cuts} function from the pooling network library and then returns a list of cuts that \galini{} adds to \texttt{linear\_problem}.
\end{itemize}

Users can activate this cut generator by adding \texttt{'pooling'} to the list of cuts generator in the \galini{} configuration file.

The cuts generator initializes the linear relaxation of the pooling problem at the root node.
The pooling flow variable bounds change throughout the branch \& bound algorithm, but our preliminary results suggest that these
changes do not significantly affect the \citet{luedtke2020strong} relaxation. Therefore, we do not recompute the relaxation equations in the branch \& bound tree.

\section{Computational Results}

All tests are performed on a machine with an Intel Core i7-6700HQ CPU and 32Gb of RAM. We use \texttt{pooling\_network} version $1.0.0$ \citep{ceccon2021pooling}, \galini{} version $1.0.1$ \citep{ceccon2021galini} using CPLEX version $12.10$ as the MIP solver and Ipopt \citep{wachter2006} version $3.13.2$ to locally solve NLPs. For comparison, we use GUROBI version $9.1$ together with the Pyomo \texttt{gurobi\_direct} solver. We set the optimality relative gap to $10^{-6}$, the optimality absolute gap to $10^{-8}$, and the time limit to $600$ seconds, all other options are left as default. Additionally, we set Gurobi \texttt{NonConvex} option to $2$. We set GUROBI absolute gap by changing the \texttt{MIPGapAbs} option, and we use the \texttt{MIPGap} option to change the relative gap.

The absolute gap is the \textit{absolute difference} between the best possible objective value and the objective value of the best feasible solution.
The relative gap is the \textit{relative difference} $d(a, b)$ between the best possible objective value $a$ and the objective value of the best feasible solution $b$:

\begin{equation}
  d(a, b) = \begin{cases}
    |a - b|/\max\{|a|, |b|\} & \mbox{if } a \neq 0, b \neq 0\\
    \infty & \mbox{if } |a| = \infty \mbox{ or } |b| = \infty\\
    |a - b|/\varepsilon & \mbox{otherwise}, \\
  \end{cases}
  \label{eq:relative_difference}
\end{equation}

\noindent
where $\varepsilon$ is a small quantity. In the instances studied in this section, the objective values are always non-zero and the choice of $\varepsilon$ is irrelevant.

The pooling network library has functions, available under the \texttt{pooling\_} \texttt{network.intances} namespace, to automatically build the network structure of standard pooling problems used in these computational results. The data set summarized in Table~\ref{tbl:pooling.instances} contains classic small-scale problems \citep{tawarmalani2013convexification}, sparse instances generated based on the Haverly instances \citep{luedtke2020strong}, and large-scale dense instances \citep{Dey2015}.


We test Gurobi with the following configurations:

\begin{itemize}
\item \textbf{Default}: default Gurobi configuration.
\item \textbf{Disable Cuts}: disable cuts generation by setting the \texttt{Cuts} option to $0$.
\item \textbf{Pooling Cuts}: generate \citet{luedtke2020strong} cuts. The test linearizes the optimization problem with the convex hulls of the bilinear terms $v_{ilj} = q_{il} y_{lj}$. This configuration solves the linear relaxation of the pooling problem to generate \citet{luedtke2020strong} cuts with at most $20$ cut loop iterations. Finally, the test reintroduces the bilinear terms and solves the problem again with Gurobi. This test considers the time spent solving the linear relaxations and generating cuts as part of the Gurobi solving time.
\item \textbf{Pooling Cuts Oracle}: this configuration generates the same cuts as the \textbf{Pooling Cuts} test, but the time spent generating the cuts and solving the linear relaxations is excluded. This test assumes Gurobi had access to an oracle providing a list of valid cuts.
\item \textbf{Warm Start}: this configuration uses the MIP restriction of the pooling problem to pass a feasible solution to Gurobi. This test considers the time spent solving the MIP restriction as part of the Gurobi solving time.
\item \textbf{Warm Start Oracle}: this configuration uses the same MIP restriction as the previous test, but the time spent solving the MIP restriction is excluded. This test is equivalent to an oracle that provides a good feasible solution to the pooling problem.
\end{itemize}

We test \galini{} with the following configurations:

\begin{itemize}
\item \textbf{Default}: default \galini{} configuration, without pooling specific cuts generator and heuristic active.
\item \textbf{Pooling Cuts Generator}: enable the pooling cuts generator, use the default initial primal heuristic.
\item \textbf{Primal Heuristic}: enable the pooling initial primal heuristic, use the default cuts generators.
\item \textbf{Pooling Cuts Generator and Primal Heuristic}: enable both the pooling cuts generator and the initial primal heuristic.
\end{itemize}

To understand the default configuration, Table~\ref{tbl:large.scale.gurobi.cuts} presents the number of each cut type generated by Gurobi, together with the average number of branch \& bound nodes visited. 
RLT cuts are the only class of cutting planes that is generated consistently for all pooling problems.

We compare configurations with performance profiles \citep{dolan2002benchmarking}, which measure performance of a set of solver configurations $S$ on a test set  $P$. For each problem $p \in P$ and solver $s \in S$, we define a performance measure $t_{ps}$. This performance measure can be the run time needed to solve the problem to global optimality, or the relative gap after a fixed amount of time. For each problem and solver, we define the performance ratio:

\begin{equation}
  \label{eq:large.scale.pp.perf.ratio}
  r_{ps} = \frac{t_{ps}}{\min\{t_{ps} : s \in S\}}.
\end{equation}

The performance profile $\rho_s(\tau)$ is the probability for solver $s \in S$ that a performance ratio $r_{ps}$ is within a factor $\tau \in \mathrm{R}$ of the best possible ratio.

\begin{equation}
  \label{eq:large.scale.pp.perf.profile}
  \rho_s(\tau) = \frac{1}{\text{card}(P)}\text{card}\{p \in P : r_{ps} \leq \tau\}
\end{equation}

Figure~\ref{fig:large.scale.pp.gurobi.dense} and Table~\ref{tbl:dense.runtime.gurobi} compare  Gurobi configurations on the dense instances \citep{Dey2015}.
Figure~\ref{fig:large.scale.pp.gurobi.dense} shows that Gurobi does not benefit, in these instances, from its own cutting planes: the default performance of Gurobi is similar to the performance when cutting planes are deactivated.
The many \citet{luedtke2020strong} inequalities damage Gurobi's performance, even when adding the cuts as an oracle. 
But the MIP restriction greatly benefits Gurobi, especially as an oracle.
Our experiments confirm the \citet{Dey2015} results: Gurobi is very competitive in solving smaller dense instances (\texttt{randstd11} to \texttt{randstd40}), but for bigger instances a specialized heuristic is beneficial.

Figure~\ref{fig:large.scale.pp.gurobi.sparse} and Table~\ref{tbl:sparse.runtime.gurobi} show the performance of Gurobi configurations on sparse instances  \citep{luedtke2020strong}.
For these instances, the Gurobi default configuration performs best. Adding the pooling cuts before solving the non-convex problem hurts performance because of the overhead of (i) computing cuts and (ii) solving the linear relaxation in the cut loop. The configuration where the pooling cuts are generated as if with an oracle performs similarly to the default Gurobi configuration. The initial primal heuristic is not as beneficial because Gurobi is able to find a good feasible solution by default.
We note that the Gurobi defaults are very strong in solving this class of sparse problems.

To consider whether the pooling cuts help solve more difficult sparse problems, we consider the $30$ sparse instances where the Gurobi default configuration requires $>10$ seconds to solve to global optimality. We plot the performance profile on these instances in Figure~\ref{fig:large.scale.pp.gurobi.sparse.big}: for these instances the pooling cuts are beneficial and provide a performance boost over the default Gurobi.
As observed by \citet{luedtke2020strong}, the pooling cuts are not as helpful on denser instances.
In our results the benefit provided by the pooling cuts is not as strong as in~\citet{luedtke2020strong} because Gurobi is more performant than the ones used in the original paper.

\begin{figure}
  \centering
  \subfloat[\centering]{
    \resizebox{0.48\textwidth}{!}{%
      \input{pp_gurobi_dense.pgf}
    }
  }
 
  \subfloat[\centering]{
    \resizebox{0.48\textwidth}{!}{%
      \input{pp_gurobi_dense_oracle.pgf}
    }
  }
  \caption{Performance profile comparing different Gurobi configurations on 50 dense pooling problem instances \citep{Dey2015}. Performance profiles \textbf{a} include the time spent generating the cutting planes and solving the MIP restriction in the total time limit. Performance profiles \textbf{b} use an oracle to generate the cutting planes and find an initial feasible solution.}
  \label{fig:large.scale.pp.gurobi.dense}
\end{figure}

\begin{figure}
  \centering
  \subfloat[\centering]{
    \resizebox{0.48\textwidth}{!}{%
      \input{pp_gurobi_sparse.pgf}
    }
  }
 
  \subfloat[\centering]{
    \resizebox{0.48\textwidth}{!}{%
      \input{pp_gurobi_sparse_oracle.pgf}
    }
  }
  \caption{Performance profile comparing different Gurobi configurations on 360 sparse pooling problem instances \citep{luedtke2020strong}. Performance profiles \textbf{a} include the time spent generating the cutting planes and solving the MIP restriction in the total time limit. Performance profiles \textbf{b} use an oracle to generate the cutting planes and find an initial feasible solution.}
  \label{fig:large.scale.pp.gurobi.sparse}
\end{figure}

\begin{figure}
  \centering
  \subfloat[\centering]{
    \resizebox{0.48\textwidth}{!}{%
      \input{pp_gurobi_sparse_big.pgf}
    }
  }
 
  \subfloat[\centering]{
    \resizebox{0.48\textwidth}{!}{%
      \input{pp_gurobi_sparse_big_oracle.pgf}
    }
  }
  \caption{Performance profile comparing different Gurobi configurations on 30 difficult sparse pooling problem instances~\citep{luedtke2020strong}. We consider an instance difficult if the Gurobi default configuration takes longer than $10$ seconds to solve it to optimality. Performance profiles \textbf{a} include the time spent generating the cutting planes and solving the MIP restriction in the total time limit, while performance profiles \textbf{b} use an oracle to generate the cutting planes and find an initial feasible solution.}
  \label{fig:large.scale.pp.gurobi.sparse.big}
\end{figure}

Figure~\ref{fig:large.scale.pp.galini.dense} and Table~\ref{tbl:dense.runtime.galini} compare the different \galini{} configurations and the default Gurobi configuration on 50 dense pooling instances. Activating the initial primal heuristic improves \galini{} performance and makes it comparable to Gurobi. When we compare \galini{} with Gurobi only on the 20 largest dense instances, we can see that \galini{} with the initial primal heuristic enabled has a very strong performance.
Table~\ref{tbl:large.scale.result.summary} compare the relative gaps between \galini{} (with pooling cuts and initial primal heuristic enabled) and Gurobi on dense instances: the \galini{} relative gaps are stable as the problem sizes increases, while Gurobi relative gap degrades with increasing problems size.

Figure~\ref{fig:large.scale.pp.galini.sparse} and Table~\ref{tbl:sparse.runtime.galini} compare different \galini{} configurations and the default Gurobi configuration on 360 sparse pooling instances. In this case, the \galini{} performance is not competitive with Gurobi. 
The GALINI branch \& bound algorithm is not as fast as Gurobi's in visiting the hundreds or thousands of nodes required to solve these problems.

\begin{figure}
  \centering
  \subfloat[\centering]{
    \resizebox{0.48\textwidth}{!}{%
      \input{pp_galini_dense.pgf}
    }
  }
 
  \subfloat[\centering]{
    \resizebox{0.48\textwidth}{!}{%
      \input{pp_galini_dense_big.pgf}
    }
  }
  \caption{Performance profile of different \galini{} configurations compared to Gurobi on \textbf{a)} 50 large scale dense instances~\citep{Dey2015}, and \textbf{b)} the 20 largest dense instances. The default and pooling cuts configurations did not find a feasible solution and so they are not visible in the plots.}
  \label{fig:large.scale.pp.galini.dense}
\end{figure}

\begin{figure}
  \centering
  \subfloat[\centering]{
    \resizebox{0.48\textwidth}{!}{%
      \input{pp_galini_sparse.pgf}
    }
  }
 
  \subfloat[\centering]{
    \resizebox{0.48\textwidth}{!}{%
      \input{pp_galini_sparse_big.pgf}
    }
  }
  \caption{Performance profile of different \galini{} configurations compared to Gurobi on \textbf{a)} 360 large scale sparse instances~\citep{luedtke2020strong}, and \textbf{b)} 30 difficult sparse instances.  The default and pooling cuts configurations did not find a feasible solution and so they are not visible in the plots.}
  \label{fig:large.scale.pp.galini.sparse}
\end{figure}


To understand the difference in relative gaps between \galini{} and Gurobi in the dense instances, Table~\ref{tbl:large.scale.bounds.difference} shows the relative difference between the lower bounds and upper bounds computed by the two solvers.
The objective lower bounds computed by the two solvers are fairly similar, while the upper bound computed by \galini{} for the larger problems is on average $28\%$ better than the upper bound computed by Gurobi. Therefore, \galini{} performs better on large scale pooling problems because it's able to exploit the pooling problem structure to compute a feasible solution using a pooling-specific heuristic.


\section{Conclusion}

This manuscript presented the pooling network library, a Python library to describe pooling problems and generate Pyomo models using the PQ-formulation. The pooling network library provides functions to add inequalities and cuts to build a convex relaxation of the non-convex problem, which helps in providing a stronger lower bound on the problem solution. The library also includes functions to build a MIP restriction of the pooling problem and find a feasible solution using a MIP solver.

We then showed that, by integrating the cuts generator and heuristic from Section~\ref{chapt:network.structure}, we can solve large scale dense pooling problems to a smaller gap than the best available commercial solver. We compare the lower and upper bounds produced by the two solvers and see that the \galini{} performance derives from a better primal heuristic that takes advantage of the pooling problem structure.
Our experiments also show that, on sparse pooling instances, the MIP heuristic is not as effective since Gurobi already has strong performance on this class of problems. We show that, on dense pooling instances, solving the MIP restriction of the pooling problem can improve Gurobi's performance.


\section{Acknowledgements}
This work was funded by an Engineering \& Physical Sciences Research Council Research Fellowship to RM [GrantNumber EP/P016871/1].

