Synthesis of distillation configurations. II: A search formulation for basic configurations

Computers and Chemical Engineering 34 (2010) 84–95

Contents lists available at ScienceDirect

Computers and Chemical Engineering journal homepage: www.elsevier.com/locate/compchemeng

Synthesis of distillation configurations. II: A search formulation for basic configurations Arun Giridhar, Rakesh Agrawal ∗ School of Chemical Engineering and the Energy Center at Discovery Park, Purdue University, 480 Stadium Mall Drive, West Lafayette, IN 47907, USA

a r t i c l e

i n f o

Article history: Received 25 October 2008 Received in revised form 4 May 2009 Accepted 5 May 2009 Available online 15 May 2009 Keywords: Multicomponent distillation Synthesis Supernetwork Optimization

a b s t r a c t It is known that there exist a combinatorially large number of multicomponent distillation column configurations for a given feed and given product specifications. In the accompanying part I of the paper, we classified configurations into basic configurations (those that use (n − 1) columns for an n-component feed mixture) and non-basic configurations (those that use more than (n − 1) columns for the same feed). As n exceeds four, the number of non-basic configurations greatly exceeds the number of basic configurations and have a huge impact on the size of the search space. However, through extensive calculations for a four-component feed, we have found that a non-basic configuration never has a lower heat duty than the lowest heat duty basic configuration. While prior researchers have proposed a mathematical formulation that could be used to generate only the set of basic configurations, we present an alternative supernetwork model to achieve this goal. We show how the supernetwork can be reduced to a binary integer program (BIP), and how all feasible configurations can be drawn for a given application using a suitable solver. We present different solution techniques, including parallel algorithms based on the fact that the evaluation of the performance for each configuration is independent of any other configuration. Through such a procedure, all applicable basic configurations can be rank listed for a given application according to a chosen criterion. We also extend the supernetwork formulation to a unified non-convex MINLP that can be used in a flexible way to find the best candidate configurations for a given application. We show the power of our search formulation in accommodating general separation techniques other than distillation, such as membranes, adsorption, etc. © 2009 Elsevier Ltd. All rights reserved.

1. Introduction The synthesis of separation networks for multicomponent feed mixtures has been a challenging problem in process design for several decades. Distillation configurations in particular have been the focus of several research efforts due to their widespread use in the process industries. The problem of finding the best configuration of distillation columns (with associated reboilers, condensers, piping, and control systems) that should separate a given multicomponent feed into individual components is usually a difficult one. This is due to two main reasons: the set of distinct feasible distillation configurations is often ill-defined, and grows combinatorially fast with the number of components that need to be separated; and, the fitness measure of a given configuration is usually a highly non-linear function of the given parameters, causing computational challenges. Due to the industrial importance of multicomponent distillation network synthesis, the problem has been examined by several

∗ Corresponding author. E-mail address: [email protected] (R. Agrawal). 0098-1354/$ – see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.compchemeng.2009.05.004

researchers in the literature (Doherty & Malone, 2001). Many early results were restricted to the case of sharp splits alone, where each column has exactly two product streams (distillate and bottom), and these two streams have no common intermediate component distributing in appreciable quantities. Thompson and King (1972) and King (1980) have examined such configurations. Heuristics for the synthesis have been proposed based on physical insights by Rathore, van Wormer, and Powers (1974), Seader and Westerberg (1977) and Douglas (1988). A need for systematic solution methods in process synthesis has long been recognized (Westerberg, 2004), and there is a need to provide some assurance that a design that has been synthesized is at least locally optimal, if not globally. Being able to provide such an assurance is not a trivial undertaking in the case of distillation synthesis, due to the highly non-linear and non-convex search space and objective functions commonly used. The broad goal of arriving at solutions which one knows to be good may be achieved in the long-term through two complementary approaches. A better search method is the first weapon, one that makes it easier to find good locally optimal solutions than before. A tighter search space formulation may also be used independently of the search method, to aid the search by filtering out solutions that may not be desirable, or by

A. Giridhar, R. Agrawal / Computers and Chemical Engineering 34 (2010) 84–95

providing a representation that makes it easier to locate regions that are likely to have optimal solutions. In this paper, we describe such a formulation and describe computational search methods that we have found to be useful in the synthesis of configurations with low heat duty requirement.

2. Prior work Generally, the problem of distillation configuration synthesis has been attacked by various researchers using two broad approaches. The first is an additive approach, involving the listing of a few configurations by hand through a discovery process, then building rules to essentially interpolate or extrapolate the discovered structures to obtain more configurations. Building configurations sequentially using heuristics is an additive approach: we are assured that the final configuration will be feasible, but it is easy to miss some configurations if the rules are not complete. The second is a subtractive approach, involving the formulation of a search space that embeds various configurations, and eliminating sections of the space that do not correspond to feasible configurations, arriving at a reduced space where every point represents some feasible configuration, including those not discovered by hand. Using a mathematical superstructure is a subtractive approach, which can potentially ensure completeness of the search space, but which can also potentially generate configurations that may not be desirable from a practical viewpoint and should preferentially be removed using constraints. Several authors have examined the concept of a superstructure, which seeks to embed all possible configurations, and which can be used to systematically explore and evaluate the space of all configurations. An early superstructure was proposed by Sargent and Gaminibandara (1976), and allowed non-sharp splits to exist through splitters, mixers and flow control. An extended space was proposed by Agrawal (1996) to include satellite configurations that were not embedded in the Sargent and Gaminibandara superstructure. Agrawal (2003) also provided a systematic rule-based method whereby individual basic configurations can be created by starting at the product end in a network and then marching towards the feed end by making several decision regarding reboilers, condensers and transfer streams. Rong, Kraslawski, and Turunen (2003), and Fidkowski (2006) published the exhaustive space for multicomponent feeds, generated through a combination of individual split possibilities. Hu and Rippin (1991) presented a method for embedding various multicomponent configurations in a single unified framework. A state task network formulation was proposed by Sargent (1998), and extended by Yeomans and Grossmann (1999). Caballero and Grossmann (2001, 2004) have proposed an elegant mathematical generalized disjunctive model for synthesis of multicomponent distillation networks, building on the observations suggested by Agrawal (1996) and the state task network of Sargent (1998). In their formulation, the separation tasks involving pseudo-columns and heat exchangers (reboilers and condensers) are considered independently. Then a set of logical relationships are used to ensure feasibility of the sequence of tasks. In this approach, the heat exchangers are not directly connected with the set of tasks. For example, in the two-step optimization procedure, all the internal heat exchangers associated with mixtures are initially taken to be absent and the resulting configurations with varying degree of thermal coupling are solved to identify an ‘optimal’ configuration. This is then followed with cost and energy tradeoff calculations by adding internal heat exchangers to the identified optimal configuration to further increase the optimality of the solution. This two-step procedure, while being useful in providing an optimal solution, cannot guarantee optimality. It presumes that post-assignment of the internal heat exchangers to a configuration found to be initially

85

optimum with the constraint of no internal heat exchangers, will result in the overall optimum configuration. There is a need for a framework that will directly allow the generation and optimization of (n − 1) column configurations with any degree of thermal coupling. This goal of our work is to introduce a set of logical constraints that can first create and enumerate all the (n − 1) column configurations without any thermal coupling (which we refer to as basic configurations), and use these constraints to find the basic configuration with the lowest heat duty. In a subsequent paper, we will demonstrate its extension to include any desired level of thermal coupling. It should be noted that there have been earlier successful attempts to elucidate basic-only configurations for n-component mixtures (Agrawal, 2003; Caballero & Grossmann, 2006). In this work, our mathematical description of the physical structures of the distillation configurations has some similarities with the prior work in the literature. However, especially when compared to Caballero and Grossmann (2006), our method does differ in modeling the component and stream flows in the distillation configuration, and in the computational methods used to find the configuration with the lowest heat duty. When an n-component feed is separated in n product streams each enriched in one component, the number of configurations increases rapidly when enumerated exhaustively. Moreover, for an n-component mixture, the set not only includes (n − 1) column basic configurations, but also non-basic configurations containing more than (n − 1) columns. In the accompanying paper we showed an example of a non-basic configuration containing nine distillation columns for a six-component feed. These additional columns in a non-basic configuration add cost not only due to additional columns, but also the associated hardware such as reboilers and condensers along with their utility supply systems. For a fourcomponent feed, through extensive calculations for various feed conditions, we found in the companion paper that the six non-basic configurations having four distillation columns never have lower heat duty than the 18 basic configurations with three columns. Furthermore, as the number of components rises beyond four, the subset of non-basic configurations grows much more quickly than the subset of basic configurations. For example, in the case of a six-component mixture, the exhaustive method leads to over four million configurations of which only 4373 are basic and the rest non-basic. This problem of large search space is further exacerbated by orders of magnitude when thermal coupling between distillation columns is considered. In order to make the problem of search space manageable and in the light of heat duty result through extensive calculations for the four-component feed mixtures and the anticipated increased capital cost due to additional columns for the non-basic configurations, we include only basic configurations in a problem formulation for finding an optimal configuration. The reduced search space of basic-only configurations provides other associated benefits. For example, when the search space is reduced in size, it is easier to enumerate several alternative configurations individually and optimize their stream flows to minimize heat duty. This ability to rank list alternative configurations according to a chosen criterion can be quite useful for a practicing engineer to assess benefits of the parameters that cannot be easily represented in mathematical form. For example, for a five or six-component feed, all the configurations within some x% of the lowest heat duty may be rank listed and then a process engineer can pick the desired configuration based on other parameters such as control and operability, manufacturability, shipability of columns, etc., using their experience and knowledge. Such an exercise can also be beneficial in identifying retrofittable configurations for debottlenecking of an existing distillation train. From a computational point of view, a search space that embeds only basic configurations can be used to exhaustively generate all feasible structures (configurations), with each then being subjected

86

A. Giridhar, R. Agrawal / Computers and Chemical Engineering 34 (2010) 84–95

to a local optimization of the objective function with respect to its parameters (molar flow rates). The optimized parameters for each of the basic configurations then provide an opportunity to not only find the configuration with the best value of the objective function but also rank list all these configurations according to their relative performance. An exhaustive search over the space of all possible structures is especially attractive when a distributed computing framework is available to exploit the parallel nature of the problem. The problem may also be posed as a single large mixed integer non-linear program (MINLP), which may be solved with suitable solvers that can handle non-convex functions and constraints. While we have formulated the problem by all methods and the ultimate choice of the method is dependent on the user’s preference, we have extensively used the parallel nature of the problem by calculating the optimized value of the parameters for each of the basic configurations from the set of n-component basic configurations. Our initial focus in this work is on obtaining a compact formulation that includes only basic configurations. We subsequently describe how to use this formulation to find distillation configurations with low heat duty requirement. The exhaustive methods of Rong et al. (2003) and Fidkowski (2006) provides no guideline on how to exclude rapidly growing search space associated with nonbasic configurations as the number of components are increased in a feed mixture. Agrawal (2003) recently provided a set of rules starting at the product end to generate only basic configurations. However, it is desirable to generate a set of equations describing basic-only configurations that are independent of whether one starts at the product end or the feed end. Furthermore, a good formulation should allow an engineer to impose process configuration constraints, such as the ones associated with the use of reboilers and condensers as well as transfer of intermediate streams at not only the feed or product ends but also at any intermediate point in the distillation sequencing. The formulation should use a simple set of logical constraints to exclude non-basic configurations and be easy to use as mathematical constraints in the problem solving exercise. While Caballero and Grossmann (2006) have recently proposed a method, we now present an alternate useful formulation to create all feasible basic configurations. The method is valid for configurations both with and without thermal coupling, is compact and easy to use in an exercise to identify candidate configurations. It is of interest to note that the formulation can be easily extended for other separation methods that employ binary partitioning of feed streams, such as in membrane cascades for multicomponent separation.

3. Formulation Our formulation, which we call a supernetwork formulation, derives from graph–theoretic representations of supply chain networks. Accordingly, as shown in Figs. 1 and 2, nodes in the supernetwork represent a stream such as feed, intermediate mixture and a final product stream in the distillation network. Edges represent column sections, with upward-sloping edges (distillate edges) representing rectifiers and downward-sloping edges (bottoms edges) representing strippers. All edges are directed in the formulation. (We show later that bidirected edges can be used to represent thermally coupled configurations.) Two nodes are connected by a directed edge if the target node can be made as a distillate or bottom product from the source node in one step. An instance of a configuration is represented as a subgraph of these nodes and edges, as shown in Fig. 3. As with other formulations, the supernetwork formulation restricts edge combinations to represent valid configurations alone. We choose to represent only the presence or absence of edges

Fig. 1. A supernetwork for the separation of a three component feed. All edges are directed from left to right.

explicitly in this formulation. Other formulations are also possible with nodes alone, or with a combination of nodes and edges, as done by Caballero and Grossmann (2006). In this respect, our formulation differs from the prior work in that we have only one type of discrete variable, representing the existence of edges, and not two distinct types that represent existence of edges as well as nodes. The representation of edge existence in our work is achieved by associating a binary variable x with an edge, such that the binary variable has a value 1 when the edge exists, otherwise its value is zero. We now list three categories of constraints that when used together ensure generation of only basic configurations: 3.1. The first constraint: ensure mass balance If a node exists (i.e., that stream is produced), and it is a nonterminal node (i.e., it represents a stream which is not a final product), then it must have exactly one distillate and one bottom. For example, in a quaternary separation, the stream ABC may either exist or not exist. If it exists, we know that it must have been produced from some other stream that could produce it in one step. This corresponds to the incoming edges of the node ABC in the supernetwork. If this edge exists, then there must be exactly one distillate and one bottom for ABC. Suppose we examine all edges

Fig. 2. A supernetwork for the separation of a four-component feed.

A. Giridhar, R. Agrawal / Computers and Chemical Engineering 34 (2010) 84–95

87

Fig. 3. Correspondence between subgraph representations and physical distillation configurations.

incident on node ABC, as shown in Fig. 4, we find that the following relations suffice. xABC→AB + xABC→A = xABCD→ABC xABC→BC + xABC→C = xABCD→ABC If ABC is produced, then the RHS of both equations becomes 1, and there is forced to be exactly one distillate edge and one bottom edge leaving ABC. If ABC is not produced, then the RHS is 0, and there is forced to be exactly 0 distillates or bottoms (which is another way of saying that if a stream is not produced in the first place, it cannot produce anything else in turn). This rule can be generalized for all non-terminal nodes in the supernetwork, giving: If there is any edge feeding into a non-product node, then it must have exactly one distillate edge leaving as well as exactly one bottom edge leaving. 3.2. The second constraint: prevent redundant configurations Given that a node has different choices of distillate product and bottom product, not all combinations of these products are valid. We assert here that if a particular component is present in a stream,

Fig. 4. All edges incident on node ABC.

then it must be present in the products from that stream, either in distillate or in bottom or in both. We choose to do this compactly, as opposed to explicit enumeration of allowed splits, by recognizing that for a split to be valid, there must be at least as many components leaving in the product streams as there were entering in the feed streams. Consider the different distillate and bottom choices for the stream ABCD in a four-component separation, as shown in Fig. 5. The distillate could be A, AB or ABC, and the bottom could be D, CD or BCD. Out of the nine combinations (3 distillate, 3 bottom) for the pair, only six are allowed: A/BCD, AB/BCD, AB/CD, ABC/BCD, ABC/CD, and ABC/D. The three disallowed choices are A/CD, AB/D and A/D. The difference between the allowed and disallowed choices is that in the allowed choices, all components that were fed in are present in the product streams. In the disallowed choices, one or more components “disappear” and violate mass balance. Note that in the six allowed cases, if we add the number of components present in the distillate and the bottoms streams, the total is at least 4 and at most 6. In general, for such configurations, this total ranges from n for sharp splits to 2(n − 1) for the most sloppy split. Given that we can represent the existence of an edge by a binary variable, the number of components being produced along a given edge is the product of two quantities: the number of components in the node to which it

Fig. 5. Stream ABCD can give one of ABC, AB or A as distillate, and one of BCD, CD or D as bottom. Not all such combinations are possible however.

88

A. Giridhar, R. Agrawal / Computers and Chemical Engineering 34 (2010) 84–95

is going, and its own binary existence variable. Using this, we get the following relationship for ABCD: xABCD→A + 2xABCD→AB + 3xABCD→ABC + 3xABCD→BCD + 2xABCD→CD

same time for the same node. This can be generalized as follows: sum(xin,distillate ) = 1,

for all nodes

sum(xin,bottoms ) = 1,

for all nodes

+ xABCD→D ≥ 4 In general, the node from which the candidate distillates and bottoms can be produced can itself not be produced in the first place. We therefore multiply the RHS by an appropriate sum of existence variables of that node’s incoming edges. sum(xout × componentsout ) ≥ sum(xin,distillate ) × componentsnode sum(xout × componentsout ) ≥ sum(xin,bottom ) × componentsnode where xin,distillate denotes all edges in the supernetwork that make the stream node as a distillate, and xin,bottom is defined analogously. 3.3. The third constraint: ensure basicity Ensure that the configuration generated by the other constraints can be implemented in the minimum number of (n − 1) columns. This can be transformed into an equivalent constraint on the location of the product streams (distillate, bottom or side draw). A distinguishing feature of all non-basic configurations, is that at least one stream is produced as two heavies (i.e., bottoms product) or as two lights (i.e., as distillate products), resulting in more columns being needed than the minimum. Consider the two intermediate components B and C for a four-component case as shown in Fig. 6. Accordingly, we prevent this using the following constraints for the quaternary case: xBC→B + xBCD→B ≤ 1 xBC→C + xABC→C ≤ 1 What we do here is to prevent more than one incoming distillate edge (or one incoming bottoms edge) from being present at the

However, to ensure that streams can be withdrawn as side draws, we also add the constraint: sum(xin,distillate ) + sum(xin,bottoms ) ≥ 1 It should be mentioned that these are not the only rules that ensure basicity. Basicity may be achieved through other constraints also, but the ones presented above turn out to be the simplest to implement for an automated search. It can also be restated in a different form, which is very useful when drawing configurations by hand from the feed end towards the product end. Given a partially drawn configuration (starting from the feed end), we check the terminal streams we have produced so far. For example, we may take a fivecomponent feed ABCDE and decide to split it as ABCD and CDE, which are our current terminal streams. Now we look at the common subsets of these two streams (C, D and CD), which further have at least two components (only CD). The rule to prevent non-basicity is to ensure that this stream CD is produced either as a side draw or not at all. Mathematically this means that the node CD in the supernetwork has either two incoming edges (one being distillate and the other being bottoms), or no incoming edges (not produced), which can be further generalized to its indegree being even (not odd). This rule should be applied recursively after each split decision is made. In the example, if the rule is violated and CD is made but not as a side draw, then it must be either a distillate from CDE or a bottom from ABCD. If CD is produced as a distillate from CDE, then D must be made as a heavy from CD and it must be made as a heavy from either ABCD or BCD (neither of which can make CD since CD is a distillate). Since D is made as two heavies, the configuration is non-basic. Alternatively, if CD is produced as a bottom from ABCD (or indirectly through BCD), then CDE and CD must both produce C as a distillate, which means that C is made as two distillates, making the configuration non-basic. We also provide extensions to the earlier work (Agrawal, 2003) in the Appendix A, to make a complete set of rules that a human designer could use to ensure basicity of the configurations. Combining all the mathematical constraints together, for a three-component supernetwork as shown in Fig. 1, the equations are: Mass balance: xABC→A + xABC→AB = 1 xABC→BC + xABC→C = 1 xAB→A − xABC→AB = 0 xAB→B − xABC→AB = 0 xBC→B − xABC→BC = 0 xBC→C − xABC→BC = 0 Prevent redundant configurations: xABC→A + 2xABC→AB + 2xABC→BC + xABC→C ≥ 3

x ∈ {0, 1} Fig. 6. Non-basicity arises when one or more streams are produced from at least two distillate edges or two bottoms edges.

For a four-component supernetwork as shown in Fig. 2, the equations are:

A. Giridhar, R. Agrawal / Computers and Chemical Engineering 34 (2010) 84–95

89

Mass balance: xABCD→A + xABCD→AB + xABCD→ABC = 1 xABCD→BCD + xABCD→CD + xABCD→D = 1 xABC→A + xABC→AB = xABCD→ABC xABC→BC + xABC→C = xABCD→ABC xBCD→B + xBCD→BC = xABCD→BCD xBCD→CD + xBCD→D = xABCD→BCD xAB→A = xABCD→AB + xABC→AB xAB→B = xAB→A xBC→B ≥ xABC→BC xBC→B ≥ xBCD→BC xBC→B ≤ xABC→BC + xBCD→BC xBC→C = xBC→B xCD→C = xABCD→CD + xBCD→CD xCD→D = xCD→C Prevent redundant configurations: xABCD→A + 2xABCD→AB + 3xABCD→ABC + 3xABCD→BCD + 2xABCD→CD + xABCD→D ≥ 4 xABC→A + 2xABC→AB + 2xABC→BC + xABC→C ≥ 3xABCD→ABC xBCD→B + 2xBCD→BC + 2xBCD→CD + xBCD→D ≥ 3xABCD→BCD

Fig. 7. A pictorial representation of the number of nodes and edges in the supernetwork. If each cube in the diagram is assumed to rest on a tower of identical cubes, the total number of cubes equals half the number of edges in the supernetwork for n components where n is the number of cubes along any coordinate direction. Projecting the stack of cubes onto any coordinate plane and counting the number of resulting squares gives the number of nodes for the same supernetwork.

4. Results and discussion

Ensure that products are made:

4.1. Configurations

xABCD→A + xABC→A + xAB→A = 1 xBCD→B + xBC→B + xAB→B ≥ 1 xABC→C + xBC→C + xCD→C ≥ 1 xABCD→D + xBCD→D + xCD→D = 1

The supernetwork formulation was tested using MATLAB and GAMS/CPLEX. The formulation was proved to be correct and complete by brute force generation of all possible configurations (basic and non-basic) for up to six components. The results are summarized in Table 1. Regarding comparison with the previous work of Agrawal (2003), please see the Appendix A. The number of basic configurations we obtain also agrees with the results of Caballero and Grossmann (2006). As described in part I of this paper, the ranges for the thermally coupled exhaustive case describe the bounds on the number of configurations obtained by disregarding multiplicity (lower bound) and considering multiplicity (upper bound). Note that for n > 4, the basic configurations are a very small subset of the total number of configurations, our supernetwork formulation is expected to be beneficial in reducing the search space. The search space with only the basic configuration should be more manageable in locating solutions that are optimal or close to optimal.

Basicity constraints: xBCD→B + xBC→B ≤ 1 xABC→C + xBC→C ≤ 1 x ∈ {0, 1} The size of the supernetwork formulation grows polynomially with respect to the number of feed components n. The number of variables is equal to the number of edges in the supernetwork, which is a cubic function of n: Number of variables =

(n − 1)n(n + 1) 3

The number of constraints is proportional to the number of nodes, since mass balance, basicity, etc., are defined with respect to streams and not sections. The number of nodes is a quadratic function of n, namely n(n + 1)/2. As a point of mathematical interest, the number of edges may be represented diagrammatically on a cube as shown in Fig. 7, with each small cube assumed to rest on a tower of identical cubes. The number of such small cubes equals half the number of edges. By projecting the tower of cubes onto any of the coordinate planes and counting the number of squares, one also gets the number of nodes. The formulation may be very easily extended to generate all thermally coupled configurations by introducing new binary variables yi , which indicate whether a particular edge is bidirected or not. All constraints to be obeyed are yi ≤ xi , for all edges i that do not make final products. What this means is that if an edge exists in the basic configuration (the corresponding xi is 1), then it has a reboiler or a condenser by default. The additional variable yi represents whether this reboiler (or condenser) has been replaced by a vapor (or liquid) stream from a downstream column. A detailed formulation incorporating thermal coupling will be presented elsewhere.

4.2. Optimal configurations After having solved the problem of finding all configurations, the problem now becomes how to find the best configurations for a given feed. This question naturally leads from a search space formulation to an objective function definition and an optimization method. As our current work is focused on finding basic configurations with low heat duty, we have chosen our objective function to be the overall vapor duty for the entire configuration (summed over all columns), assuming that the columns are all thermodynamically pinched and at minimum reflux. We evaluate the overall vapor duty at minimum reflux using Underwood’s method (1946). Our overall approach is shown in Fig. 8, and described in greater detail in the following paragraphs. The problem therefore reduces to finding the configurations and their component molar flows, that give the lowest overall vapor duty with respect to both the physical structure and the parameters describing that structure. The classic approach to the problem of simultaneously optimizing both the structure and its parameters is to pose the problem as an MINLP, which is one of our approaches. The discrete variables are used to represent the structure (binary variables associated with the existence of sections, and

90

A. Giridhar, R. Agrawal / Computers and Chemical Engineering 34 (2010) 84–95

Table 1 Number of feasible configurations. Number of components

Configurations Exhaustive

3 4 5 6

3 24 1,352 4,074,863

Thermally coupled Basic only 3 18 203 4,373

grouping sections into columns), while additional continuous decision variables zij and qi (where i represents the edge as before, and j represents the component) are introduced. The variable zij represents the molar flow of component j associated with edge i, and qi represents its thermal quality. (It is conceptually easier to visualize qi , but in actual practice we find it more computationally robust to use liquid and vapor molar flows of each component, and compute total molar flows and thermal quality from them as needed. For conceptual clarity, however, we will continue to refer to qi as though it were a decision variable instead of a variable calculated from other decision variables.) Further continuous variables, specific to Underwood’s method, are also used, representing the roots of Underwood’s equations for each stream. We add the following constraints:

• Linking constraints to ensure that zij is finite only when edge i exists. These are of the big M type, and impose an upper boundbased on feed molar flow of the component and a lower bound to prevent minuscule composition ratios of components in any stream. For instance, if the edge xABCD→ABC exists, then the corresponding molar flows zABCD→ABC,A , zABCD→ABC,B and zABCD→ABC,C must all be finite. Moreover, we can assert that zABCD→ABC,A must be (in this case) exactly equal to the molar feed rate of component A, while zABCD→ABC,B and zABCD→ABC,C are constrained to be less than or equal to their corresponding component feed rates. All three variables zABCD→ABC,A , zABCD→ABC,B and zABCD→ABC,C are also constrained to be greater than or equal to some small value (we use 10−3 ) to prevent the structure asserting a certain set of components in one stream and the actual molar flows representing

Exhaustive 5 226–288 72,769–191,240 914,944,696–9,368,313,073

Basic only 5 134 5,925 502,539

some other smaller set of components. • Mass balance and composition constraints over zij . Mass balance is straightforward, and composition constraints refer to certain splits such as ABCD → ABC/BCD where at least two components distribute to both distillate and bottoms. In such cases, we find it computationally prudent to assert that the B:C ratio in the distillate ABC should not be less than the B:C ratio in the feed ABCD. This is to prevent the optimization solver from assigning values to the variables which would render the split closer to AC/BD, which is not physically possible. In other words, it is mathematically possible (but not physically) to send most or all the B in the feed to the bottom stream and simultaneously send most or all the C in the feed to the distillate stream, while still obeying mass balance. To prevent such oddities, we constrain the B:C ratio in the distillate to be at least as large as the B:C ratio in the feed. For three components, such composition constraints are not required. For four components, they result in one linear equation, for those configurations where both B and C distribute to both distillate and bottom in the first column. For five or more components, they result in a set of multilinear equations involving the product of two variables, which in turn results in a non-linear and non-convex space. • Equilibrium constraints for two-phase streams. When drawing the stream BC as a two-phase side draw, we require that the compositions of B and C in the liquid and vapor phases be in equilibrium. This constraint is simplified considerably by using separate variables for liquid and vapor component flows than a total flow and a thermal quality factor. • Heat balances to account for different liquid and vapor flows. Two types of heat balance are required. The first type is required in our

Fig. 8. A flowchart of the overall optimization strategy used in our work. The method is essentially an exhaustive search over all configurations, with the subproblem of minimizing the overall vapor duty of each configuration with respect to its stream molar flows being solved as a local NLP search. Useful initial guesses are found by minimizing a dummy objective function, namely a (convex and unimodal) paraboloid centered between the upper and lower bounds.

A. Giridhar, R. Agrawal / Computers and Chemical Engineering 34 (2010) 84–95

formulation only due to modeling liquid and vapor molar flows separately. A model that treats thermal quality as an actual decision variable would not require these constraints. A second type (needed irrespective of the model) is a vapor balance around each section of a column, which is used to compute the reboil duty of the column after offsetting any vapor feeds to the column. For example, if a column accepts ABC and BCD as feed streams, and produces AB, BC and CD, then the vapor present in ABC can be used to make part of the vapor needed in stream AB. Similarly the vapor component of feed BCD can be used to make part of vapor products AB and BC, depending on how much vapor is available. Note that only feed streams that are below a product stream (i.e., higher temperature) may contribute vapor to that product stream. Any balance vapor duty not met by incoming feed streams is then provided by the reboiler. • Any additional constraints specific to a certain process instance, governed by process or market factors, which we term “user constraints”, in addition to the “model constraints” listed above. Such constraints may restrict the structural search space of configurations, or also restrict the molar flows or compositions of certain streams, if present. Given the discrete variables describing the structure of the configuration and the continuous variables describing the parameters of the structure, such as molar flows, etc., we posed the entire problem as a single monolithic MINLP (simultaneously optimizing the structure and the continuous parameters to minimize vapor duty) and were able to solve problems globally for three and four components (Giridhar, 2008) using GAMS/BARON (Tawarmalani & Sahinidis, 2002). The runtimes for four components were found to be typically between 5 and 7 min on a single processor to find the best one configuration; additional configurations take more time. Running the MINLP model as a sequence of NLP models by freezing integer variables in place improved the speed to anywhere from a few seconds to 4 min per configuration (depending on both the structure and feed parameters). However, the MINLP approach was found to be very sensitive to the initial guess values, and in some cases failed to find a solution that satisfied all the constraints, or reported a solution that was known to be suboptimal. For five components, the unified MINLP was found to take inordinate time, indicating a need to make the MINLP more robust and more efficient in the context of global optimization. Accordingly, we have worked around this problem using other methods as described in the following paragraphs. Given that for four and five components there are not more than a few hundred feasible configurations present in the search space, and given that the optimization process for each configuration is completely independent from that of every other configuration, we recognize this problem is amenable to parallel computation. This is made particularly easy by the availability of multicore and multi-threaded processors, and by the fact that we require only process-level parallelization and not operation-level parallelization. In other words, if we decide to do an exhaustive search over all configurations, the order of computation is completely arbitrary, and no configuration needs the results of any other configuration. An exhaustive search over all configurations on a single computer may still take prohibitive time, but our computational experience indicated that this was solely due to the number of configurations and not due to the complexity of the parametric optimization, particularly if the optimum is not required to be guaranteed global. Accordingly, we split the problem into an outer binary integer program (BIP) over the configurations and an inner non-linear program (NLP) over the parameters (molar flows, etc.) within a configuration. Using a four-component example, there are 18 basic configurations for a four-component feed ABCD. In the MINLP approach described earlier the objective was to simultaneously find the

91

best configuration among the 18 by varying the discrete variables (xABCD→ABC , etc.) representing the edges, and also assign the most suitable values to its molar flow and thermal quality variables. In the parallel computation approach, we split up the single MINLP into 18 smaller NLPs, each corresponding to one of the 18 configuration structures. For each NLP, the discrete variables representing the structure are fixed to the corresponding set of values (1 if the corresponding edge is present, 0 if absent), and the optimization is done with respect to only the molar flows and thermal quality variables. The best configuration then is the one that has the lowest value among all 18 NLP results. Similar procedures are carried out for feeds with three and five components as well. The constraints for each NLP instance are generated programmatically based on the physical structure of the corresponding configuration. We solved the NLP problems using the standard Optimization Toolbox routines in MATLAB, namely the default local search method provided by the constrained minimization routine fmincon. The performance for this distributed exhaustive approach was found to be faster than the time taken to obtain the first optimized configuration by solving the problem as one unified MINLP using a similar local search-based solver. We believe that while such an approach using multiple independent NLPs may not be as mathematically pleasing as a single MINLP, it is nevertheless one that we believe may have become more computationally relevant than previously. Two reasons for this include much easier parallelization due to multicore processors, and also that many large industrial distillation problems can be made more amenable to computation by being reduced to five lumped components: the light key (LK), the heavy key (HK), the lighter-than-light-key (LLK), the heavierthan-heavy-key (HHK), and the distributing components. We have also encountered problem instances where several configurations were clustered near the configuration with the lowest heat duty, all with very similar heat duties, making it difficult to separate the truly global minimum configuration from the others. Therefore, we believe that our approach of optimizing each configuration separately as an independent NLP and rank listing them is attractive. It should be mentioned that obtaining subsequent solutions in the MINLP takes additional time, while the distributed exhaustive computing approach gave the results for all configurations in comparable time. It should also be mentioned that using a global solver in any stage of the problem (whether on the discrete variables or the continuous variables) vastly increases the time required to find the optimum solution. An added advantage of this approach is that we have performance rank lists of all the configurations and then evaluate those that are “similar” in heat duty through other criteria, such as manufacturability, operability, heat integrability with the rest of the plant, etc. One problem with using a local search method for an NLP is that if the NLP is not unimodal, no mathematical guarantee can be made about whether a solution reported is global or not. We have chosen to work around this problem by restarting the local search from different initial guesses, and selecting the best answer from all the runs. While this is not a mathematical guarantee of globality, the probability of not finding the global optimum asymptotically decreases to zero with the number of restarts. For practical considerations of CPU time (i.e., to avoid trying an inordinate number of initial guesses for each NLP), we find that good initial guesses can be obtained by solving a subproblem with the same constraints as the main NLP to be solved, but with the objective function replaced by various paraboloids instead of the overall vapor duty. In other words, we find a point that is as close as possible to the “center” of the feasible search space by centering a paraboloid over the midpoint of the upper and lower bounds and finding its minimum feasible point. This process takes much less time than minimizing the vapor duty over the same search space. The solution to the paraboloid-minimization problem thus obtained serves

92

A. Giridhar, R. Agrawal / Computers and Chemical Engineering 34 (2010) 84–95

Table 2 Feed composition and lumping for cyclopentane case study. Feed

Feed, wt%

Feed, mol%

BPT (◦ C)

Relative volatility

A

n-Butane i-Pentane n-Pentane

0.040 1.791 14.816

0.053 1.909 15.795

−0.6 27.95 36.3

12.723 5.629 4.265

3.413

B

2,2-Dimethylbutane Cyclopentane

0.219 36.716

0.195 40.261

49.7 49-50

2.863 2.789

2.146

C

2,3-Dimethylbutane 2-Methylpentane 3-Methylpentane n-Hexane

1.483 13.502 4.687 11.373

1.323 12.051 4.183 10.151

58.0 60.2 63.2 69.0

2.174 2.029 1.813 1.520

1.398

D

Methyl-cyclopentane Benzene Cyclohexane

13.174 0.428 1.771

12.038 0.421 1.618

71.8 80.1 80–81

1.345 1.170 1.000

1.000

as a feasible initial guess for the vapor duty minimization. There is no mathematical guarantee of this procedure yielding the global minimum for general non-convex problems, but we have found it to yield good results in practice. To ensure that the solutions reported by the distributed NLP approach were of good quality, we examined them using two different crosschecks: a grid search routine and a particle swarm routine. This was done only in a research context and not as potentially suitable methods for practitioners. The grid search routine was simply to discretize the space of the continuous variables into a sufficiently fine grid (typically 1000 points per dimension), and evaluate the objective function at all legal grid points. Performance is of course improved by choosing the decision variables and calculated variables carefully. In the particle swarm approach, we searched the continuous space using a swarm of 20 particles connected as a 4 × 5 von Neumann grid, initialized randomly in a uniform distribution inside the space. Each particle in the swarm changed its trajectory based on its own history and those of its topological neighbors, iterating to convergence in both function value and in swarm diameter. In all the cases that we have seen, we have found the local search NLP routines to give the same global solution as reported by the grid search routine, which is also in agreement with the values reported by the particle swarm routine. We have not investigated the use of hybrid methods in this problem, but we believe that the use of metaheuristic optimization methods (particularly particle swarms, simulated annealing and possibly genetic algorithms) followed by a local search method may be a useful complement to solving the monolithic MINLP. Even in the case of larger instances, where six or more components could be separated into individual products, the full search space of basic configurations may not be realized due to external user constraints. Examples of such user constraints could be “Product B should not be produced as a bottom stream”, possibly because B may be required at high purity without reboiler residue, or “Product C should be recovered as early as possible” and so on. Each of these external constraints when incorporated into the supernetwork reduces the search space significantly; cumulatively they may reduce the space drastically. Accordingly, the parallel computation approach could still be relevant for scales larger than toy problems.

4.3. Case study We ran the parallelized version of the algorithm on a hydrocarbon separation problem in the production of cyclopentane (Kataoka, Noda, Nakaiwa, & Horiuchi, 2006). The feed has 12 components, as shown in Table 2, which have been lumped into four components using the rules proposed by Glinos and Malone (1984), whereby lumped mole fractions are the sum of mole fractions of individual components, and lumped relative volatilities are the

mole fraction weighted average of individual relative volatilities. The lumped cuts are largely the lights, the cyclopentane product, heavy straight chain alkanes and heavy ring hydrocarbons. We evaluated the minimum vapor duty using Underwood’s method for pinched columns. Using a feed flow rate of 0.4 mol per time, we are able to rank order all 18 configurations on the basis of total network vapor duty as shown in Fig. 9. The results indicate a strong variation of vapor duty requirement with the configuration chosen, and a large difference (on the order of 70%) between the sharp split configurations and the sloppy split configurations. In this case, the difference between the lowest and second lowest heat duty configurations is about 13%, indicating that only one configuration may be a strong candidate for further consideration. The rank listing also provides a designer a means of evaluating tradeoffs between heat duty and competing performance measures, such as capital cost or exergy. We evaluated the lowest vapor duty sharp split configuration and the lowest vapor duty non-sharp configuration in ASPEN, using rigorous thermodynamics simulations, for a feed rate of 1650 kg/h. We obtained a result of 2.585 million Btu per hour for the sharp split configuration, and 1.232 million Btu per hour for the same feed in the non-sharp split configuration. These results are in remarkable agreement with what was predicted using Underwood’s method, especially when one considers the approximate nature of lumping and constant relative volatility assumptions used with the method.

5. Extensions beyond distillation Since our supernetwork formulation is completely based on binary arithmetic and the constraints are completely independent of distillation, we find a powerful extension for our supernetwork. Consider a different separation method, such as membrane separation. We could argue that the synthesis of multicomponent membrane separation networks is similar to the synthesis of distillation networks, with the differences being that the feed components are separated on the basis of their permeability rather than their volatility, and the performance measure is compressor duty rather than heat duty (Agrawal, 1996b). This analogy is made clear in Fig. 10. For example, the membrane cascade sections from feed ABC to mixtures AB and BC can be thought of as edges associated with binary variables xABC→AB and xABC→BC in Fig. 1, respectively. Similarly, other edges of the ternary supernetwork in Fig. 1 can be located in the membrane cascade of Fig. 10. Clearly for a multicomponent mixture separation using membranes, one can use the supernetwork formulation to generate all the equivalent basic configurations using membranes alone. As long as a separation method separates a mixture into two fractions one “lighter” and one “heavier”, we can denote the path to the “lighter” fraction as one edge and the path to the “heavier”

A. Giridhar, R. Agrawal / Computers and Chemical Engineering 34 (2010) 84–95

93

Fig. 9. Rank order of the 18 four-component basic configurations for the cyclopentane case study. The first five configurations shown are sharp split configurations. Numbers represent total vapor duty of a configuration.

fraction as another edge in the supernetwork formulation. Thus, multicomponent adsorption processes based on adsorptivity, centrifuge separation based on centrifugal forces, solvent extraction processes based on absorptivity, etc., are amenable to the process synthesis procedure described here. It is also possible to draw feasible hybrid separation processes. In such a supernetwork formulation, an edge will represent all feasible separation processes and the designer will have a choice to pick the optimal one. For example, in a ternary separation, either all or some of the edges in Fig. 1 may simultaneously represent a distillation column section or a membrane cascade as discussed

for Fig. 10. Therefore, ternary hybrid configurations using membranes and distillation columns can be easily generated. This will provide a process designer a larger array of choices. Clearly, possibilities grow very rapidly if additional separation alternatives besides membrane and distillation were to be incorporated in the search for an optimal hybrid scheme. However, while generating hybrid schemes, one expects that process and operating constraints may significantly prune the search space. The number of hybrid configurations from a supernetwork formulation is extremely large. Each basic configuration can be seen to have anywhere from 2(n − 1) to n(n − 1) edges from the supernet-

Fig. 10. Analogy between membranes (a) and distillation (b) showing the similarity of large-scale component separation decisions, and the opportunity to design hybrid configurations that use more than one separation method.

94

A. Giridhar, R. Agrawal / Computers and Chemical Engineering 34 (2010) 84–95

work representation, corresponding to distillation column sections or their equivalent. In the case of a hybrid separation network, each of these edges could be one of several different process units based on the number of the separation techniques available. If a process designer has K different techniques available to choose from, then each configuration as represented in the supernetwork could have anywhere from K2(n−1) to Kn(n−1) realizations. Multiplying the number of such realizations by the number of valid basic configurations results in a very large number of configurations that must be pruned, evaluated and optimized. However, the supernetwork formulation does allow them to be represented for further manipulation. 6. Conclusion We present here an alternate method for the generation of a complete set of only (n − 1) distillation column configurations for the separation of an n-component mixture into n product streams. It compares well to certain exhaustive methods proposed in the literature that generate non-basic configurations containing more than (n − 1) columns for values of n > 3. As the number of components in the feed is increased beyond four, the number of such non-basic configurations severely outweighs the number of basic configurations. The impact of this difference is even further exacerbated when other aspects such as thermal coupling between the columns

are included in the search space containing non-basic configurations. In our previous paper, we showed that additional columns in the non-basic configurations do not lead to reduction in overall heat duty when compared to a low heat duty basic configuration. This presents a need for a mathematically defined search space that would include all feasible basic configurations and exclude the non-basic configurations. In this work, we presented such a mathematical formulation which ensures that only basic configurations are generated. The formulation is compact, and can be used for any number of feed components. We used this formulation to generate all basic configurations in an efficient manner for up to six components, and found the results to agree with prior results available in the literature. We found that the supernetwork formulation is easily extensible as part of an optimization framework to identify the best configurations on the basis of chosen performance measures. The method is reliable when minimizing the total vapor flow using Underwood’s method at minimum reflux. The mathematical structure of the problem lends itself to a clean decomposition between listing all feasible structures and finding the optimal structures. We found that this problem could be solved either as a traditional MINLP problem, or through the suitable use of parallel computation. We find that parallel computation is extremely suitable for this problem due to its inherently parallel structure. We successfully applied this method to identify a three

Fig. A1. Configurations missing from prior published work (Agrawal, 2003), for a five-component feed. The common feature of these configurations is that they all have at least 1 four-component stream that is split sharply into 2 two-component streams.

A. Giridhar, R. Agrawal / Computers and Chemical Engineering 34 (2010) 84–95

column configuration to recover cyclopentane from a feed mixture. In this case, we were able to rank list various configurations based on their total heat duty requirements. A powerful consequence of the proposed supernetwork method is that it can be used to generate feasible networks for other separation techniques, such as membranes, adsorption, absorption, etc. An interesting aspect of our work is that hybrid separation networks, which use more than one separation technique, can also be easily generated using the supernetwork formulation. Acknowledgement The authors would like to acknowledge the Department of Energy (DOE grant DE-FG36-06GO16104) for providing financial support. Special thanks are also due to Dr Z.T. Fidkowski, Professor Venkat Venkatasubramanian and Dr Jeff Siirola. Appendix A. It should be mentioned here that the method outlined in Agrawal (2003) is found to be complete only up to distillations of four components. For five components in the feed, it is found to omit configurations that split quaternary streams directly to two binary streams without any ternary streams present in the system, and direct analogs of such systems, as shown in Fig. A1. Such configurations for mixtures containing five and more components can be drawn through an addendum to Observation No. 8 of the original method: For an n-component mixture, whenever in the final modified network (i.e., after application of Step 1 and Step 2 cases a–f), an internal submixture at depth m (where m < n) has submixtures containing more than one component as its immediate neighbors at the height above and below (but at the same depth m), then the internal submixture can be eliminated from the network. Furthermore, for m < n − 1 where such submixtures contain more than two components and include a submixture from the upper branch (and/or lower branch) then it may be feasible to eliminate the submixture from the upper branch (and/or lower branch) along with its immediate neighboring internal submixture . This is certainly true when the submixture on the upper branch (lower branch) at a depth one greater than the upper branch (lower branch) submixture being eliminated can be separated into submixtures that contain all its components. Thus in Fig. A1a, submixtures ABC from the upper branch and CDE from the lower branch are eliminated along with the internal submixture BCD. In Fig. A1c, only the submixture ABC from the upper branch is eliminated with the internal submixture BCD. Other configurations in Fig. A1 are similarly generated.

95

References Agrawal, R. (1996). Synthesis of distillation column configurations for a multicomponent separation. Industrial & Engineering Chemistry Research, 35, 1059. Agrawal, R. (1996b). Membrane cascade schemes for multicomponent gas separation. Industrial & Engineering Chemistry Research, 35, 3607. Agrawal, R. (2003). Synthesis of multicomponent distillation column configurations. AIChE Journal, 49, 379. Caballero, J. A., & Grossmann, I. E. (2001). Generalized disjunctive programming model for the optimal synthesis of thermally linked distillation columns. Industrial & Engineering Chemistry Research, 40, 2260. Caballero, J. A., & Grossmann, I. E. (2004). Design of distillation sequences: From conventional to fully thermally coupled distillation systems. Computers & Chemical Engineering, 28, 2307. Caballero, J. A., & Grossmann, I. E. (2006). Structural considerations and modeling in the synthesis of heat integrated thermally coupled distillation sequences. Industrial & Engineering Chemistry Research, 45, 8454. Doherty, M. F., & Malone, M. F. (2001). Conceptual design of distillation systems. Boston: McGraw-Hill. Douglas, J. M. (1988). Conceptual design of chemical processes. New York: McGraw-Hill. Fidkowski, Z. T. (2006). Distillation configurations and their energy requirements. American Institute of Chemical Engineers Journal, 52(6), 2098. Giridhar, A., (2008). Synthesis of multicomponent distillation configurations. PhD thesis, Purdue University. Glinos, K., & Malone, M. F. (1984). Minimum reflux, product distribution, and lumping rules for multicomponent distillation. Industrial & Engineering Chemistry Process Design and Development, 23, 764. Hu, Z. B. C., & Rippin, D. W. T. (1991). Synthesis of general distillation-based separation systems. In AIChE Annual Meeting Los Angeles, CA, (Paper 155b). Kataoka, K., Noda, H., Nakaiwa, M., & Horiuchi, K. (2006). Recent energy-saving distillation technology for prevention against global warming. In InterAmerican Congress of Chemical Engineering Buenos Aires. King, C. J. (1980). Separation processes, 2/e. New York: McGraw-Hill. Rathore, R. N. S., van Wormer, K. A., & Powers, G. J. (1974). Synthesis of distillation systems with energy integration. AIChE Journal, 20, 940. Rong, B.-G., Kraslawski, A., & Turunen, I. (2003). Synthesis of functionally distinct thermally coupled configurations for quaternary distillations. Industrial & Engineering Chemistry Research, 42, 1204. Sargent, R. W. H. (1998). A functional approach to process synthesis and its application to distillation systems. Computers & Chemical Engineering, 22, 31. Sargent, R. W. H., & Gaminibandara, K. (1976). Optimum design of plate distillation columns. In L. W. C. Dixon (Ed.), Optimization in action (pp. 267–314). New York: Academic Press. Seader, J. D., & Westerberg, A. W. (1977). A combined heuristic and evolutionary strategy for synthesis of simple separation sequences. AIChE Journal, 23, 951. Tawarmalani, M., & Sahinidis, N. (2002). Convexification and global optimization in continuous and mixed-integer nonlinear programming: Theory, Algorithms, software, and applications. Boston: Kluwer Academic Publishers. Thompson, R. W., & King, C. J. (1972). Systematic synthesis of separation schemes. AIChE Journal, 18(5), 941. Underwood, A. J. V. (1946). Fractional distillation of ternary mixtures. Part II. Journal of the Institute of Petroleum, 32, 614. Westerberg, A. W. (2004). A retrospective on design and process synthesis. Computers & Chemical Engineering, 28, 447. Yeomans, H., & Grossmann, I. E. (1999). A systematic modeling framework of superstructure optimization in process synthesis. Computers & Chemical Engineering, 23, 709.