Advances in MINLP to Identify Energy-efficient Distillation Configurations

In this paper, we describe the first mixed-integer nonlinear programming (MINLP) based solution approach that successfully identifies the most energy-efficient distillation configuration sequence for a given separation. Current sequence design strategies are largely heuristic. The rigorous approach presented here can help reduce the significant energy consumption and consequent greenhouse gas emissions by separation processes, where crude distillation alone is estimated to consume 6.9 quads of energy per year globally. The challenge in solving this problem arises from the large number of feasible configuration sequences and because the governing equations contain non-convex fractional terms. We make several advances to enable solution of these problems. First, we model discrete choices using a formulation that is provably tighter than previous formulations. Second, we highlight the use of partial fraction decomposition alongside Reformulation-Linearization Technique (RLT). Third, we obtain convex hull results for various special structures. Fourth, we develop new ways to discretize the MINLP. Finally, we provide computational evidence to demonstrate that our approach significantly outperforms the state-of-the-art techniques.


Introduction
Separation of mixtures of chemical components is ubiquitous in all chemical and petrochemical industries. Among the numerous technologies available for separation of multicomponent mixtures (three or more components) into almost pure components, distillation is the predominant choice. A few well-known applications include fractionation of crude oil (Sholl and Lively, 2016), production of ultra pure nitrogen and oxygen from air (Agrawal and Woodward, 1991), Natural Gas Liquid (NGL) recovery, Benzene-Toluene-Xylene (BTX) separation, etc. It is estimated that distillation accounts for 90 -95% of the liquid phase separations in the US (Humphrey, 1997). With the increased potential to harness shale reserves (Siirola, 2014;Ridha et al., 2018), the use of distillation is projected to increase further. Industrial distillations are energy intensive, and the energy consumed constitutes about 40 -60% of the total operating cost (Humphrey, 1997). Since energy consumed affects the effective fuel consumption, suboptimal configurations also tend to release more CO 2 . In this article, we develop the first tractable approach to solve distillation sequencing via Mixed-integer Nonlinear Programming (MINLP) techniques.
Mixtures are separated in a distillation configuration consisting of a series of distillation columns/towers (see Figure 1) arranged to carry out the separation in a specific order (see Figure 2 for an example). The number of admissible configurations grows rapidly with the number of components in the given mixture. For example, over half a million alternative configurations are admissible for separating a six component mixture ( Figure 1: Schematic of a conventional distillation column. Double-lines in brown, known as trays/stages, establish contact between the vapor (red arrows) and liquid (blue arrows) for mass transfer. Here, and in the rest of the article, condensers and reboilers are denoted by filled and open ellipses (or circles), respectively. Agrawal, 2010). Besides the number of choices, the nonconvex fractional terms used to model the minimum energy requirement make this problem hard to solve. Prior to this work, the state-of-the-art methods are unable optimize the material flows even for a specific configuration. As such, conventional design practices are based on heuristics and intuition of the process engineer, and they often result in suboptimal solutions.
To address these challenges, we develop a new mixed-integer nonlinear programming (MINLP) based approach and make several advances. First, we develop a new model for the space of admissible configurations, by incorporating convex hulls of various substructures. The resulting formulation is provably tighter than those in the literature (Caballero and Grossmann, 2006;Giridhar and Agrawal, 2010b;Tumbalam Gooty et al., 2019). Second, we show that polynomial division and partial fraction decomposition can significantly improve the quality of relaxations developed using the classical Reformulation-Linearization Technique (RLT) (Sherali and Alameddine, 1992) when the constraints involve fractional terms. Third, we develop simultaneous convex hulls of multiple nonlinear terms over a polytope obtained by intersecting bounds on variables with material balance equations. Fourth, we derive the first rigorous relaxation for distillation sequencing. The governing equations for this problem involve fractions, whose denominator can approach zero. To sidestep this issue, the literature has imposed an ad hoc bound on the denominator. This approach, however, can prune optimal solutions when some component flows are small. Instead, using the cuts from the RLT variant described above, our approach derives a rigorous relaxation for the problem.
In §2, we briefly describe the key concepts of multicomponent distillation, and survey the current literature. §3 defines the problem statement and introduces relevant notation. We formulate the MINLP in §4, and outline the overall relaxation and solution procedure in §5. We report on computational experiments in §7. Finally, we make concluding remarks in §8.

The Distillation Process
Distillation is a way to separate mixtures, consisting of two or more components with different relative volatilities, by boiling the mixture so that the vapor produced is rich in more volatile (or light) components, while the residual liquid is enriched in less volatile (or heavy) components. Industrial distillation is carried out in a staged-tower/column (see Figure 1), where each stage establishes liquid-vapor contact for mass transfer. The feed (mixture of components) is introduced at an intermediate location of the column. The sections above and below the feed stream are known as rectifying and stripping sections, respectively. Conventional columns have a condenser (resp. reboiler) at the top (resp. bottom) which condenses (resp. vaporizes) the vapor (resp. liquid), and feeds a portion of it back to the column, known as liquid (resp. vapor) reflux. The liquid flowing from the top to bottom strips away heavy components from the vapor, while the vapor flowing from bottom to top gets enriched with lighter components. The net outflow from the rectifying and stripping sections, respectively, are known as distillate and residue. In short, distillation enriches the distillate with light components, and the residue with heavy components. Figure 2: Example of a distillation configuration. C i . . . C j denotes a mixture, and each C p corresponds to a distinct chemical component. C 1 . . . C 4 is the process feed, and intermediate mixtures C 1 C 2 C 3 , C 2 C 3 C 4 and C 2 C 3 are referred as submixtures. (Nallasivam et al., 2016;Mathew et al., Working paper).

Remark 1. The recovery of a lighter component in distillate (ratio of component flowrate in distillate to flowrate in feed) is higher than the recovery of a heavier component, and the converse is true for residue
A given N´component mixture, referred to as the process feed, is separated into N constituent components using a sequence of distillation columns (see Figure 2, for example). Let C i . . . C j denote an intermediate stream (referred to as submixture), where components are sorted in a decreasing order of relative volatilities, and C p denotes the p th component in the process feed. Each column splits a feed submixture into two product submixtures, each of which has at least one component less than the feed. The composition of the product submixtures governs the threshold vapor requirement for the column. This requirement can be determined using the classical Underwood method (Underwood, 1948), as long as relative volatilities are constant, each section has infinite stages, and there is constant molar overflow. As shown in Figure 2, condensers and reboilers can be replaced with two-way vapor-liquid transfer streams known as thermal couplings, so that the required liquid/vapor reflux is borrowed from other columns. Since thermal couplings allow vapor to be transferred between two or more columns, a column may be operated above its threshold vapor requirement to supply the vapor to another column. We remark that configurations with many thermal couplings may be hard to control (Agrawal, 2000) and require hot and/or cold utilities at extreme temperatures. Although we do not explicitly model these issues, configurations with few thermal couplings can easily be found by simple changes to our formulation.
For above ambient distillation, the required vapor flow is generated at reboilers by a hot utility. By adding these vapors, we obtain the vapor duty of the configuration, which is often used as a proxy for its energy consumption and operating cost. The vapor duty indirectly affects the capital cost as well, since internal vapor flows dictate column diameters. For these reasons, we will minimize vapor duty, an objective that has also been used in previous studies (Fidkowski and Królikowski, 1987;Fidkowski and Agrawal, 2001;Nallasivam et al., 2016). Industrial practitioners may instead be interested in minimizing the total annualized cost (capital plus operating costs), or maximizing the thermodynamic efficiency. The model we propose can be tailored to the desired objective by appending the relevant constraints and modifying the objective as in Jiang et al. (2019a,b).
Given its importance, this problem has been studied extensively, but has resisted formal solution guarantees. Grossmann (2004, 2006) formulated an MINLP to identify configurations with lowest total annualized cost, but did not certify global optimality. Given the non-convexity, these local approaches do not always find optimal designs (Nallasivam et al., 2013;Jiang et al., 2019b). Giridhar and Agrawal (2010b) proposed an alternate MINLP formulation to minimize vapor duty, and solved it using BARON (Tawarmalani and Sahinidis, 2005) for three and four-component mixtures. However, their methodology does not scale to five component mixtures. Nallasivam et al. (2016) enumerated all configurations, and solved a nonlinear program for each using BARON. However, some configurations fail to converge. The current state-of-the-art formulation of Tumbalam Gooty et al. (2019) still fails to converge on 36% of the MINLP instances. Figure 3 shows all possible streams and heat exchangers in a distillation configuration that separates a fourcomponent mixture into pure components. We represent streams as squares, condensers as filled circles and reboilers as open circles. Each condenser/reboiler is associated with a process stream, that is not the process feed C 1 . . . C N . Throughout the formulation, we denote a stream C i . . . C j as ri, js, and heat exchangers as pi, jq, so that condenser pi, jq (resp. reboiler pi, jq) represents the heat exchanger through which ri, js is withdrawn as distillate (resp. residue). By Remark 1, a configuration cannot contain streams of the form C i . . . C k C k`l . . . C j , where l ą 1.

Problem Definition
[  , reboilers (ρ i,j ) and condensers (χ i,j ) present in a four-component system. Section variables τ i,k,j and β i,l,j are defined in (1).
We denote the set of streams as T , the set of condensers as C, and the set of reboilers as R (see Table 1 for definition). For convenience, we create a set containing streams that are mixtures P " T ztr1, 1s, . . . , rN, N su, and a set containing submixtures S " Pztr1, N su. Note that every stream in P is a mixture, and must undergo a split in order to produce products.
The required input to the problem consists of (1) composition of the process feed tF p u N p"1 either in terms of mole fractions or molar flowrates of the components in the stream, (2) relative volatilities tα p u N p"1 (such that α N ă¨¨¨ă α 1 ) of its constituent components; and (3) liquid fraction (fraction of the total flow in liquid phase) of the process feed Φ 1,N and that of the pure components tΦ i,i u N i"1 . We write tpu N p"1 or tpu 1ďpďn to denote the set t1, 2, . . . , N u, and p N 1 to denote @ p P t1, . . . , N u. Given a process feed, the problem is then to identify the best distillation configuration, along with its optimal operating conditions, that requires least vapor duty.

Set Symbol Definition
Streams T tri, js : 1 ď i ď j ď N u Reboilers R tpi, jq : 2 ď i ď j ď N u Table 1: Definition of sets.

Problem Formulation
We formulate the MINLP in this section. Before proceeding further, we introduce the definition of parents and children of a stream. By top (resp. bottom) parents of ri, js: we refer to streams tri, nsu N n"j`1 (resp. trm, jsu i´1 m"1 ) which can produce ri, js as distillate (resp. residue). Analogously, by top (resp. bottom) children of ri, js, we refer to streams tri, ksu j´1 k"i (resp. trl, jsu j l"i`1 ) which can be produced as distillate (resp. residue) from ri, js. For conciseness, we write ri, js Ò ri, ks (resp. ri, js Ó rl, js) to denote stream ri, ks (resp. rl, js) is produced as the distillate (resp. residue) from ri, js, and ri, ks{rl, js to denote ri, ks and rl, js are produced as the distillate and residue from ri, js.

Objective Function
The objective is to determine the configuration(s) which minimizes the total vapor duty: where FR i,j is the vapor flow generated in reboiler pi, jq. The MINLP we develop will be denoted as MINLP (A), and the constraints will be numbered as (A#).

Space of Admissible Configurations
We define column/stream binary variables so that @ ri, js P T , ζ i,j " 1 if ri, js is present and 0 otherwise. Further, we define binary variables associated with the presence/absence of condensers and reboilers so that @ pi, jq P C (resp. @ pi, jq P R), χ i,j " 1 (resp. ρ i,j " 1) if condenser (resp. reboiler) pi, jq is present and 0 otherwise (See Table 1 for set definitions). Although these variables suffice (Tumbalam Gooty et al., 2019), we introduce auxiliary variables to derive a tighter representation.
For every ri, js P P, we define section variables tτ i,k,j u j´1 k"i and tβ i,l,j u j l"i`1 , such that τ i,k,j " t1, if ri, js Ò ri, ks; 0, otherwiseu and β i,l,j " t1, if ri, js Ó rl, js; 0, otherwiseu. In other words, section variables model distillate and residue streams from a mixture. Figure 3 shows all the section variables for a four-component mixture. We now relate column and section variables. Consider the split of stream ri, js. In configurations of interest, known as regular-column configurations, if ri, js Ò ri, ks, for any i ď k ď j´1, then ri, js and ri, ks must be present and tri, nsu j´1 n"k`1 must be absent (Caballero and Grossmann, 2006;Giridhar and Agrawal, 2010b). Analogously, if ri, js Ó rl, js, for any i`1 ď l ď j, ri, js and rl, js must be present, while trm, jsu l´1 m"i`1 must be absent. Therefore, section variables are defined as A PREPRINT -OCTOBER 26, 2020 We introduce variables tν i,k,j : 1 ď i ď k ď j ď N u and tω i,l,j : 1 ď i ď l ď j ď N u to linearize (1): for ri, js P P # τ i,k,j " ν i,k`1,j´1´νi,k,j´1´νi,k`1,j`νi,k,j , k j´1 i β i,l,j " ω i`1,l´1,j´ωi,l´1,j´ωi`1,l,j`ωi,l,j , where ν i,k,j " ś j n"k p1´ζ i,n q and ω i,l,j " ś l m"i p1´ζ m,j q. Note that ν i,k`1,j´1 (resp. ω i`1,l´1,j ) are defined as one if k`1 " j (resp. i`1 " l). Clearly, ν i,k,j " ω i,l,j " 1´ζ i,j if k " j and l " i. Besides this relationship, the introduced variables ν i,k,j and ω i,l,j are linearly independent. To see this, note that ś jPJ x j , where J Ď t1, . . . , nu are linearly independent and, therefore, so are ś jPJ p1´y j q, where y j " 1´x j . Since ν i,k,j and ω i,l,j are of the latter form, they are linear independent.
We now relax ν i,k,j and ω i,l,j variables for k ‰ j and l ‰ i as follows. Since ζ i,j is binary, p1´ζ i,j q 2 " p1´ζ i,j q. We use the definition of ν i,k,j and ω i,l,j , to derive the following: In the above, for n ď m`1, ν i,n,m (resp. ω n,m,j ) is a common factor for both ν i,k,m and ν i,n,j (resp. ω i,m,j and ω n,l,j ). we regard ν i,n,m and ω n,m,j as one if n " m`1. Thus, 0 ď ν i,k,m ď ν i,n,m , 0 ď ν i,n,j ď ν i,n,m , 0 ď ω i,m,j ď ω n,m,j , and 0 ď ω n,l,j ď ω n,m,j . Using these bounds, we relax (2) as: where we used ν i,k,m " ν i,k,m ν i,n,m , ν i,n,j " ν i,n,j ν i,n,m , ω i,m,j " ω i,m,j ω n,m,j , and ω n,l,j " ω n,l,j ω n,m,j .
Proof. See §A in the appendix.
We remark that the result in Proposition 1 also follows from Theorem 10 in Del Pia and Khajavirad (2018). Our proof is, however, different and elementary. We mention that this proof shows a previously unobserved connection to the recursive McCormick procedure. Our proof can be used to show that the recursive McCormick procedure, with a few additional linearization variables, yields the convex hull of the multilinear polytopes for γ-acyclic hypergraphs, as obtained in Del Pia and Khajavirad (2018).
We now describe the constraints to model the space of admissible distillation configurations.
To restrict the search to a subset of configurations, for example, in order to retrofit an existing design, we may explicitly include (resp. eliminate) a specific submixture ri, js by setting ζ i,j " 1 (resp. ζ i,j " 0). We show next that ζ i,j variables are affinely related to τ i,k,j and β i,l,j variables.

Proposition 2.
Let, x P r0, 1s n , y i,j " p1´x i qx i`1 . . . x j´1 p1´x j q for 1 ď i ă j ď n, z i,j " ś j r"i x r for 1 ď i ď j ď n, and x n " 0, which in turn implies that z i,n " 0 for 1 ď i ď n. Then, there is an invertible affine transformation between ty i,j u 1ďiăjďn and tz i,j u 1ďiďjďn , given by y i,j " z i`1,j´1´zi`1,j´zi,j´1`zi,j , z p,q " 1´q ÿ r"p n ÿ s"q`1 y r,s .
Proof. First, we show that y i,j can be written as an affine transformation of z i,j . By definition, y i,j " p1´x i qx i`1 . . . x j´1 p1´x j q " ś j´1 r"i`1 x r´ś j r"i`1 x r´ś j´1 r"i x r`ś j r"i x r " ś j´1 r"i`1 x r´zi`1,j´zi,j´1`zi,j . Substituting the first term in the last equality, ś j´1 r"i`1 x r , with 1 if i`1 " j, and z i`1,j´1 if i`1 ă j, yields the required affine transformation.
Next, to obtain the inverse affine transformation, we define w k,l " p1´x k qx k`1 . . . x l for 1 ď k ď l ď n. We show the affine transformation between tw k,l u 1ďkďlďn and ty i,j u 1ďiăjďn variables to be using induction on n´l. For l " n, (3) is trivially satisfied because w k,n " 0 as x n " 0. Now, assuming that (3) holds for l`1, i.e., w k,l`1 " ř n r"l`2 y k,r , we show that it holds for w k,l as well: w k,l " p1´x k qx k`1 . . . x l p1x l`1`xl`1 q " y k,l`1`wk,l`1 " y k,l`1`ř n r"l`2 y k,r " ř n r"l`1 y k,r . In a similar vein, we show for 1 ď p ď q ď n, the affine transformation between tz p,q u and tw k,l u variables to be using induction on q´p. For q " p, (4) follows because z q,q " x q " 1´p1´x q q " 1´w q,q . Next, assuming (4) holds for p`1 i.e., z p`1,q " 1´ř q r"p`1 w r,q , we show that it holds for z p,q as well: z p,q " ś q r"p x r " r1´p1´x p qs ś q r"p`1 x r " z p`1,q´wp,q " 1´ř q r"p`1 w r,q´wp,q " 1´ř q r"p w p,q . Finally, substituting (3) in (4) leads to the required inverse affine transformation given below: Indeed, the correctness of (5) can be checked via direct verification using y r,s " z r`1,s´1´zr`1,s´zr,s´1`zr,s , z i,n " 0 for 1 ď i ď n, and z i`1,i " 1 for 1 ď i ď n.
ř j l"i`1 β i,l,j " ζ i,j for all ri, js P P.

Conservation of components
Corollary 1 has the physical interpretation that the stream ri, js, when present, produces exactly one stream as distillate and one stream as residue. However, the distillate and residue streams cannot be chosen arbitrarily. They must be chosen such that, all components are conserved when ri, js undergoes a split. In other words, for k j´1 i (resp. l j i`1 ), if ri, js Ò ri, ks (resp. ri, js Ó rl, js), then for conservation of components, the residue (resp. distillate) from ri, js must be one of trl, jsu k`1 l"i`1 (resp. tri, ksu j´1 k"l´1 ). Consider the digraph shown in Figure 4 for stream ri, js.
Complete bipartite graph with ψ i,N +1,0,j = 0 We partition the nodes into four sets D 1 through D 4 , where D 1 " tiu (resp. D 4 " tju), and D 2 " tku j´1 k"i (resp. D 3 " tlu j l"i`1 ) contains the heaviest (resp. lightest) component in the top (resp. bottom) children of ri, js. The edges in D 1ˆD2 (resp. D 3ˆD4 ) correspond to all plausible distillate (resp. residue) streams from ri, js. Edges in D 2ˆD3 correspond to feasible splits of ri, js, i.e., each node k P D 2 connects to ti`1, . . . , k`1u P D 3 . We associate these edges with auxiliary variables Ť j´1 k"i tσ i,k,l,j u k`1 l"i`1 , referred as split variables hereafter (see Figure 4). We let σ i,k,l,j " t1, if ri, ks{rl, js; 0, otherwiseu, and write mass balances on the network by interpreting stream, section and split variables as material flows along the respective edges of the graph.
For ri, js P P Mass balances around the nodes in D 1 and D 4 , and non-negativity constraint on section variables are implied from (A2)-(A4) (see Corollary 1 and Remark 3), so it is not required to impose them explicitly. We show below that, for any ri, js P P, the relaxation (A2)-(A5) is the best possible for the substructure represented by the digraph in Figure 4.
As in §4.2.2, we partition the nodes into four sets D 5 through D 8 (see Figure 5), where D 5 " tiu (resp. D 8 " tju), and D 6 " tnu N`1 n"j`1 (resp. D 7 " tmu i´1 m"0 ) contains the heaviest (resp. lightest) component in the top (resp. bottom) parents of ri, js. Recall that m " 0 and n " N`1 have a special meaning as described in the previous paragraph. The edges in D 5ˆD6 (resp. D 7ˆD8 ) correspond to all plausible ways ri, js can be produced as distillate (resp. residue), and the edges in D 6ˆD7 indicate whether ri, js is produced only as distillate or only as residue or both. We introduce variables for edges in D 6ˆD7 such that ψ i,n,m,j " 1 iff ri, ns Ò ri, js and rm, js Ó ri, js.
A PREPRINT -OCTOBER 26, 2020 We require that ψ i,N`1,0,j " 0, which, otherwise, would mean that ri, js can be present even if it is neither produced as distillate nor as residue. Now, we write mass balances on the network.
for ri, js P T ztr1, N su Mass balances around the nodes in D 5 and D 8 , and non-negativity constraint on section variables are implied from (A2) and (A3), so it is not required to impose them explicitly.
Proposition 4. The constraints (A2), (A3), (10) and 0 ď ζ i,j ď 1 define a set such that, for every ri, js P T ztr1, N su, pτ, β, ζ, ψq is contained in the convex hull of Proof. We use a similar argument as the one used to prove Proposition 3. We recognize that (10), ř N`1 n"j`1 τ i,j,n " ř i´1 m"0 β m,i,j " ζ i,j , 0 ď ζ i,j ď 1 and non-negativity requirement on section variables together constitute a network flow polytope, whose extreme points are integral and precisely those in S i,j .

Constraints on Heat Exchanger Variables
For every pi, jq P C, condenser pi, jq is present only if the stream ri, js is not produced as residue, i.e., β 0,i,j " 1 (Tumbalam Gooty et al., 2019). Similarly, for every pi, jq P R, reboiler pi, jq is present only if the stream ri, js is not produced as distillate, i.e., τ i,j,N`1 " 1. Further, a condenser (resp. reboiler) must be present with a pure component ri, is, if ri, is is not produced as residue (resp. distillate) i.e. β 0,i,i " 1 (resp. τ i,i,N`1 " 1).

all the variables are zero.
Proof. We modify the graph in Figure 5 to accommodate (A6) and (A7), and combine it with the graph in Figure 4. The resulting graph is shown in Figure 6. Next, observe that (A5), (10), n"j`1 τ i,j,n " ř i´1 m"0 β m,i,j " ζ i,j , 0 ď ζ i,j ď 1 (which are implied from (A2)-(A4)), and non-negative constraint on all variables together constitute a network flow polytope. The extreme points this polytope are integral, and are precisely those mentioned in the Proposition.
Since ψ variables are not used elsewhere, we project (10) to the space of section variables (τ, β). Figure 6: Digraph for the proof of Proposition 5 Proof. See §C in the Appendix.
Proof. See §D in the Appendix.
The fact that our formulation is strictly tighter will follow from numerical examples.

Mass Balance Constraints
We model the problem as a network flow problem. Figure 7 shows the representative nodes and arcs in the network, and variable definitions are in Table 2. Each split ri, ks{rl, js is performed in a distillation column Q iklj (see Figures 7(a) and 7(b)). Material flows to and from the column Q iklj only when σ iklj " 1. The material balances across each column Q iklj are as follows for ri, js P tr1, N su, k j´1

liquid and vapor flowsu
, for ri, js P P, k j´1 The constraints in (A9) model component, vapor, and liquid mass balances across column Q iklj . In the above δ p¨q is 1 if p¨q is true and 0 otherwise. (A10) handles the case where the feed stream is the process feed, r1, N s. F p and Φ 1,N are as defined in §3. The last constraint in both (A9) and (A10) suppresses material flows to column Q iklj when σ iklj " 0. We use p¨q up to denote the upper bound on p¨q, and discuss how these are obtained later. The first (resp. second) constraint in (A11) models that the net distillate (resp. residue) flow Q iklj as the difference between the vapor and liquid (resp. liquid and vapor) flows in the rectifying (resp. stripping) section. Vapor out-flow from Q iklj to reboiler pi, jq   Figure 7: (a) Representative column for splits of process feed i.e., ri, js P tr1, N su, k j´1 i , l k`1 i`1 (b) Representative column for the remaining splits ri, js P S, k j´1 i , l k`1 i`1 (c) Representative condenser for pi, jq P Cztri, isu N´1 i"1 (see (A12) for domain of indices m, n, k, l) (d) Representative reboiler for pi, jq P Rztri, isu N i"2 (see (A13) for domain of indices m, n, k, l) (e) Representative arrangement for pure product withdrawals (see (A14) for domain of indices m, n, and (A15) for domain of indices m 1 , n 1 ) (f) Representative arrangement for overall component mass balance for ri, js P S (see (A17) for domain of indices m, n, m 1 , n 1 , k, l)

Variable Definition
Column Q iklj receives feed from the associated condenser pi, jq and/or reboiler pi, jq (see Figures 7(c) and 7(d)). Further, condenser (resp. reboiler) pi, jq regulates vapor-liquid traffic from all the splits producing ri, js as distillate (resp. residue), and distributes flows to all the splits of ri, js. Material balances across these condensers and reboilers are given below: We are interested in configurations that either have heat exchangers or thermal couplings, but not both. The last two constraints in (A12) and (A13) suppress flows in appropriate arcs if the heat exchangers are absent. The above constraints are written only for heat exchangers associated with mixtures. For heat exchangers associated with pure component products, the vapor and liquid flows are further constrained to produce Φ i,i F i and p1´Φ i,i qF i of component i in liquid and vapor phases, respectively (see Figure 7(e)). Mass balances around these heat exchangers are given below.
For pi, iq P C : For pi, iq P R : where p¨q lo denotes the lower bound on p¨q. From (A16) and (A15) (resp. (A14)), pK rs iiii q lo "´F p Φ i,i (resp. pU ss iiii q lo "´F p p1´Φ i,i q). For each submixture ri, js P P, the net inflow of component p equals the sum of component flows from all the splits that produce ri, js as distillate or residue. The net inflow is distributed among all splits of ri, js (see Figure 7(f)).
For ri, js P S : Finally, modeling the problem in the above manner requires rigorous bounds on all material flows. The net component inflow to and outflow from any column cannot exceed in steady-state the component flow in the process feed. Therefore, the upper bound on all flows of component p is chosen to be F p i.e., pf in ikljp q up " pf rs ikljp q up " pf ss ikljp q up " F p . However, although required for deriving rigorous relaxations, there is no simple upper bound on vapor and liquid flows in the columns and heat exchangers. For deriving a bound, we use optimality-based bound tightening, where we find feasible flows for an admissible configuration using the technique of Nallasivam et al. (2013). This technique can also be replaced with a local nonlinear programming solver. Let this upper bound be VD˚. Then, we solve the following linear programs (LP) to derive bounds: ( 13) We choose φ " 1, if only the optimal solution is desired. Since the model does not capture all operability concerns, such as controllability and suitability w.r.t heat integration with the rest of the plant, and vapor flow predictions are based on shortcut methods rather than rigorous simulations, industrial practitioners are often interested in identifying a ranklist of a few best solutions for this MINLP. Such a ranklist allows them to a posteriori incorporate such considerations. Therefore, to allow construction of such a ranklist, we choose φ " 1.5. With this choice, any configuration that consumes at most 50% more energy than the feasible solution remains in the search space. Our numerical experiments show that each LP can be solved in a fraction of a second using solvers such as Gurobi (Gurobi Optimization, 2018), and the computational time taken to solve all the LPs for a five-component mixture is typically negligible.

Underwood Constraints
As mentioned in §2, for a given split, there is a minimum threshold vapor requirement in each section of a column, below which the products are not produced with the desired purity. A column can, however, carry more vapor than the threshold, and the excess vapor can, if transferred to other columns, be utilized in those columns. The threshold vapor requirement can be computed using Underwood constraints included below: where Υ rs iklj and Υ ss iklj denote the threshold vapor flow in rectifying and stripping sections, respectively. Note that, for the process feed r1, N s, f in ikljp and U rs iklj δ jăN´U ss 1ăiklj δ 1ăi in (A18) and (A22) are replaced by F p σ iklj and´ř N p"1 F p¯p 1´Φ 1,N qσ i,k,l,j , respectively. (A18) is commonly known in the literature as the Underwood feed equation, and it computes Underwood roots tθ ijq u k q"l´1 , which satisfy α q`1 ď θ ijq ď α q (Underwood, 1948). (A19) governs the minimum vapor requirement in rectifying and stripping sections as a function of the distillate and residue compositions. (A20) ensures that the minimum vapor constraints are binding for tθ ijq u k´1 q"l . These constraints are required for the model to have the correct degrees of freedom as described in Tumbalam Gooty et al. (2019). (A22) models vapor balance at the feed location in terms of minimum vapor flows. (A23) ensures that the actual vapor in each section is at least as high as the threshold vapor flow. Remark 5. Since the process feed is always present i.e., ζ 1,N " 1, and the net component and vapor inflow to columns Q 1klN where 1 ă l ď k`1 ď N are known, we solve the Underwood feed equation (A18) a priori to determine the Underwood roots tθ 1N q u N´1 q"1 , and fix these variables to the calculated values.
Next, using (14), component mass balance f in ikljp " f rs ikljp δ pďk`f ss ikljp δ pěl , and (A18), it can be shown that Since the vapor flows are bounded, we have finite upper and lower bounds on all nonlinear expressions in (A18)-(A20).

Exploiting Monotonicity of Underwood Equations
These cuts are inspired from Carlberg and Westerberg (1989) and Halvorsen and Skogestad (2003a). Although these relations are implicit in the model, they are not implied in the relaxation, when Underwood constraints are relaxed. We refer to Tumbalam Gooty et al. (2019) for a derivation.
When ri, js is produced as distillate from one of its top parent ri, ns where j`1 ď n ď N i.e., τ i,j,n " 1, but not produced as residue from any of its bottom parents i.e., β 0,i,j " 1, and the associated condenser pi, jq is absent, then θ inq lower bounds θ ijq for q j´1 i . Similarly, when ri, js is produced as residue from one of its bottom parent rm, js where 1 ď m ď i´1 i.e., β m,i,j " 1, but not produced as distillate from any of its top parents i.e., τ i,j,N`1 " 1, and the associated reboiler pi, jq is absent, then θ mjq upper bounds θ ijq for q j´1 i . These constraints are imposed as follows: where M q " pα q´αq`1 q corresponds to the upper bound on the difference of Underwood roots (see (A21)). Numerical examples in Tumbalam Gooty et al. (2019) illustrate that these cuts help branch & bound converge faster. Given that our formulation has been developed in a lifted space, we use τ and β variables to give a tighter representation of the constraint in (A24). Moreover, if the variables ψ 1,m,n,j are not eliminated using Proposition 6, they can be used to further tighten the above constraints. For example, in the first constraint, p1´τ i,j,n q`p1´β 0,i,j q can be replaced with p1´ψ i,n,0,j q. This concludes the formulation of MINLP (A).

Relaxation and Solution Procedure
Apart from integrality requirements on stream pζ i,j q and heat exchanger variables (ρ i,j and χ i,j ), the remaining source of nonconvexity in the MINLP is the Underwood constraints. In this section, we describe the construction of a convex relaxation of Underwood constraints ((A18)-(A21)), referred to hereafter as the relaxation, defined using convex constraints that admits all feasible solutions. One of the challenges in constructing a valid relaxation is that the denominator of certain fractions in Underwood constraints can approach arbitrarily close to zero (see (A18)-(A21)). Consequently, off-the-shelf global solvers, such as BARON (Tawarmalani and Sahinidis, 2005), report an error and are not able to solve the problem. The common strategy used in the literature is to add/subtract θ (typically 10´2´10´3) from the bounds of θ ijq to prevent it from approaching either α q`1 or α q (see (A21)). However, this ad-hoc strategy has been adopted without a rigorous proof. Our numerical experiments suggest that the choice of this θ is not straightforward, and varies from one instance to another. In the following, we show that a rigorous relaxation for the fraction can be constructed although the denominator may approach close to zero.
In the following, we drop indices iklj. This is because, Underwood equations apply to a column, say Q iklj , and these indices are easily gleaned from the column specification or the associated split ri, ks{rl, js. Moreover, for notational convenience, we describe the relaxation using U " tpf, U, Υ, θq | (16); pf in p , f rs p , f ss p q P r0, F p s 3 , p " 1, 2; 0 ď p¨q ď p¨q up , @ p¨q P tU rs , U ss , Υ rs , Υ ss uu, where Here, we assume that column Q iklj performs the split of a binary mixture. Observe that (16a), the second inequality in (16b) and (16c) are simplified versions of (A18) and (A19) for binary mixtures. We ensure that all fractions are non-negative by factoring out a negative sign from the fractions whose denominator is negative (see (16)). Next, E rs and E ss denote lower bounds on nonlinear expressions in (16b) and (16c), respectively. We choose E rs (resp. E ss ) to be Υ rs (resp. Υ ss ) if the second inequality in (16b) (resp. (16c)) needs to be binding, as in (A20). Else, we choose the lower bound derived in (14) and (15). (16d), (16e), and (16f) correspond to (A21), (A22), and (A9), repectively. Lastly, we remark that, in (A19) and (A20), f rs 2 " f ss 1 " 0 for a split of a binary mixture. Since our purpose in restricting to the binary case is to illustrate the mathematical structure of relaxations, we do not consider this restriction. In general splits, one or more components may distribute between the distillate and residue.
The first step in standard approaches to relax U is to linearize Underwood constraints by introducing an auxiliary variable representing the graph of each fraction. Then, the restriction that this variable take the value of the fraction is replaced with the less stringent restriction that the variable lies in a convex set containing the graph of fraction. Instead, we reformulate U as described in §5.1 before linearizing the Underwood constraints.

Reformulation
We adapt classical Reformulation-Linearization Technique (RLT) (Sherali and Alameddine, 1992) to fractions, and reformulate U by appending RLT cuts derived using Underwood constraints. For clarity, we present the derivation of RLT cuts with Underwood minimum vapor constraint in the rectifying section (second inequality in (16b)), and describe the entire reformulated set towards the end. We multiply each Underwood constraint with the bound factors of θ, pθ´θ lo q, and pθ up´θ q. A naive approach would then disaggregate the product, leading to following which auxiliary variables are introduced to linearize each nonlinear term: H rs p " f rs p {|α p´θ |, Hθ rs p " f rs p θ{|α p´θ |, for p " 1, 2, and Υθ rs " Υ rs¨θ . Here, and in the rest of the article, the variables introduced to linearize a product will be written by underlining the concatenation of symbols, as in Υθ rs " Υ rs¨θ . Instead, we use polynomial long division prior to linearization, which transforms (17) to Next, we introduce auxiliary variables to linearize nonlinear terms: H rs p " f rs p {|α p´θ |, for p " 1, 2, and Υθ rs " Υ rs¨θ . We shall refer to the proposed variant as the Reformulation-Division-Linearization Technique (RDLT) of fractional terms, in order to easily distinguish and emphasize the use of polynomial division as an intermediate step. Clearly, RDLT cuts require fewer variables than those derived by naive application of RLT as described above. In addition, RDLT cuts lead to a tighter relaxation of U, which we demonstrate below.
The sets C and C 1 in Proposition 8 are typically created by relaxing the nonlinear expressions. We illustrate, via an example, that the relations in Proposition 8 can be strict. Example 1. Let, α 1 " 15, α 2 " 9, f lo 1 " f up 1 " 0.6, f lo 2 " f up 2 " 0.4, Υ lo "´10, Υ up " 10, θ lo " 9.1, θ up " 14.9. The sets C and C 1 are constructed by under-and over-estimating the nonlinear terms with their respective convex and concave envelopes. Figure 8(a) depicts the projection of sets S, S std , S RLT and S RDLT in Υ´θ space. It is clear that proj pf,Υ,θq pS std q Ą proj pf,Υ,θq pS RLT q Ą proj pf,Υ,θq pS RDLT q Ą S. Besides improving the quality of relaxation by introducing fewer auxiliary variables, RDLT has another benefit in our context that we describe next. θ´α2 is nonconvex (see Figure 8(b)), because it is a difference of two convex functions. When this function is multiplied by pθ´θ lo q (resp. (θ up´θ )), it becomes convex (resp. concave) (see 8(b)). In the naive RLT approach, where each fraction is relaxed independently, the product pf 1 {pα 1´θ qq¨pθ´θ lo q is disaggregated and relaxed as a difference of the convex envelope of f 1 θ{pα 1´θ q with the concave envelope of f 1 {pα 1´θ q. Whereas, the polynomial division step makes the convexity apparent revealing better ways to construct the relaxation.

Generalizations
We remark that RDLT can be used for problems with constraints that have the form ř r i"1 xigipyq hipyq ď x 0 , g i pyq ( r i"1 and h i pyq ( r i"1 are some polynomials of y. We follow the steps below to derive RDLT cuts.
1. We multiply the constraint by some ratio of polynomials of y, npyq{dpyq, such that the sign of the ratio does not change over the domain of y. Here, we assume, w.l.o.g, that npyq{dpyq ě 0 over the domain of y.
2. We use polynomial long division to express each gipyq¨npyq hipyq¨dpyq " m i pyq`k ipyq lipyq such that degpk i q ă degpl i q, where degpk i q denotes degree of polynomial k i pyq. 3. We factorize l i pyq and express it as a product of polynomials tq ij pyqu si j"1 that are non-factorizable over real numbers (e.g., y`2 or y 2`y`1 ). 4. We use the general theorem of partial fraction decomposition to express each fraction k i pyq{l i pyq as ř si j"1 p ij pyq{q ij pyq, where degpp ij q ă degpq ij q. This transforms the constraint to ř r i"1´x i¨mi pyq`ř si j"1 x i¨pij pyq{q ij pyq¯ď x 0¨n pyq{dpyq. 5. We linearize the constraint by introducing auxiliary variables for each nonlinear term.
The reformulation described earlier is a specific case, where we chose to multiply each Underwood constraint by pθ´θ lo q and pθ up´θ q. By changing the factor used in the reformulation step, we can derive alternative RDLT cuts by following the steps described above. As an illustration, we derive two types of additional RDLT cuts for reformulation of U. While we do not use these cuts for our extensive computational experiments, we demonstrate with numerical examples in §6 that they further improve the relaxation for some instances.

RDLT cuts with quadratic polynomials:
Here, we choose the product of bound factors of θ, viz. pθ´θ lo q 2 , pθ´θ lo q¨pθ up´θ q and pθ up´θ q 2 , for reformulation. As an illustration, we derive the RDLT cut by multiplying the second inequality in (16b) with pθ´θ lo q¨pθ up´θ q. The remaining RDLT cuts are derived in a similar fashion. Steps 1 and 2 lead to 2 ÿ p"1˜α p f rs p¨p θ`α p´θ lo´θup q´α p pα p´θ lo qpα p´θ up qf rs p α p´θ¸ď Υ rs¨p θ´θ lo q¨pθ up´θ q.
Since (24) is already in the form attained in Step 4, we do not need Steps 3 and 4. Finally, we disaggregate the products of f rs p and Υ rs with polynomials of θ, and linearize (24) by introducing auxiliary variables for f rs p {pα p´θ q, f rs p¨θ , Υ rs¨θ2 and Υ¨θ. RDLT cuts with inverse bound factors: Here, we use inverse bound factors`1 θ´1 θ up˘a nd`1 θ lo´1 θ˘f or reformulation. Since`1 θ´1 θ up˘" θ up´θ θ up¨θ , inverse bound factors are essentially ratios of first-degree polynomial to another first-degree polynomial. As before, for illustration, we derive the RDLT cut obtained by multiplying the second inequality in (16b) with`1 θ´1 θ up˘. The remaining RDLT cuts are obtained in a similar fashion.
Step 1 leads to ř 2 p"1 αpf rs p pαp´θqθ´α p f rs p pαp´θqθ up ď Υ rs θ´Υ rs θ up , which is already in the form described in Step 2. Further, the denominator of each fraction is already expressed as product of non-factorizable polynomials. Next, we use partial fraction decomposition (Step 4) to obtain The nonconvexity in U ref is due to F 1 , F 2 , and V. We convexify these sets to construct a convex relaxation of U ref . However, we first assume that α 2 ă θ lo and θ up ă α 1 , and relax this assumption later in §5.3. This assumption prevents the denominator of fractions in F 1 and F 2 from becoming zero. This discussion is needed for two reasons: (i) it will guide us in deriving additional valid cuts needed to strengthen the relaxation when θ " α 2 and θ " α 1 are admissible (ii) it is needed to construct a piecewise relaxation in §5.4, where we discretize the domain of θ such that every partition excluding the extreme partitions satisfy α 2 ă θ lo ď θ ď θ up ă α 1 .
The standard approach to create a relaxation is to replace each equality H p " f p¨Tp pθq in F p (resp. Υθ " Υ¨θ in V) with a less stringent restriction that H p (resp. Υθ) lies in the convex hull of f p¨Tp pθq (resp. Υ¨θ) over a rectangle defined by the ranges of f p (resp. Υ) and θ. However, this approach does not take advantage of the fact that the component (resp. vapor) flows are constrained by mass balances (see (16e),(16f)) and, thus, results in a weaker relaxation. Instead, we use Proposition 9, which describes the construction of simultaneous hull of multiple nonlinear terms over a polytope (not necessarily a hyperrectangle), to construct a tighter relaxation of U ref .
Proof. Since S is compact, its convex hull is compact and, by Krein-Milman theorem, is the convex hull of its extreme points. Therefore, we determine the extreme points of S, and take their convex hull to obtain ConvpSq. When y is restricted to y P ry lo , y up s, the set S " tpx, y, z, xyq | z " gpyq x, xy " y x, x P X, y " yu can be expressed as an affine transform of X. Thus, the extreme points of S project to the set of extreme points of X and we may restrict attention to these points in order to construct ConvpSq. Let S i , for i " 1, . . . , m, denote the set S where x is restricted to v i i.e., S i " tpx, y, z, xyq | z " v i gpyq, xy " v i y, x " v i , y P ry lo , y up su. Then, ConvpSq is given as the convex hull of disjunctive union of S i , i " 1, . . . , m, i.e., ConvpSq " ConvpS 1 Y¨¨¨Y S m q " ConvpConvpS 1 q Y¨¨¨Y ConvpS m qq.
To determine ConvpS i q, we reformulate each S i as S i " tpx, y, z, xy, wq | z " v i w, xy " v i y, w " gpyq, x " v i , y lo ď y ď y up u, which is an affine transform of the set tpy, wq P ry lo , y up sˆR | w " gpyqu. This implies that it suffices to convexify the latter set to obtain ConvpS i q " proj π tpπ, wq | (28)u, where w ě gpyq, w ď gpy lo q`ˆg py up q´gpy lo q y up´ylo˙p y´y lo q, The disjunctive union of ConvpS i q, i " 1, . . . , m, leads to (26), where w i and y i are to be regarded as linearization of λ i w and λ i y, respectively.
Remark 7. In Proposition 9, if ConvpS i q (see proof for definition) is bounded, closed and cone-quadratic representable (CQR), for i " 1, . . . , m, then ConvpSq is CQR (see Proposition 3.3.5 in Ben-Tal and Nemirovski (2001)). This result also applies to other conic representations. Let P δ,1´δ ě |x 3 |u where 0 ă δ ă 1 is the power-cone, and K exp " tx 1 ě x 2¨e xppx 3 {x 2 q, x 2 ą 0u Y tpx 1 , 0, x 3 q | x 1 ě 0, x 3 ď 0u is the exponential-cone. It is known that various elementary functions have cone representations (MOSEK, 2020). For example, let gpyq " |y| δ where δ ą 1 (resp. gpyq " y δ where δ ă 0). Then, w i ě g˚pλ i , y i q in Proposition 9 can be replaced with pw i , λ i , y i q P P 1{δ,1´1{δ 3 (resp. pw i , y i , λ i q P P 1{p1´δq,´δ{p1´δq 3 ). For this work, we are interested in δ "´1 and δ " 2 (for reformulation with quadratic polynomials described in §5.1.1). Next, let gpyq "´lnpyq, y ą 0 (resp. gpyq " exppyq), which arises in formulations for identifying thermodynamically efficient distillation configurations (see Jiang et al. (2019a)). Here, we replace w i ě g˚pλ i , y i q in Proposition 9 with py i , λ i ,´w i q P K exp (resp. pw i , λ i , y i q P K exp ). Remark 8. In Proposition 9, when gpyq is nonlinear, the convex hull description has nonlinear constraints (see (26a)). To capitalize on LP solvers, we derive a polyhedral outer-approximation of ConvpSq by outerapproximating the convex hull of each S i before taking their disjunctive union. Let y r P ry lo , y up s for r " 1, . . . , R. Then, an outer-approximation of the convex hull of S i is given by Conv OA pS i q " proj π tpπ, wq | w ě maxtgpy r qg 1 py r qpy´y r qu R r"1 ; (28b)´(28d)u, where g 1 pyq denotes the first derivative of gpyq. The disjunctive union of Conv OA pS i q, i " 1, . . . , m, yields an outer-approximation of the convex hull of S, given by Conv OA pSq " proj π tpπ, y 1 , . . . , y m , w 1 , . . . , w m , w, λ 1 , . . . , λ m q | w i ě maxtgpy r qλ i`g1 py r qpy i´yr λ i qu R r"1 ; (26b)´(26e)u. Now, consider the set F p . We lift F p to a higher dimensional space by appending bilinear terms of the form f p¨θ i.e., F p " tpf, θ, H, f θq | (22), f θ in p " f in p¨θ , f θ rs p " f rs p¨θ , f θ ss p " f ss p¨θ u. Observe that the fractions and bilinear terms in F p are defined over the polytope obtained by the intersection of hyperplane f in p " f rs p`f ss p with the hypercube r0, F p s 3 (see (22)). We now use Proposition 9 to obtain ConvpF p q " tpf p , θ, H p , f θ p q | (29)u (see §E for a detailed derivation), where H rs p ď f rs p T p pθ lo q`" H ss p ď f ss p T p pθ lo q`" and the positively homogeneous function Tp pλ, θq is defined as in (27) from T p pθq. Note that the convex hull description does not require introduction of auxiliary variables. This is in contrast to the typical application of disjunctive programming, where new variables are introduced to derive the convex hull in a lifted space. We remark that the above yields a tighter relaxation of F p compared to the one obtained by relaxing each fraction and bilinear term separately over the bounds of f in p , f rs p , f ss p , and θ. This is because the first two equations in (29f) are not implied in the latter set. Although, these relations can be obtained using RLT, appending these constraints does not result in (29). This is because, the set described in (29) is the simultaneous convex hull of the fraction and bilinear terms. It is known that the simultaneous hull of these functions is strictly contained in the intersection of their individual hulls (see Example 3.8 in Tawarmalani (2010)). In particular, (29b) and (29c), which are linearizations of´f rs p |αp´θ|¨p θ up´θ q¨pθ´θ lo q ď 0 and´f ss p |αp´θ|¨p θ up´θ q¨pθ´θ lo q ď 0 respectively, are not implied in the intersection of individual convex hulls.
The convex hull description in (29) is cone-quadratic representable (see Remark 7), since the constraints in (29a) can be expressed as second-order cones. For example, H rs H 1¨b pα 1 f rs 1´f θ rs 1 q ě |f rs 1 | (Note that 0 ď θ lo f rs 1´f θ rs 1 ă α 1 f rs 1´f θ rs 1 ). However, we use the conequadratic representation only in §6. For our computational experiments in §7, we use its outer-approximation given by Conv OA pF p q " tpf p , θ, H p , f θ p q | H rs p ě maxtf rs p T p pθ r q`T 1 p pθ r qpf θ rs p´θ r f rs p qu R r"1 , H ss p ě maxtf ss p T p pθ r q`T 1 p pθ r qpf θ ss p´θ r f ss p qu R r"1 , (29b)´(29f)u for some θ r P rθ lo , θ up s, r " 1, . . . , R, where T 1 p pθ r q denotes the derivative of T p pθq at θ r ; see Remark 8.
Next, consider the set V, which contains bilinear terms defined over a polytope obtained by the intersection of a hyperrectangle in the positive orthant with the hyperplane U rs´U ss " Υ rs´Υss . Clearly, Proposition 9 can be used to construct the convex hull of V (only (26c), the equation with xy as the left-hand-side in (26d), and (26e) are needed to construct the hull). However, Proposition 9 requires enumeration of the extreme points of X. Instead, in this context, it is more convenient to directly use Proposition 10, which is a special case of Proposition 2.2 in Davarnia et al. (2017), to obtain ConvpVq " tpU, Υ, θ, U θ, Υθq | (30)u, where pU θ rs´θlo U rs q´pU θ ss´θlo U ss q " pΥθ rs´θlo Υ rs q´pΥθ ss´θlo Υ ss q, pθ up U rs´U θ rs q´pθ up U ss´U θ ss q " pθ up Υ rs´Υ θ rs q´pθ up Υ ss´Υ θ ss q, (30b) 0 ď p¨qθ´θ lo p¨q ď p¨q up pθ´θ lo q, @ p¨q P tU rs , U ss , Υ rs , Υ ss u, (30c) 0 ď θ up p¨q´p¨qθ ď p¨q up pθ up´θ q, @ p¨q P tU rs , U ss , Υ rs , Υ ss u.
Proof. See §F in the Appendix.
Remark 9. We remark that ConvpF p q and ConvpVq in (29) and (30) imply the convex envelope of ř 2 p"1 "´α p pαp´θ up qf rs p αp´θ`α p f rs p ı`Υ rs¨θ´Υrs¨θup over bound constraints on f rs p , θ, and Υ rs (see (18b)). This is because, when all f rs p and Υ rs are fixed, the function is concave in θ. Then, by Theorem 1.4 in Rikun (1997), it follows that the convex envelope is obtained by replacing´α p pαp´θ up qf rs p αp´θ for all p and Υ rs¨θ by their convex envelopes.
We comment on the construction of convex relaxations of U when additional RDLT cuts described in §5.1.1 are appended to U ref . Reformulation of Underwood constraints using quadratic polynomials of θ introduces nonconvex terms of the form f p¨θ , Υ¨θ 2 (see (24)), in addition to the existing f p¨Tp pθq and Υ¨θ terms in U ref . We relax f p¨Tp pθq and f p¨θ using the simultaneous hull description in (29). Although Proposition 9 yields the simultaneous hull of Υ¨θ 2 and Υ¨θ terms over the polytope in V, we do not implement this relaxation. This is because the hull description does not project onto the space of problem variables in a striaghtforward manner. Instead, we convexify each pair of Υ¨θ 2 and Υ¨θ terms over a box using Corollary 2, and append the RLT cuts U rs¨θ2´U ss¨θ2 " Υ rs¨θ2´Υss¨θ2 and U rs¨θ´U ss¨θ " Υ rs¨θ´Υss¨θ .
On the other hand, reformulation of Underwood constraints using inverse bound factors introduces nonconvex terms of the form f p¨θ´1 and Υ¨θ´1 (see (25)), in addition to the existing f p¨Tp pθq and Υθ terms in U ref . We relax f p¨Tp pθq and f p¨θ using the simultaneous hull description in (29). We use a similar hull description, obtained using Proposition 9, to relax f p¨θ´1 and f p¨θ . Finally, for the same reason mentioned above, we convexify each pair of Υ¨θ´1 and Υ¨θ terms using Corolloary 2, and append RLT cuts U rs¨θ´1´U ss¨θ´1 " Υ rs¨θ´1´Υss¨θ´1 and U rs¨θ´U ss¨θ " Υ rs¨θ´Υss¨θ . Corollary 2. Let B " rx lo , x up sˆry lo , y up sˆR 2 , where we assume 0 ď x lo , gpyq : ry lo , y up s Ñ R is convex, and S " tpx, y, z, xyq P B | xy " x¨y, z " x¨gpyqu. Then, ConvpSq " tpx, y, z, xyq | (31)u, where z ď gpy lo q¨x`" gpy up q´gpy lo q y up´ylo ´x y´y lo x¯, px up´x qy lo ď x up y´xy ď px up´x qy up , px´x lo qy lo ď xy´x lo y ď px´x lo qy up , and g˚pλ˚, y˚q is defined as in (27). Further, the outer-approximation of the convex hull is Conv OA pSq " tpx, y, z, xyq | (32), (31b)´(31c)u, where z ě x lo x up´xlo maxtgpy r qpx up´x q`g 1 py r qpx up y´xy´px up´x qy r qu R r"1x up x up´xlo maxtgpy r qpx´x lo q`g 1 py r qpxy´x lo y´px´x lo qy r qu R r"1 , for some y r P ry lo , y up s, r " 1, . . . , R, and g 1 pyq denotes the first derivative of gpyq w.r.t y.
Proof. See §G in the Appendix.

Valid
Relaxation for θ lo " α 2 and/or θ up " α 1 In the previous subsection, we have assumed that α 2 ă θ lo and θ up ă α 1 . Instead, if α 1 and/or α 2 is an admissible value of θ, we cannot directly use (29) to convexify F p , because T 1 pα 1 q and T 2 pα 2 q are not well-defined. To construct a valid relaxation, we first restrict the admissible values of θ to a subset of the interval rα 2 , α 1 s by recognizing that each fraction in U is bounded. pΥ rs q up pα 1´α2 q`α 1 F 1`α2 F 2 α 1 pα 1´α2 q , (33a) pH rs 2 q up "

Proposition 11. (i) Valid upper bounds on
(ii) the admissible region of θ in the interval rα 2 , α 1 s is given by Proof. (i) Consider the second inequality in (21c). Since this inequality holds for any θ lo less than θ, if we substitute θ lo with α 2 , the inequality remains valid. Then, we obtain (33a) from α 1 pα 1´α2 qH rs 1 ď pΥθ rs´α 2 Υ rs q`α 1 f rs 1`α 2 f rs 2 ď pΥ rs q up pα 1´α2 q`α 1 F 1`α2 F 2 , where the last inequality is because f rs 1 ď F 1 , f rs 2 ď F 2 , and Υθ rs´α 2 Υ rs ď pΥ rs q up pα 1´α2 q. Similarly, we substitute θ up " α 1 in the first inequality in (21d), and rearrange to get α 2 pα 1´α2 qH rs 2 ď´E rs¨p α 1´θ q`α 1 f rs 1`α f rs 2 . We maximize the right hand side by substituting f rs 1 " F 1 , f rs 2 " F 2 , and pE rs¨p α 1´θ qq by its lower bound which is computed using the bounds on E rs and θ. This leads to the bound in (33b).
We remark that the bounds on H in p and H ss p for p " 1, 2 can be computed in the same manner as in the proof of (i) in Proposition 11. Even when additional fractions are present in the Underwood constraints, each fraction can be bounded, since the remaining fractions are strictly bounded in the interval of θ. We revisit the A PREPRINT -OCTOBER 26, 2020 argument on bounds of θ in light of Proposition 11. As mentioned before, the common approach used in the literature to overcome the singularity arising due to θ approaching one of the adjoining relative volatilities has been to restrict θ to belong to rα 2` θ , α 1´ θ s. However, observe that our bounds in (34) depend on f rs 1 and f rs 2 . This explains the difficulty we encountered in choosing a value for θ in our computations with prior formulations. We have found that there are instances when θ is fairly close to one of the relative volatilities, particularly when the corresponding flow is small. We will provide a rigorous approach to addressing this singularity using (34). Our approach will be to construct a relaxation of F 1 as the intersection of simultaneous convex hulls of f 1¨T1 pθq and f 1¨θ . For brevity, we only discuss the relaxation for F 1 in detail, and remark that a similar result is easily derived for F 2 . Proposition 12. Let H 1 " tpf 1 , θ, H Proof. See §H in the Appendix.
The convex hull in Proposition 12 requires several additional variables. To avoid the introduction of these additional variables, we use its relaxation, H 1,Relax , derived in §I and shown below: where θ r P rθ lo , α 1 q, r " 1, . . . , R. Here, we argue from first principles that (36) is a valid relaxation. To derive the first inequality in (36a), observe that H 1 ě f 1¨T1 pθq ě f 1¨m axtT 1 pθ r q`T 1 1 pθ r qpθ´θ r qu R r"1 . Disaggregating the product and linearizing the bilinear term yields (36a). To derive the second inequality in (36a), we begin with H 1¨p α 1´θ q ď f 1 , and replace the bilinear term on the left hand side with its convex envelope. (36b) is the convex hull of f θ 1 " f 1¨θ over r0, F 1 sˆrθ lo , α 1 s, and (36c) is the same as (34). Using (36), we obtain a valid relaxation of F 1 given by F 1,Relax " tpf 1 , θ, H 1 , f θ 1 q | pf in 1 , θ, H in 1 , f θ in 1 q P H in 1,Relax , pf rs 1 , θ, H rs 1 , f θ rs 1 q P H rs 1,Relax , pf rs 1 , θ, H rs 1 , f θ rs 1 q P H ss 1,Relax , H in 1 " H rs (29), the last two equations in the relaxation are derived by multiplying the component mass balance, (16f), with T 1 pθq and θ, respectively.

Discretization and Solution Procedure
In this work, instead of using convex relaxations of U in a spatial branch-and-bound framework to solve the MINLP, we construct a piecewise relaxation (see Definition 1) that is iteratively improved until we prove r´o ptimality. This approach capitalizes on state-of-the-art MIP solvers, such as Gurobi. Definition 1 (Piecewise Relaxation). Let x " px 1 , . . . , x n q, B " rx lo , x up sˆry lo , y up s Ă R n`1 , S " tpx, yq P B | g i px, yq ď 0, i " 1, . . . , mu, and S Relax " tpx, yq P B | q g i px, yq ď 0, i " 1, . . . , mu be its convex relaxation , where tq g i u m i"1 denote convex underestimators of tg i u m i"1 over B. Let, the domain of y be partitioned as I " trY 0 , Y 1 s, . . . , rY |I|´1 , Y |I| su with Y 0 " y lo , Y |I| " y up and Y 0 ď Y 1 ď . . . Y |I| . By piecewise relaxation of S, we refer to Ť |I| t"1 S t,Relax , where S t,Relax " tpx, y, zq P B t | q g i,t px, yq ď 0, i " 1, . . . , mu, B t " rx lo , x up sˆrY t´1 , Y t s, and q g i,t is the convex under-estimator of g i over B t .
A PREPRINT -OCTOBER 26, 2020 Piecewise relaxation of U can be constructed by partitioning the domain of Underwood root as I " trΘ 0 , Θ 1 s, . . . , rΘ |I´1| , Θ |I| su, where Θ 0 " α 2 , Θ |I| " α 1 , and Θ 0 ď Θ 1 ď¨¨¨ď Θ |I| , and taking the union of sets Ť |I| t"1 U t,Relax , where U t,Relax denotes the convex relaxation of U restricted to θ P rΘ t´1 , Θ t s. The set U t,Relax is constructed as outlined in §5.2 and §5.3. Next, using standard disjunctive programming techniques, the piecewise relaxation can be expressed as a Mixed Integer Program (MIP). While this approach leads to a locally ideal formulation, it leads to a bigger problem size, because of which the computational time required is higher. Thus, in favor of smaller problem size, we do the following.
Instead of reformulating U in each partition using the local bound factors of θ, we reformulate with the overall bound factors of θ: pθ´α 2 q and pα 1´θ q. Next, we require that pf p , θ, H p , f θ p q, p " 1, 2, and pU, Υ, θ, U θ, Υθq lie in piecewise relaxations of F p and V, respectively. We choose piecewise relaxation of F 1 to be Ť |I|´1 t"1 Conv OA pF 1,t q Y F 1,|I|,Relax , piecewise relaxation of F 2 to be F 2,1,Relax Y Ť |I| t"2 Conv OA pF 2,t q, and piecewise relaxation of V to be Ť |I| t"1 ConvpV t q. Here, the additional subscript t denotes that the set is restricted to θ P rΘ t´1 , Θ t s. Observe that if zero is not an admissible value to the denominators of the fractions, we use outer-approximation of convex hulls derived in §5.2 to relax F p . Otherwise, we use a relaxation of the convex hull description, such as the one derived in §5.3. We use disjunctive programming to express the piecewise relaxations as mixed-integer sets (see §J for description of the sets). In §6, we illustrate through numerical examples the impact of various aspects described in this section in strengthening the overall relaxation of MINLP (A).
Algorithm 1 outlines our approach to solve the MINLP. We start with a coarse discretization and use an adaptive partitioning scheme to iteratively refine the partitions until r´o ptimality is achieved. To avoid numerical issues, we maintain that each partition, pΘ t ijq´Θ t´1 ijq q, is at least MinPrtSize in length.

Effect of Individual Cuts on Relaxation
This section illustrates, through numerical examples, the impact of of various aspects described in §5 in strengthening the overall relaxation of MINLP (A). We highlight the individual effect of RDLT cuts derived from Underwood constraints, simultaneous hulls derived in §5.2, and discretization on the overall relaxation. In all the scenarios below, stream and heat exchanger variables are considered to be binary.
Scenario 1 : (BARON's root node relaxation) Here, we use BARON 18.5.8, on GAMS 25.1, to construct and solve the relaxation of MINLP (A). This is achieved by specifying BARON option MaxIter = 1, which terminates the branch-and-cut algorithm after processing the root node. We let θ ijq P rα q`1` θ , α q´ θ s, with θ " 10´7, for every q j´1 i , ri, js P S to avoid a possible division by zero. We use BARON's root node relaxation as a reference for comparison. We remark that BARON solves MIP relaxations as needed (Kılınç and Sahinidis, 2018). We also verified that the bound obtained is close to solving a factorable MIP relaxation.
Scenario 2 : (Simultaneous hull of fractional terms) This scenario illustrates the improvement in relaxation due to the use of simultaneous convexification techniques. We linearize all Underwood constraints in the MINLP by introducing auxiliary variables for each fraction. To relax fractional terms, we use (29), or (36) if zero is an admissible value for the range of the denominator of fractions. The nonlinear constraints in (29) are expressed as second-order cones, and the resulting Mixed Integer Second-order Cone Program (MISOCP) is solved with Gurobi 8.0 using Gurobi/MATLAB interface.
Scenario 4 : (RDLT with quadratic polynomials of θ) To the relaxation in Scenario 3, we add cuts derived by reformulating Underwood constraints with quadratic polynomials of θ (see (24)), as described in §5.1.1. This introduces additional nonlinear terms of the form Υ¨θ 2 , which we relax in the manner described towards the end of §5.2.
Scenario 5 : (RDLT with inverse bound factors of θ) To the relaxation in Scenario 3, we add cuts derived by reformulating Underwood constraints with inverse bound factors (see (25)). This introduces additional nonlinear terms of the form f 1 {θ and V {θ, which we relax in the manner described towards the end of §5.2.
The relative tolerance r is achieved.
end 32 return I ijq 33 end Algorithm 1: Algorithm to solve MINLP using an adaptive discretization scheme Scenario 6 : (Discretization) Finally, to illustrate the potential of discretization, we construct piecewise relaxation of Scenario 3. We discretize the domain of each Underwood root into two partitions, and choose the roots of columns performing the split of the process feed, tθ 1N q u N´1 q"1 , as the partition points. In other words, we let I ijq " trα q`1 , θ 1N q s, rθ 1N q , α q su for i ď q ă j and ri, js P S. As pointed out in Remark 5, these roots can be computed prior to solving the optimization problem. We construct the piecewise relaxation of MINLP (A) as outlined in §5.4. Table 3 reports the percentage gap value, defined as % Gap " 100ˆˆ1´O ptimal value of relaxation Optimal value of (A)˙( 37) on a set of cases evaluated for all the Scenarios. To compare against BARON, we also report % gap closed (numbers in parenthesis in Table 3), defined as % Gap Closed " 100ˆˆ1´O ptimal value of (A)´Optimal value of relaxation Optimal value of (A)´Optimal value in Scenario 1˙ (  38) We refer to a particular combination of parameter settings: N , tF p u N p"1 , tα p u N p"1 , Φ 1,N and tΦ p,p u N p"1 , as a case. The parameter settings for the cases considered in Table 3 are listed in the caption. It is worth noting that Case-A (Caballero and Grossmann, 2004), Case-B and Case-C (Nadgir and Liu, 1983) correspond to physical mixtures: mixture of alcohols, mixture of light paraffins and mixture of light olefins and paraffins. The remaining cases do not directly correspond to physical mixtures, but are representative of specific classes of separations (see Giridhar and Agrawal (2010a) for more details). Under Scenario 2, we report % Gap value, and % Gap closed for all cases when simultaneous hulls are used to convexify fractions. It can be observed that, this approach closes on an average 45.8% of the gap. In particular, in Case-E, implementation of simultaneous hull completely closes the gap at root node. Next, under Scenario 3, we report the combined effect of our RDLT approach and simultaneous hulls. This approach closes on an average 74.1% of the gap. Under Scenarios 4 and 5, we report further improvement in relaxation due to addition RDLT cuts discussed in §5.1.1 to the relaxation in Scenario 3. RDLT cuts with quadratic polynomials of Underwood roots closes the gap completely in Case-B. Finally, the gap can be completely closed for all the cases considered in Table 3 by discretizing the domain of Underwood root into two partitions, as described in Scenario 6.

Computational Results
We conducted computational experiments on a test set of 496 cases, taken from (Giridhar and Agrawal, 2010a;Nallasivam et al., 2013), which is a representative of a majority of separations. Parameter settings for the test set are listed in §M in e-companion. In this section, we demonstrate that our proposed approach is able to solve MINLP (A) within a relative tolerance of 1%. We also compare the performance of our approach with prior approaches in the literature (Caballero and Grossmann, 2004;Nallasivam et al., 2016;Tumbalam Gooty et al., 2019). Since the prior approaches develop an (MI)NLP model, we use BARON 18.5.8 via GAMS 25.1 to solve these (MI)NLPs, where all BARON options are set at their default values. For the adpative partitioning scheme described in Algorithm 1, we use Gurobi 8.0 (Gurobi Optimization, 2018) to solve the resulting MIPs, and use IPOPT (Wächter and Biegler, 2006) as a local solver. The model is loaded into Gurobi using the MATLAB/Gurobi interface, while IPOPT is used via MATLAB/GAMS interface and GAMS 25.1. We used single CPU thread to solve the MIPs so as to keep the comparison with BARON fair. Besides the setting of number of threads, the remaining options for Gurobi and IPOPT were left at their defaults. All computations were done on a Dell Optiplex 5040 with 16 GB RAM, which has Intel Core i7-6700 3.4 GHz processor and is running 64-bit Windows 7.

Comparison with Prior Approaches
Here, we compare the performance of three approaches, namely those of Caballero and Grossmann (2006);Tumbalam Gooty et al. (2019), and the one proposed here. For all the computations, we set the relative tolerance for convergence ( r ), defined as where BestLB and BestUB are the best-known relaxation bound and feasible solution, to 1% i.e., r " 0.01. We impose a CPU time limit of five hours as the termination criterion.

Case-B and
Case-C from Nadgir and Liu (1983), and Case-D through Case-H from Giridhar and Agrawal (2010a).
A PREPRINT -OCTOBER 26, 2020 Approach 1 : We solve MINLP (A) using the adaptive partitioning approach described in Algorithm 1. We begin with four partitions for each Underwood root i.e., I ijq " trα q`1 , pα q`1`θ1N q q{2s, rpα q`1θ 1N q q{2, θ 1N q s, rθ 1N q , pα q`θ1N q q{2s, rpθ 1N q`αq q{2, α q su for every q j´1 i , ri, js P S. We compute the Underwood roots for the splits of the process feed tθ 1N q u N´1 q"1 prior to solving the MINLP (see Remark 5). For all but 4 cases, we set MinPrtSize " 10´3. For the remaining cases, we reduced MinPrtSize to 10´4 in order to achieve the relative tolerance of 1%. Finally, we point out that the upper bounds on material flows are computed by solving (13) and φ " 1.5. We note that (40) is the objective function value corresponding to a feasible point of one of the admissible configurations, commonly known in literature as Fully Thermally Coupled or Petlyuk configuration (see Fidkowski and Krolikowski (1986); Halvorsen and Skogestad (2003b)). Approach 2 : We obtained the GAMS code of the model proposed in Caballero and Grossmann (2006) from the MINLP library (Caballero and Grossmann, 2009). There, the authors were interested in identifying the configuration minimizing the total annual cost. For our computations, we modify their code in the following manner. First, as mentioned in (Tumbalam Gooty et al., 2019), the model of Caballero and Grossmann (2009) admits solutions that are physically infeasible. This is because the constraints corresponding to (A19) in their model should be tight for certain Underwood roots, and their model does not impose this requirement. We have added these missing constraints to their GAMS code. Second, the authors employed the BigM approach in order to transform certain disjunctions into a set of inequalities. Unfortunately, the BigM value used for vapor and liquid bypass in their GAMS code made a few test cases infeasible. Therefore, we specified 2.5VD˚as the BigM value for the vapor and liquid bypasses. This number was found by choosing the smallest BigM value for which we found a feasible solution. Third, the authors use a parameter θ and restrict θ ijq P rα q`1` θ , α q´ θ s for q j´1 i , ri, js P P in order to avoid the singularity associated with θ ijq approaching α q or α q`1 . Their choice of θ , in some cases, made the optimal solution infeasible. Empirically, we found that θ " 10´4 does not cut off the optimal solution, so we set θ " 10´4. Fourth, the cost equations required for the evaluation of the objective function were removed from the model, and the objective function was modified to compute the total vapor duty instead. The resulting MINLP is then solved with BARON. Approach 3 : Here, we consider the MINLP proposed in Tumbalam Gooty et al. (2019). For a consistent comparison, we set the upper bound on all vapor flows to be 1.5VD˚. Further, we restrict θ ijq P rα q`1` θ , α q´ θ s, where θ " 10´4, for q j´1 i , ri, js P S in order to avoid the singularity associated with θ ijq approaching α q or α q`1 . The resulting MINLP is then solved using BARON.
Figure 9(a) shows the percentage of cases solved to 1%-optimality against time, with Approach 1 (solid blue curve), Approach 2 (dotted black curve), and Approach 3 (dashed red curve). Observe that Approach 2 solves about 10% of cases to 1%-optimality within five hours. This is not surprising because Grossmann (2004, 2006) also reported difficulties in convergence. To overcome the challenges, the authors architected an algorithm by modifying logic-based outer-approximation. While the method resulted in good solutions, optimality was not guaranteed. Approach 3 solves 64% of the cases in the test set.
We remark that Tumbalam Gooty et al. (2019) introduced a new search-space formulation, derived cuts that exploit monotonicity of Underwood constraints, and modeled the absence/presence of a column using disjunctions. Nevertheless, this approach fails to solve the problem to 1%-optimality for 36% of the cases. The progress of lower bound for a majority of these cases is either stagnant or very slow. Figure 9(b) depicts the cumulative percentage of cases as a function of the remaining duality gap at the end of five hours. In contrast, our approach, for the first time, solves all 496 cases from this test set within an optimality tolerance of 1%. Figure 10 depicts cumulative percentage of cases as a function of the remaining duality gap at specific time instances for Approach 1. This graph demonstrates that our solution approach, with a CPU time of twenty minutes, already outperforms the best prior MINLP based approach allowed to run for a CPU time of five hours. Further, within 1800 s (green curve), 3600s (magenta curve) and 7200s (black curve), the proposed approach solves all 496 cases to less than 5.5%, 3.5% and 2.5% gap, respectively. Since (A) is primarily designed as a screening tool for an otherwise highly cumbersome search of optimal distillation configuration, practicing engineers can use Approach 1 to quickly identify near optimal solutions that are worthy of further exploration. Although we do not provide specific configurations found using our procedure, the potential

Comparison with Nallasivam et al. (2016)
Recently, Nallasivam et al. (2016) proposed an alternative technique that relies on explicit enumeration for identifying distillation configuration requiring the least vapor duty. After enumerating all the configurations, an NLP is formulated for each configuration and solved to 1%´optimality with BARON. We refer to this as Approach 4. We compare the performance of Approach 4, with Approaches 1 and 3 by fixing the discrete decisions to a specific configuration. We choose Fully Thermally Coupled (FTC) configuration, characterized by ζ i,j " 1 @ ri, js P T , χ i,j " 0 @ pi, jq P Cztp1, 1qu, χ 1,1 " 1, ρ i,j " 0 @ pi, jq P RztpN, N qu, and ρ N,N " 1, for comparison. This comparison ignores the advances in the search space formulation discussed in §4.2 and other advances that relate Underwood constraints with stream variables, since we fix the binary variables a priori. We set the time limit as one hour and a relative gap of 1% ( r " 0.01) as termination criteria.
Time (s)  Figure 11 depicts the percentage of cases solved as a function of computational time for the three approaches.
Clearly, BARON solves more number of cases to 1%´optimality with Approach 3 than with Approach 4. Despite the improvement, only 82% of the cases are solved to 1%´optimality using Approach 3. In contrast, our approach solves all cases in this test set within 100 s.

Concluding Remarks
This work addressed the optimal design of distillation configurations, which are widely used in all chemical and petrochemical industries, and are significant consumers of energy in the world economy. We proposed a novel MINLP that identifies energy-efficient configurations for a given application. Given the complexity from combinatorial explosion of the choice set and nonconvex Underwood constraints, this problem has resisted solution approaches. In this paper, we report on the first successful approach and solve this problem to global optimality for five-component mixtures. The key contributions that make this possible are (i) new formulation for discrete choices that is strictly tighter than the previous formulations, (ii) new valid cuts to the problem using RDLT, and various other convexification results for special structures, and (iii) discretization techniques and an adaptive partitioning scheme to solve the MINLP to ´optimality. On a test set that is a representative of a majority of five-component separations, we demonstrated that our approach solves all the instances in a reasonable amount of time, which was not possible using existing approaches. In summary, this paper describes the first solution approach that can reliably and quickly screen several thousands of alternative distillation configurations and identify solutions that consume less energy and, thereby, lead to less greenhouse gas emissions. This approach has the potential to reduce the carbon footprint and energy usage of thermal separation processes.

Disclaimer
The information, data, or work presented herein was funded in part by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

B Proof of Remark 3
We show the proof for τ i,k,j variables, and the proof for β i,l,j variables is similar. By Remark 2, the convex hull of ν i,k,j " ś j n"k p1´ζ i,n q over pζ i,k , . . . , ζ i,j q P r0, 1s j´k`1 , given by is implied from (A3), for every ri, js P P, k j´1 i . We use the above inequalities, in addition to (A2) and (A3), for the proof. We consider two cases: k`1 ă j and k`1 " j. When, k`1 ă j, the convex hull of τ i,k,j " ζ i,k p1´ζ i,k`1 q . . . p1´ζ i,j´1 qζ i,j over pζ i,k , . . . , ζ i,j q P r0, 1s j´k`1 is given by (Crama, 1993) On the other hand, when k`1 " j, the convex hull of τ i,k,j " τ i,j´1,j " ζ i,j´1 ζ i,j over pζ i,j´1,ζi,j q P r0, 1s 2 is given by In the following, we present the proof only for k`1 ă j, and point out that the proof for the case k`1 " j is similar.

C Proof of Proposition 6
Definition 2. Let, D " pV, Aq be a digraph and b P R |V | . A function f : A Ñ R is called as b´transshipment if excess f pv i q :" f @ δ in pv i q D´f xδ out pv i qy " bpv i q @ v i P V , where δ in pv i q Ď A (resp. δ out pv i q Ď A) is the set of all arcs entering (resp. leaving) the vertex v i , and f xδpv i qy :" ř aPδpviq f paq. In our case, the function f paq evaluates the flow along the arc a. Lemma 1 (Rado (1943)). Let D " pV, Aq be a digraph, and let b : V Ñ R with b xV y " 0. Then there exists a b-transshipment f ě 0 if and only if b xU y ď 0 for each U Ď V with δ in pU q " H.

A-2
A PREPRINT -OCTOBER 26, 2020 We now use Lemma 1 to prove Proposition 6.
Consider the digraph D " pV, Aq, where V " D 6 Y D 7 and A " pD 6ˆD7 qztpN`1, 0qu (see §4.2.3 and Figure 5 for definition of D 6 and D 7 ). We have discarded the arc from N`1 P D 6 to 0 P D 7 , because the flow along that arc is zero (see (10)). Observe that, for every n P D 6 , bpnq " excess ψ pnq "´ř mPD7 ψ i,n,m,j " τ i,j,n (see Figure 5). Similarly, for every m P D 7 , bpmq " excess ψ pmq " ř nPD6 ψ i,n,m,j " β m,i,j . Then, b xV y " ř N`1 n"j`1 bpnq`ř i´1 m"0 bpmq "´ř N`1 n"j`1 τ i,j,n`ř i´1 m"0 β m,i,j " 0 (from definition of S i,j ). From Lemma 1, a b´transshipment ψ ě 0 exists if and only if b xU y ď 0 for each U Ď V with δ in pU q " H. For every U Ď D 6 Ă V , b xU y ď 0 is satisfied trivially. On the other hand, U cannot be chosen to be a subset of D 7 , because for every U Ď D 7 , δ in pU q ‰ H. Therefore, in order to derive non-trivial inequalities, we must choose subsets of V containing vertices of both D 6 and D 7 .
Let U " pD 6 ztN`1uq Y t0u. Note that δ in pU q " H. Then, a b´transshipment ψ ě 0 exists if and only if b xU y "´ř N n"j`1 τ i,j,n`β0,i,j ď 0, or It can be verified that for every other subset U Ď V satisfying δ in pU q " H, the inequality ensuring b xU y ď 0 is implied from ř i´1 m"0 β m,i,j " ř N`1 n"j`1 τ i,j,n . Therefore, proj pτ,βq pS i,j q " Indeed, ψ i,n,m,j can be defined to verify that (48) is the projection of S i,j .

D Proof of Proposition 7
In addition to binary variables associated with the presence/absence condensers and reboilers, CG06 has variables for the presence of heat exchanger, which we denote as η i,j . To our model, we add η i,j " χ i,j`ρi,j .
Further, we remark that for i ď k ď j´1, ri, js P P, In Tables 5, 4 and 6, we prove that the set defined by (A2)-(A8), ζ i,j P r0, 1s, @ ri, js P T , ρ i,j P r0, 1s, @ pi, jq P R and χ i,j P r0, 1s, @ pi, jq P C is tighter than CG06, GA10 and TAT19, respectively. We point out that, in GA06, the authors did not consider thermally coupled configurations. Thus, we show the proof only for the constraints they reported.

F Proof of Proposition 10
When x is restricted to x P X, the set S " tpx, y, xyq P D | xy " x¨y, x " xu can be expressed as an affine transformation of y lo ď y ď y up , whose extreme points are y P ty lo , y up u. Therefore, the extreme points of convex hull of S are contained in the set of points where y P ty lo , y up u. Let S 1 " tpx, y, xyq | xy " y lo x, Bx ď b, y " y lo u and S 2 " tpx, y, xyq | xy " y up x, Bx ď b, y " y lo u. Then, by Krein-Milman theorem, convex hull of S is obtained by taking the disjunctive union of S 1 and S 2 , i.e., ConvpSq " proj px,y,xyq tpx, y, xy, x 1 , x 2 , λ 1 , λ 2 q | Bx i ď bλ i , i " 1, 2, (52), λ 1 ě 0, λ 2 ě 0u, where xy " x 1 y lo`x2 y up , x " x 1`x2 , y " y lo λ 1`yup λ 2 , λ 1`λ2 " 1, (52b) Solving the above equations leads to x 1 " y up x´xy y up´ylo , x 2 " xy´y lo x y up´ylo , λ 1 " y up´y y up´ylo , λ 2 " y´y lo y up´ylo .
Using the above relations, we substitute out x 1 , x 2 , λ 1 and λ 2 to obtain the convex hull description in the proposition.
We remove the equality w " w 1`w2 to project out w. Solving the linear equations yields λ 1 " px upx q{px up´xlo q, λ 2 " px´x lo q{px up´xlo q, y 1 " px up y´xyq{px up´xlo q and y 2 " pxy´x lo yq{px up´xlo q. Using these equations, we substitute out auxiliary variables y 1 , y 2 , λ 1 and λ 2 . Finally, eliminating variables w 1 and w 2 using Fourier-Motzkin elimination yields the convex hull description in the Proposition. The outer-approximation of the convex hull follows directly from Remark 8.