Using Functional Programming to Recognize Named Structure in an Optimization Problem: Application to Pooling

Branch-and-cut optimization solvers typically apply generic algorithms, e.g., cutting planes or primal heuristics, to expedite performance for many mathematical optimization problems. But solver software receives an input optimization problem as vectors of equations and constraints containing no structural information. This article proposes automatically detecting named special structure using the pattern matching features of functional programming. Speciﬁcally, we deduce the industrially-relevant nonconvex nonlinear Pooling Problem within a mixed-integer nonlinear optimization problem and show that we can uncover pooling structure in optimization problems which are not pooling problems. Previous work has shown that preprocessing heuristics can ﬁnd network structures; we show that we can additionally detect nonlinear pooling patterns. Finding named structures allows us to apply, to generic optimization problems, cutting planes or primal heuristics developed for the named structure. To demonstrate the recognition algorithm, we use the recognized structure to apply primal heuristics to a test set of standard pooling problems. V C 2016 The Authors AIChE Journal published by Wiley Periodicals, Inc. on behalf of American Institute of Chemical Engineers AIChE J , 62: 3085–3095, 2016


Introduction and Literature Review
The pooling problem is an industrially-relevant nonconvex nonlinear optimization problem minimizing cost on a feedforward network of input nodes, intermediate storage or pool nodes, and output nodes. 1 Variants of the pooling problem have applications including 2 : crude oil scheduling, [3][4][5][6] water networks, 7 natural gas production, 8,9 fixed-charge transportation with product blending, 10 hybrid energy systems, 11 and multiperiod blend scheduling. 12 The difficulty of the pooling problem arises from the nonconvex bilinear terms representing linear blending in the pools. These bilinear terms, which are necessary for tracking quality across the network, typically multiply flow rates and concentrations or flow fractions. It is possible to integrate additional complexity into the pooling problem, e.g., allowing mutable topological decisions 13,14 or nonlinear blending rules. 15 A wide variety of pooling variants with generic process networks applications can be found in MINLPLib. 16 In solving process networks optimization problems there is a common theme: it is much easier to solve large-scale instantia-tions of the standard, archetypal pooling problem than it is to solve variants including mutable topology, nonlinear blending, or temporal aspects. For example, a recent primal heuristic performs consistently well on the order of 10k variables and constraints, 17 but the approach exploits the standard pooling network structure and does not apply to pooling variants. Other recent advances specially designed for the standard pooling problem include: developing cutting planes, 18,19 discretizing the problem, 20,21 and analyzing computational complexity. 22,23 All of these approaches require knowledge of global problem structure including network topology and node identity.
Previous work in mixed-integer optimization (MIP) has found that using network structure can significantly help generate strong cutting planes. 24 Automatically identifying these embedded networks in large-scale optimization models is NPhard, 25 but there exist several polynomial time approximation algorithms to find good networks. [25][26][27] State-of-the-art MIP solver software, including MINTO 28 and CPLEX use preprocessing heuristics to automatically find these network patterns. More recent work has considered detecting more complex structures such as multi-commodity flow 29 and recovering variable meaning from flat MIP structures. 30 Due to nonlinearities, deterministic global optimization solvers do not currently infer pooling network structure and therefore will solve the pooling problem and its variants via branch-and-cut on the McCormick 31 bilinear convex hull. [32][33][34][35][36][37][38] Recent advances in solving pooling problems are therefore unavailable to solver technology; without automatically deducing global network structure, the solvers cannot incorporate primal heuristics or cutting planes based on process network structure. But, for generic process networks problems, Additional Supporting Information may be found in the online version of this article.
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Correspondence concerning this article should be addressed to R. Misener at r.misener@imperial.ac.uk.
The copyright line in this article was changed on 16 June 2016 after original online publication. even locally-relevant, single equation reformulation toward a pooling problem may expedite deterministic global optimization solvers. [39][40][41] Solver ANTIGONE 37,38 will, in preprocessing, eliminate variables, disaggregate bilinear terms, and add Reformulation-Linearization Technique equations. Although these local transformations never consider the whole optimization problem, global solvers ANTIGONE, BARON, COU-ENNE, and LINDO can solve an additional 10% of the MINLPLib2 process networks problems after applying the ANTIGONE preprocessing. 41 Other best practices for solving process networks problems with off-the-shelf software includes asserting mass balances with the Reformulation-Linearization Technique, [42][43][44][45] disaggregating bilinear terms, 46 and developing cuts from disaggregated bilinear terms. 47,48 This manuscript proposes automatically recognizing pooling structure within a mixed integer nonlinear optimization problem (MINLP). The pooling structure inside of a generic process optimization problem is a subset of the entire problem, so specialized, pooling-specific, cutting planes will also be valid bounds for the entire process networks problem. The Results section of this article shows that primal heuristic solutions to the standard pooling problem would be a good starting point for primal heuristics for the entire optimization problem. The wider goal may be stated as detecting named structures within an MINLP optimization problem; Liberti 49 posed this broader challenge to researchers after seeing this work presented at the Oberwolfach MINLP workshop.
Identifying pooling problem structure hinges on pattern matching. Patterns are defined by which variables and coefficients are expected in a constraint and by constraint bounds. The implementation is in F#, a strongly typed functional programming language targeting the .NET runtime environment. Pattern matching, one of the most distinctive F# features, has many uses, from decomposing data to control flow. The core concept is defining how data is expected to look and acting accordingly.
This manuscript begins by defining the standard pooling problem and introduces a canonical PQ-formulation. We propose the canonical formulation to avoid developing algorithms specialized to each known formulation since there are many similar formulations. 1, 43,44,[50][51][52] The recognition algorithm simultaneously converts the MINLP to canonical pooling form and recognizes network structure. The article describes the algorithm implementation; the online supplementary material justifies why the algorithm is implemented in F#. After presenting the results, the article concludes by generalizing the work to other optimization problems.

Statement of the Archetypal Standard Pooling Problem
This section defines the standard pooling problem using the Table 1 notation. 2 The first formulation, due to Haverly, 1 is the P-formulation: Product Quality X l2L p l; k y l; j 1 X i2I C i; k z i; j P L j; k ð X l2L y l; j 1 X i2I z i; j Þ 8j 2 J; k 2 K X l2L p l; k y l; j 1 X i2I C i; k z i; j P U j; k ð X l2l y l; j 1 X i2I z i; j Þ 8j 2 J; k 2 K 8 > > > < > > > : Ben-Tal et al. 50 develop a mathematically equivalent Qformulation that introduces fractional flow rates q i; l 5x i; l = Pî 2I xî ; l for arcs between an input i and pool l: min y; z; q X i2I; l2L; j2J c i q i; l y l; j 2 X l2L; j2J d j y l; j 2 X i2I; j2J Product Quality X i; l C i; k q i; l y l; j 1 X i2I C i; k z i; j P L j; k ð X l2L y l; j 1 X i2I z i; j Þ 8j 2 J; k 2 K X i; l C i; k q i; l y l; j 1 X i2I C i; k z i; j P U j; k ð X l2L y l; j 1 X i2I z i; j Þ 8j 2 J; k 2 K The P-and Q formulations have variable bounds: Flow from input i to pool l q i; l Flow from input i to pool l, as a fraction of total flow into l v i; l; j Flow from input i to output j through pool l (v i; l; j 5q i; l y l; j ) y l; j Flow from pool l to output j z i; j Bypass flow from input i to output j p l; k Level of quality attribute k in pool l Parameters c i Unit cost of raw material feed stock i d j Unit revenue of product j Combining the P-and Q-formulations, Quesada and Grossmann 43 develop the PQ Cut and Tawarmalani and Sahinidis 44 justify its utility. The PQ-formulation is equivalent to the Qformulation with the addition of the PQ Cut: X i2I q i; l y l; j 5y l; j 8l 2 L; j 2 J (PQ Cut)

Canonical Representation for the PQ-Formulation
Beyond the P-, Q-, and PQ-formulations, there are several other mathematically equivalent pooling problem formulations. 51,52 But we want to automatically translate every pooling problem into just one formulation (here, the proposed canonical PQ-formulation) and save implementation time by specializing methodologies such as cutting planes to just one formulation.
The proposed canonical PQ-formulation is equivalent to a PQ-formulation that introduces auxiliary variables v i; l; j and expands all product terms q i; l y l; j into v i; l; j . 17 The formulation also uses the PQ Cut to replace each flow y l; j with the relative sum of flows v i; l; j in every equation except for the PQ Cut and the definition v i; l; j 5q i; l y l; j . After rearranging terms, the canonical PQ-formulation becomes: The canonical PQ-formulation of the Figure 1 example pooling problem is visualized in Figure 11. Rows represent constraints and columns denote variables. The canonical PQformulation enables the pattern matching algorithm by allowing the algorithm to focus on linear terms; bilinear terms only appear in the PQ Path Definition constraints. Reformulating products q i; l y l; j using v i; l; j implies that identifying pooling constraints and building the network can be completed by primarily analyzing linear equations. Expressing the material capacity and product quality constraints with respect to only variables v and z simplifies the network building step and allows us to deduce parameter values for C i; k ; P L j; k ; P U j; k . The canonical PQ-formulation would perform badly in a branch-and-bound framework due to possibly weak relaxations and poor interval arithmetic. 34,37,46 But this article aims to uncover structure in preprocessing, so it is sensible to  Steps required to create a canonical PQformulation when the initial problem is expressed as a P-or Q-formulation; formulations detected via pattern matching are transferred to canonical form.
reformulate to a representation allowing easy network topology discovery. Once the algorithm deduces underlying structure, it may reformulate toward representations designed for whatever algorithm will best solve the problem.

Developing a Pattern Matching Method for Uncovering Pooling Problem Structure
This section describes how to: (1) reformulate a pooling problem to the canonical PQ-formulation, (2) identify the pooling problem constraints in an optimization problem, and (3) build a pooling problem network. We introduce the three steps with respect to finding patterns; and subsequently describe implementing the pattern matching algorithm using functional programming. Note that the pattern matching algorithm creates canonical pooling problems at the same time as it identifies pooling problem constraints, so the manuscript sections are coupled.
We take, as input, optimization problems with flat structure, that is, optimization problems where the variables and constraints are represented only by positions in a vector and have no associated identifier. The following discussion develops as if every constraint and/or variable will fit into some pooling meaning or another. But, under the hood, the implementation splits the MINLP into two parts, a pooling problem and a remaining optimization problem. The algorithm first assumes that every constraint and variable does not belong to the pooling problem and then, if it proves that the constraint or variable fits a known pattern, moves the variable or constraint to the pooling problem.
Converting to the canonical PQ-formulation Figure 2 outlines steps transforming an optimization problem toward the canonical PQ-formulation. The original problem formulation is not known a priori, so when the algorithm finds equations consistent with either a P or Q-formulation, it transforms toward an equivalent PQ-formulation. One difference between the P-and PQ-formulations is that the former uses flow rate variables nl whereas the latter introduces fractional flow rate variables q i; l for arcs between an input i and pool l such that nl and q i; l are related: To replace flow rates nl with q i; l , the algorithm first identifies which variables are nl; it does this by looking for equality constraints with constant coefficient 0 and all linear coefficients equal to either 11 or 21. Using the material balance constraints, i.e., P i2I x i; l 2 P j2J y l; j 50 8l 2 L, it labels variables with 11 and 21 coefficients as nl and y l; j , respectively. Then it introduces new variables q i; l for each variable nl and replaces each nl using Eq. 1.
The algorithm subsequently eliminates variables p l; k by noticing that, if it replaces nl in the Eq. 2 quality balance constraints, it can map each p l; k to the equivalent P i2I C i; k q i; l : The Q-formulation is attained by replacing variables p l; k with the relevant P i2I C i; k q i; l . From the Q-formulation, we obtain the PQ-formulation by taking all bilinear terms q i; l y l; j and, for each y l; j , adding the PQ Cut; this step is similar to existing algorithms. 37,39,47,48 The last step in transforming to a canonical PQ-formulation is replacing, using the PQ Cut, each occurrence of y l; j in the material capacity and product quality constraints with the sum P i2I v i; l; j .

Identifying pooling problem constraints
Coupled to transforming the optimization problem into the canonical PQ-formulation is the challenge of understanding each constraint and variable with respect to the pooling problem. This identification step enables the network extraction. Similar to the previous section, constraint and variable identification focuses on pattern matching. We show how to use the constraints to gain more problem information.
Identifying Reduction and Path Definition Constraints. . Path Definitions, v i; l; j 5q i; l y l; j , are easy to find. The algorithm uses the Path Definitions to find which y l; j is assigned to each v i; l; j and then identify Reduction Constraints: Specifically, we find Reduction Constraints by first grouping v i; l; j using Path Definitions, i.e., aggregate v i; l; j with the corresponding y l; j , and then look for constraints where path flows v i; l; j have coefficients 11 and y l; j has coefficient 21. Finding these Reduction Constraints helps identify variable subscripts since every v i; l; j has the same ðl; jÞ indices as the y l; j variable.
Identifying Input and Output Capacities. Next, the algorithm identifies what input i or output j indices correspond to each variable by differentiating the Input Capacity (Eq. 5) from the Output Capacity (Eq. 6) constraints. Recall that the algorithm has already identified the flow variables v i; l; j , Reduction Constraints (Eq. 3), and Pool Capacity constraints (Eq. 4), so the remaining v i; l; j sums with all coefficients 11 must be either Input or Output Capacity constraints: To identify Input Capacity constraints we define set vðl;ĵÞ5 v i; l; j jl5l; j5ĵ È É where all the variables v i; l; j are grouped by ðl; jÞ pairs. The algorithm easily constructs sets vðl; jÞ by using the Path Definitions to aggregate v i; l; j corresponding to the same y l; j . To find the Input Capacity constraints, observe that all v i; l; j variables in Eq. 5 are sums over ðl; jÞ and therefore each v i; l; j in the sum will correspond to a different ðl; jÞ pair. So the algorithm tags input capacity constraints as those where every variable v i; l; j belongs to a different set vðl; jÞ.
Identifying Output Capacity constraints (Eq. 6) is analogous to tagging the Input Capacity constraints; we define set vðî;lÞ We do not know the original network diagrammed in (a) but we can use, (b), the Path Definition constraints to group each flow y l; j with corresponding q i; l . All the sub-networks connected to pool l, (b), have identical y l; j network connections, so we obtain a sub network around each pool, (c). 5 v i; l; j ji5î; l5l È É and construct all the sets vði; lÞ by using the Path Definitions to aggregate v i; l; j corresponding to the same q i; l . Output Capacity constraints are sums over v i; l; j where every variable v i; l; j belongs to a different set vði; lÞ.
After finding the Input and Output Capacity constraints, the algorithm labels the extra variables appearing in Eqs. 5 and 6 as z i; j . If there are input or output nodes which are not connected to any pool, there may be unidentified inputs i or outputs j; the extra z i; j variables identify these input i or output j nodes which transmit or receive, respectively, bypass flows.
Identifying Quality Constraints. The only remaining equations to identify in the canonical PQ-formulation are the Product Quality constraints. Product Quality equations have the same variables as the Output Capacity constraints (Eq. 6), but the Output Capacities always have coefficient 11, while the Product Quality constraints each have coefficient C i; k 2 P L j; k or C i; k 2P U j; k . The algorithm tags the k quality constraints associated with output j by looking for Product Quality constraints where the same variables appear with different coefficients, an example is illustrated in Figure 5. The jKj3jJj product quality constraints are used to deduce input and output qualities, C i; k ; P L j; k ; P U j; k ; Section describes this procedure. We assume that each of the jKj qualities are given in the same order for every one of the Product Quality constraints.

Building the pooling problem network
Once all constraints and variables have been labeled, building the pooling network is almost trivial. We tag each network node with respect to an input i, pool l or output j; we also associate each node with (possibly infinite) capacity bounds. The algorithm also determines network flow arcs using path variables v i; l; j together with Eqs. 4-6 to find all inter-node connections. Figure 6 shows a path from input 2 to output 7 to pool 4 from the Figure 1 example.
Finding Input Cost and Output Revenue. The objective of a pooling problem is to minimize cost. After transformation to the canonical PQ-formulation, the objective has the form: where c i; j 5c i 2d j . We already know the variables v i; l; j and z i; j and their associated node indices, Linear Program (8) Finding Input and Output Qualities. Just as parameters c i ; d j are not directly in Eq. 7, pooling problem parameters C i; k ; P L j; k ; P U j; k are not directly present in the canonical PQformulation. To find these input and output quality parameters, define D L i; j; k 5C i; k 2P L j; k and D U i; j; k 5C i; k 2P U j; k and re-write the Product Quality constraints: Linear Program (11) deduces equivalent numerical values for C i; k ; P L j; k , and P U j; k . Here we assume that the indices of the input nodes i, output nodes j, and tracked qualities k are known. We also assume that a Product Quality constraint for Solving LP (11) may produce different values for C i; k ; P L j; k , and P U j; k than in the original problem, but the feasible space defined by the constraints is identical.

Implementation Sketch
This section sketches implementing the algorithm; the online supplementary material gives more details. We describe data structures, discuss IO facilities, show how to implement the algorithms.

Data structures
A pooling network is characterized by its inputs, pools, outputs, and arcs. Mathematical optimization problems are defined by a set of variables, an objective function, and a set of constraints; MINLP may also have an associated pooling network. The online supplementary material defines several modules and types which generically represent MINLP optimization problems. The most general module, Problem, has a record field to contain a PoolingNetwork. When the algorithm initiates, the PoolingNetwork field is None. After the implementation terminates, PoolingNetwork contains the standard pooling network and the constraints field holds additional equations irrelevant to a standard pooling framework.

Reading and writing problems
The implementation reads input problems in OSiL format; OSiL is a widely supported optimization format and is easy to parse since it is based on XML. 53 We provide two distinct output formats: the first writes the pooling network to a Graphviz file and the second writes the pooling network to a Dey and Gupte 17 formatted AMPL data file. The Graphviz file output, which may be subsequently converted to an image, generated the results shown in this article.

Transforming and identifying constraints
The implementation concurrently: (1) reformulates a pooling problem to the canonical PQ-formulation and (2) identifies the pooling problem constraints in an optimization problem. All functions transforming and identifying constraints follow a common pattern: scan every constraint and, if it matches a given pattern, execute an action. Figure 7 illustrates the procedure and the online supplement gives more specific details. Actions can be as simple as updating the constraint type or as complex as replacing a variable in the whole problem.

Building the network
The final implementation phase is building the network. We begin by building inputs, pools, and outputs from the constraints. As shown in Figure 8, node capacities are either the constraint upper bound or otherwise infinite. Finally, we need to find the price and quality specifications for both the inputs  and outputs; we find these parameters by solving LPs (8) and (11).

Results and Discussion
We tested the implementation on three sets of input These are the largest pooling problems solved to date and we therefore have effectively no size limitation for the detection methodology. Table 2 completely documents the results on the standard and extended pooling problems from Misener et al. 54 The This network is mathematically equivalent to the original Lee1 network which has 5 rather than 20 input nodes. The original network has five input nodes with infinite capacities; the implementation cannot distinguish this from five input nodes for each of the pools (all of which still have infinite capacities).  implementation has no problem excluding the nonlinear Environmental Protection Agency 15 blending rules from consideration. The pooling problems in this test set are less stylized and not consistently defined as the Dey and Gupte 17 test set, so observe that, although we get a pooling network for all the problems, it is not necessarily the original pooling network. Running our implementation on the 1342 test cases in MIN-LPLib2, three produce complete pooling problems and 78 more produce a pooling-like network. The first two problems are from Lee and Grossmann. 55 The original optimization problems are generalized pooling problems with mutable network topology, but the implementation finds an equivalent standard pooling problem. Figure 9 shows the network produced for Lee1; this is mathematically equivalent to the original Lee1 network because the input capacities are infinite. The third completely reconstructed problem is Haverly. The Haverly test case did not do as well in the Misener et al. 54 test suite; this result reminds us that our hunt for pooling problems within a larger optimization problem is ultimately a heuristic.
As an example of the 78 additional problems where the implementation finds a pooling-like network, we discuss the wastewater02m2 wastewater network problem. 56 As illustrated in Figure 10a, the wastewater problem no longer has a feed-forward network because flows exiting the treatment units (T1 and T2) may be sent through the treatment units again. The pooling problem we find in Figure 10b is exciting because it has flattened the wastewater treatment problem into a pooling problem. The two extra inputs and outputs illustrated in Figure 10b could be made equivalent to the original topology by using some of the additional constraints in the auxiliary problem, see Figure 10c. Note in Figure 10c that the number of pipes entering/exiting pools is the same as the number of pipes entering/exiting the treatment units in Figure 10a.
The 16 standard and extended Misener et al. 54 examples and the 78 MINLPLib2 test cases with network structure can be addressed with off-the-shelf solver software, 37 but the 70 large-scale standard pooling Dey and Gupte 17 examples are currently out of reach for solver software. To test the potential of using the Dey and Gupte 17 heuristic in tandem with our pattern-finding approach, we compare using performance profiles 57 : t p; s Performance metric; i:e:; best feasible solution; for problem p by solver s 2 S r p; s t p; s 2 min p2P t p; s 0 : s 0 2 S È É min p2P t p; s 0 : s 0 2 S È É ; s 2 S q s ðsÞ5 1 n p size p 2 P : r p; s s È É After finding the structure in the Dey and Gupte 17 test set, we initialized ANTIGONE, and several of the Dey and Gupte 17 approximation algorithms (U1, U2, U4, U5) and ran each for 30 min. Figure 12 shows the advantage of using the new heuristic in comparison to the built-in ANTIGONE heuristics for this problem class. Although our pattern-matching algorithm coupled with the Dey and Gupte 17 heuristic is much more powerful than the heuristics currently in ANTIGONE, we have not added the new pattern matching code to ANTIGONE because, at this point, the pattern-matching time significantly slows down ANTIGONE performance on non-pooling problems. But we are still very interested in these new results because we know that there are several cutting plane approaches for the pooling problem under development. 18,19 Once both cutting plane and heuristic approaches are well-developed for the pooling problem, finding large pooling problems will be even more advantageous. Furthermore, this article proves that optimization researchers can easily continue studying the more stylized standard pooling problem since we can find pooling structure within larger, flat MINLP.

Conclusions
This article has explored, for the first time, the possibility of finding special named structure within an MINLP optimization problem; this significantly extends work by the MIP community that has shown how to find network structure for better primal heuristics and cutting planes. We are specifically interested in the range of novel techniques which have been developed for the standard pooling problem and wanted to see how far we could get in recognizing pooling structure within a general MINLP. We have shown that we can detect all standard and extended pooling problems in the literature and that we can find pooling networks in 6% of MINLPLib. We have also shown that, after detecting the pooling problem in the flat MINLP, we can apply a good heuristic approach to get a good approximation solution.