A Revised Method of Attainable Region Construction Utilizing Rotated

Ind. Eng. Chem. Res. 2010, 49, 10549–10557

10549

A Revised Method of Attainable Region Construction Utilizing Rotated Bounding Hyperplanes David Ming,* Diane Hildebrandt,* and David Glasser School of Chemical and Metallurgical Engineering, UniVersity of the Witwatersrand, 1 Jan Smuts AVenue, Braamfontein 2000, Johannesburg, South Africa

An improved method of candidate attainable region (AR) construction, based on an existing bounding hyperplanes approach, is presented. The method uses a plane rotation about existing extreme points, rather than a plane translation, to eliminate unachievable regions from an initial bounding set. The algorithm is shown to be faster, but has currently only been implemented for two-dimensional constructions and has been extended to include the construction of candidate ARs that involve nonisothermal kinetics in concentration and concentration-time space. 1. Introduction In the past two decades, knowledge of the attainable region (AR) has allowed one to approach the problem of optimal reactor network design from a geometric viewpoint. Originally introduced by F. J. M. Horn in 1967,1 but later pioneered by the work of Glasser et al.,2 Hildebrandt,3 Hildebrandt and Glasser,4 Feinberg and Hildebrandt,5 and Feinberg,6,7 the AR has addressed the problem of optimal reactor design by determining all possible outputs for a given set of feed conditions, reaction kinetics, and operating constraints. In doing so, the AR approach has allowed for the solution of a large variety of reactor network optimization problems. Much of the most recent developments in the field have focused toward the advancement of algorithms for the automatic construction of candidate ARs. Current methods of construction fall into two categories. Methods that attempt to generate the AR from a known starting feed condition and grow the region outward8-10 are referenced as inside-out methods, and those that bound an initial and rather large superset in which the AR is known to reside, and then progressively shrink the region inward,11-13 are referenced here as outside-in methods. In this paper, we shall be concerned with a variation of an existing outside-in algorithm, developed by Abraham and Feinberg,13 for the automated construction of the AR using a large number of bounding hyperplanes. Our chief objective will be to present the necessary changes made to the original and demonstrate its application to problems in concentration and concentration-time space, particularly in R2. We aim to show that, via a change in hyperplane placement, AR construction efficiency may be improved by a significant margin. We shall only be concerned with the construction of candidate ARs in which density and, unless specified otherwise, temperature may be assumed to be constant. The consequences of such choices allow us to assert special geometric characteristics of the process of mixing and, hence, the shape of the AR boundary. Sections 2 and 3 provide the necessary mathematical background, as well as general results required, for the proper understanding of the modified algorithm. These examples allow for the determination of the AR for a wider variety of more-generalized constructions. Sections 4 and 5 introduce the reader to the original and modified construction methods, respectively, where the main * To whom correspondence should be addressed. Tel.: +27 11 717 7557. E-mail addresses: [email protected] (D.M.), [email protected] (D.H.).

motivations for improvement are detailed at the end of section 5. Section 6 gives several examples for classical isothermal problems, as well as nonstandard problems that have not been implemented by a bounding hyperplanes algorithm. Section 6.1 in particular provides a comparison between the two approaches. Finally, section 7 provides several remarks regarding higher dimensional constructions, before concluding remarks are supplied in section 8. 2. Mathematical Preliminaries 2.1. Convex Hulls. We use Rn to denote the space of n-tuples of real numbers. If S ⊂ Rn is a subset such that for any distinct pair of points x,y ∈ Rn, where the line segment joining x and y is completely contained in S, then S is said to be conVex. The conVex hull of a set S ⊂ Rn is the intersection of all the convex sets in Rn that contain S, and this is denoted as conv(S). The convex hull represents the smallest convex set that contains S. Geometrically, the convex hull can be envisioned as a convex polytope in Rn enclosed by a finite number of hyperplanes of which the facets are composed. 2.2. Hyperplanes. A hyperplane H is defined to be one that obeys the following equation: H(C, C0) ) {C ∈ Rn:(C - C0) · n ) 0}

(1)

where vectors C and C0 are points lying in H and where n is a vector perpendicular to all vectors lying in H. We will usually work with unit normals, that is, normal vectors with a magnitude of 1. We define the ( · ) operator as the standard dot product between two vectors, namely, (C - C0) and n. A hyperplane separates a space Rn into two half-spaces. The positiVe and negatiVe closed half spaces are then defined by Hg ) {C ∈ Rn:(C - C0) · n g 0} He ) {C ∈ Rn:(C - C0) · n e 0}

(2)

Similarly, for the open half spaces, H> ) {C ∈ Rn:(C - C0) · n > 0}

H< ) {C ∈ Rn:(C - C0) · n < 0}

(3)

Clearly, n forms an orthogonal subspace to H and may be used as a test for tangency to the plane.

10.1021/ie1004493  2010 American Chemical Society Published on Web 06/29/2010

10550

Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010

2.3. Tangency. Two vectors x and y are tangent, if not all elements in x and y are nonzero and where x·y ) 0

(4)

2.4. Extreme Points. A point x ∈ Rn is an extreme point if it is a vertex of the convex hull. Extreme points may not lie in the interior of the convex hull, nor in the interior of the line segment bounding the facets of the convex polytope. 2.5. Facets and Edges. A hyperplane is said to a n-face or n-facet of the convex polytope conv(S), if n linearly independent points in S lie on H. An n-edge of the convex polytope conv(S) is one that contains n - 1 linearly independent points that comprise the n-face.14 3. The Geometry of Reaction and Mixing 3.1. Reaction. For a prescribed set of reaction kinetics and concentration vector C ∈ Rn, we define the rate of formation of each species i as the function ri(C) describing the instantaneous change of the species as a function of C. The reaction vector r(C) ) [r1(C) r2(C) · · · rn(C)]T, then, is the column vector in Rn that describes the instantaneous change in each component at a point C in concentration space. Hence, for every point in concentration space, we are able to assign a unique rate vector described by r(C) and, consequently, associate a given set of reaction kinetics in terms of a uniquely defined vector field residing in Rn. All vectors used throughout this paper are, unless specified otherwise, assumed to be column vectors. 3.2. Mixing. It is easy to show that the process of mixing two concentrations, C and C*, under the assumption of constant density, results in a final concentration C′ lying between the straight line segment joining C and C*. The mixing vector v ∈ Rn is defined in terms of these two concentrations as v ) C - C*

(5)

and having the property that the final mixture concentration C′ must lie on v. 3.3. Idealized Reactors. Let us investigate the geometric nature of two idealized reactors. We begin with the plug-flow reactor (PFR), which has, as its defining equation, dC ) r(C) dτ

(6)

Here, τ is a non-negative scalar that represents residence time. It is clear from the above discussion that the concentrations mapped out by the solution of eq 6 travel along a trajectory that is tangent to the rate vector at every point along its path. Since r(C) is uniquely defined at C, the trajectories produced by a PFR are uniquely defined by their initial feed concentrations and may never cross each other.3,5 The continuously stirred tank reactor (CSTR) has the unique property that the mixing vector v is collinear with the rate vector evaluated at the reactor exit concentration r(C) and is given by C - C0 ) τr(C)

(7)

Unlike the PFR, the locus of concentrations produced by a CSTR may not be single-valued, because multiple exit CSTR concentrations that satisfy eq 7 may exist for a given feed concentration and residence time. Although the locus is uniquely defined for a given C0, there are infinitely many values of C0 that exist for a given rate vector, which allow CSTR loci to cross.3

3.4. Stoichiometric Constraints. Given a system of chemical reactions and associated kinetics, we aim to determine the limits to which a concentration, for the given set of kinetics, may be achieved by stoichiometry alone. The knowledge of an upper bound on attainability allows one to establish that the AR will lie in a smaller subspace in Rn.7 Clearly, non-negativity constraints restrict the subspace to lie within the positive orthant in Rn; however, it can be shown that even tighter bounds on attainability may be determined through systematic means. The method of determining such limits, in which possibly several reactions and initial feed concentrations are available, is given by Feinberg.6,7 Although the method is not lengthy, it will not be considered in detail here. For a detailed discussion, the authors recommend the papers of Feinberg6,7 and Abraham.15 A brief description of the theory is given below. With n denoting the number of species in a reaction network, a standard basis in Rn is constructed. For each reaction in the network, a vector equation written with respect to the basis vectors is established. The resulting system of linear equations is termed the reaction Vector matrix A, and it is the span of A that represents the stoichiometric subspace S. The rank of A gives the dimension of the stoichiometric subspace and, hence, the dimension in which the AR may be constructed. Knowledge of this subspace forms the basis of the method of bounding hyperplanes,13 which is considered in the following sections. 4. The Method of Bounding Hyperplanes Abraham and Feinberg’s method of bounding hyperplanes is an AR construction technique that utilizes the successive addition of bounding hyperplanes to eliminate unattainable regions from the stoichiometric subspace S. With each successive iteration, unattainable regionssthat is, regions that form part of the stoichiometric subspace but are physically unachievable through any conceivable reaction networksare cut away, resulting in a progressively smaller and tighter bounding set. The hyperplanes are oriented such that division of the two regions results in one of the two halves containing only unachievable concentrations. The orientation of such a hyperplane that agrees with the above condition is done so that the following two requirements are met: • All rate vectors in the open negative half space possess the property that n · r(C) > 0. • The closed positive half space contains the feed point and equilibrium points. Once a hyperplane with the above orientation is achieved, the region of the polytope lying on the negative half space of the hyperplane is discarded and the newly introduced hyperplane is added to the current set of constraints, forming the bounding set. Repeated stages of refinement transform the original stoichiometric subspace into a smaller polytope that contains only those output compositions achievable through reaction and mixing (i.e., the AR). The reason for such a choice of hyperplane orientation satisfying the above two conditions is an important one, because it forms the mechanism by which unattainable regions are eliminated from S. Such a hyperplane orientation is one that may be associated with a steady-state process of a unique arrangement, specifically, one that cannot produce an exit composition residing on the negative half space of the hyperplane. However, the full proof is a rather unnecessary justification in the context of this discussion, and we again direct the reader to the works of Abraham and Feinberg13 and Abraham,15 where a deeper examination of the above claim can be found. We only present a generalized description of the method below.


10551

Figure 1. Unattainable region elimination via rotating hyperplanes about an edge. After each hyperplane addition, the bounded region is smaller than the last and represents a more-accurate approximation of the AR.

The algorithm begins with the initialization of stoichiometric constraints; these constraints are in the form of a list of hyperplanes that bound a region in concentration space. The region bounded by these hyperplanes represents the stoichiometric subspace S. New hyperplanes are then introduced at the corners of S and the trimming process begins. The newly introduced hyperplanes are moved into S until the above conditions are no longer satisfied. At this point, the region has been refined by the addition of additional constraints. The extreme points of the resulting polytope are calculated and the process is repeated. Additional hyperplanes are introduced and the refining process continues until no further hyperplanes can be added without violation of the above conditions to within a specified tolerance. The convex polytope resulting from the bounding process is now a region that represents the set of achievable concentrations through reaction and mixing, and, as a result, it is an approximation to the AR. The method is found to be a robust AR construction technique, and it is particularly successful in constructions that belong to kinetics involving multiple steady states. This success is, above all, a consequence of its approach to infeasible region elimination rather than feasible region addition. What is unclear from the above description though is the large overhead involved for intermediate computations between trimming stages. These stages hamper the total construction time and make the algorithm uncompetitive for the construction of simpler kinetics. A modification to the above technique is presented below, which aims to maintain the virtues of infeasible region elimination while reducing the volume of unnecessary intermediate computations. 5. A Revised Method of AR Construction The original method of bounding hyperplanes seeks to eliminate unattainable regions by introducing a new hyperplane at a corner of the polytope, with a fixed orientation calculated on a weighted average of the hyperplane orientations that shape the corner. The plane is then moved into the current bounding set until stopping criteria are met. In this way, the original method by which hyperplanes are orientated in space may be pictured as a translation through space with fixed orientation. A variation to this approach is proposed, which consequently affect the way in which intermediate computations are arranged. In the same manner as the original, the revised method uses several hyperplanes to successively eliminate regions from an initial bounding set. What differs in the revised approach, however, is the choice of hyperplane orientation and positioning in space. The idea is presented as follows: Instead of fixing the hyperplane orientation and moVing it about through space, fix the hyperplane position and Vary its orientation. This is done through a rotation about an existing edge of the current polytope. (See Figure 1.) Existing extreme points of the polytope generated from previous iterations may then be used as edges

from which new hyperplanes may be introduced and rotated about. New extreme points to the AR may be combined with existing ones obtained from previous iterations, building the polytope face-by-face. The revised method is given as follows: (1) Input the initial stoichiometric constraints in the form of a list of hyperplanes that form the stoichiometric subspace S. These hyperplanes are calculated by the methods described above for the determination of S, and they are specific to the reaction stoichiometry under consideration and feed vector Cf. (2) Determine the hyperplane that passes through the feed point Cf. This is the first hyperplane to be considered for rotation. The feed point will also act as the edge about which the first hyperplane is to be rotated. (3) Using the known extreme point, rotate the chosen hyperplane by a small angle θ. If it is found that any feed or equilibrium points lie outside of the current feasible region due to the new rotation, then unrotate and stop. Output the list of extreme points and terminate the algorithm. (4) Discretize concentrations lying on the plane. (a) For each discretized point C*, compute n · r(C*). (i) If n · r(C*) e 0, then condition 1 is no longer satisfied and a tangent point has been found. Record this concentration C*. (ii) Unrotate the hyperplane back to its previous position. (iii) In addition to this, unrotate the vector (C0 - C*) by angle θ used above. This is done so that the extreme point found in step (a)(i) maintains its position on the newly unrotated hyperplane. Record this rotated C* and go to condition (6). (5) If all discretized points in the list have been considered and it is found that none of these points satisfy the tangency condition, then the current hyperplane orientation may be rotated into the polytope further. Hence, return to condition 3. (6) Combine the new rotated hyperplane with the current list of hyperplanes currently bounding the region. Add the newest extreme point found in condition (4)(a)(iii) to the current list of extreme points. Then go to condition (3). Rotations are achieved through the use of a rotation matrix R ⊂ Rn. The direction of rotation is chosen such that after each rotation, the resulting polytope is smaller than the last. Rotations performed in the special two-dimensional setting are easy to compute, involving a single matrix multiplication. In addition to this, there is only one plane through which rotations may occur, and the edges rotated about are single points. The benefit afforded to us from rotations, as opposed to translations through space, is subtle, and only arises upon implementation of the two methods. Further examination reveals that the original undergoes two computationally intensive, yet entirely necessary, stages. The first, which is common and unfortunately necessary to both, is the act of hyperplane discretization. Since a hyperplane represents an (n - 1)dimensional subspace, hyperplane discretization may be viewed as a computational process of degree n - 1. Hence, by adopting

10552


Big-O notation, we have that the number of floating point operations (flops) for such a process is proportional to approximately O(Mn-1), where M is the number of discretized points for each dimension and n is the dimension of the computational space (the dimension of S). The second intensive computation is altogether suppressed during the early stages of refinement but quickly emerges as the number of hyperplanes bounding the current polytope increases throughout construction. For new hyperplanes to be introduced, we must know the locations of the extreme points of the polytope defined by the current list of hyperplanes. The question arises how we may achieve this in the most efficient manner possible. We may, rather naively, attempt to find the intersection of all planes with all other planes in space and check that the resulting point of intersection is satisfied by all current constraints. However, this is not recommended, because L hyperplanes in Rn result in, at most, L!/[n!(L - n)!] total intersections. The problem is termed the vertex enumeration problem and has long been studied in the areas of economics, computational geometry, and game theory.16-18 Methods have been established that solve the problem in polynomial time;19 however, it is evident that the necessary stage of vertex enumeration is an unavoidable expense. A final consideration is the issue of redundant hyperplanessthat is, hyperplanes that play no useful role in defining the feasible region. We may again, for example, choose to ignore these hyperplanes as their existence in space offer no contribution to the bounding process. However, it is wise to exclude them as quickly as possible, because their numbers can quickly swamp the number of useful planes that describe the feasible region, crippling the vertex enumeration computation in the process. As a result, the additional expense of redundant hyperplane removal contributes further to the total computation time of the bounding process. The revised method utilizing rotations does not require the enumeration of vertices, and, as a result, the costly exercise of redundant hyperplane removal and vertex enumeration stages may be avoided. This is accomplished by the addition policy of the revised method. Since rotations are done about existing, known attainable points, the hyperplane tangency condition guarantees that the new extreme point found (in R2) is attainable by CSTR from eq 7. Hence, at every stage of construction, boundary points of the AR are found, and intermediate points never enter the problem. Tangent points to the plane, which are found for constructions in higher dimensions, generally may not be satisfied by eq 7 and, hence, intermediate points may be found; however, the enumeration of vertices is still an unnecessary computation. We refer the reader to section 7 for a discussion on how constructions may be performed in higher dimensions. 6. Examples 6.1. Isothermal Constructions: Comparison with the Original Method. We now demonstrate how the method of bounding hyperplanes may be improved by plane rotations for various examples. In all the examples that are presented, we have assumed that the systems are operated isothermally and under constant density. Although we shall later present an example in which the isothermal assumption is relaxed, we honor this restriction for the purposes of comparison to the original method. The assumption of constant density allows for the claim that mixing is necessarily a convex process. Although this may be found to be most suitable for isothermal liquid phase reactions, it is merely used for the convenience of illustration. However, reactions such as those which occur in the gas phase

Table 1. Rate Constants for the Van de Vusse Reaction parameter

value

k1 k2 k3 k4

1.0 0.0 1.0 20.0

such that

r(C) )

[

-k1cA + k2cB - 2k4cA2 k1cA - k2cB - k3cB

]

are clearly those in which the assumption will rarely be honored and, thus, the act of mixing may not be convex in general. Such a limitation may be overcome if, instead, we recast the equivalent problem in a space in which mixing is necessarily convex. Therefore, for problems where constant density may not be suitable, we choose to work with mass fractions instead,20,21 where the parallel between mixing using units of mass is easily proved to always obey linear mixing laws. We begin with the well-understood Van de Vusse reaction scheme, which is given by the following system: k1

k3

A {\} B 98 C k2 k4

(8)

2A 98 D Several variations of the above kinetics exist, but only the simplest scenario is presented here: a two-dimensional (2D) system in cA-cB space with mass action kinetics given by the constants shown in Table 1. The AR is then found to be produced by a CSTR followed in series by a PFR. Unfortunately, this arrangement cannot be established by the algorithm, but it has been included to demonstrate the accuracy to which constructions may agree with their theoretical predictions. The results of the construction are given by Figure 2 below, where a comparison is made between the original (translated) and revised (rotated) methods. The total running time for the original method is 10 s and is bounded by 32 hyperplanes. The shape of the boundary produced by the original method agrees with the analytical solution, although a distinct difference between the two is noticeable near the CSTR-PFR junction. In particular, the maximum concentration of component B is overestimated slightly with the original method. The results achieved with the revised approach, by comparison, are in close agreement with the theory. The time required to construct the new region is 6 s and is constructed with 116 hyperplanes. Although running time is of similar order to the former, the region is assembled with a significantly larger number of hyperplanes, allowing for a more precise bound. The results highlight an important consideration associated with both methods: that is, there is no claim that the resulting polytope obtained from construction is indeed the true AR. Rather, the results provide an upper bound on kinetic attainability. Certainly, there is evidence to suggest that, in the limit, as more hyperplanes are used to bound the region, we approach the true AR boundary; however, the regions obtained by either method are necessarily always larger than the true solution. If both an inside-out and outside-in method are available, then comparison between the two candidate ARs may allow for the determination of the true AR boundary. However, if the results of the two methods do not correspond, then a closer bound on


10553

Figure 2. Comparison of construction methods performed for Van de Vusse kinetics. The CSTR locus (denoted by crosses, ×) and PFR trajectories (denoted by circles, O) have also been drawn.

Figure 3. Nonlinear Van de Vusse kinetics.

the true boundary is reached, because it must lie between the constructions. Next, we study a fairly artificial example: one in which highly nonlinear behavior is observed. This may be produced by the following rates of formation: rA(C) ) f(cA)2 + cAcB rB(C) ) cAf(cA)

df(cA) dcA

where f(cA) ) 6cA6 - 6cA5 + 9cA4 - 16cA3 + 9cA2 - 2cA The reaction scheme is identical to that posed for regular Van de Vusse kinetics by eq 8, although the actual rate expressions are more complicated now. A PFR trajectory associated with the kinetics given above show multiple concavities, indicating that both reaction and mixing processes will feature prominently in the formation of the AR boundary. Indeed, it is found that, for a feed concentration of Cf ) [1 0]T, the AR is constructed from a four-reactor series network of CSTR-PFR-CSTRPFR.2 Figure 3 provides results obtained from both methods of construction. Once again, both the original and revised methods allow for an accurate approximation of the AR boundary, even for highly nonlinear kinetics; however, the efficiency by which constructions are performed differ greatly between the two. To construct the AR via the original method in a time of 30 s, 77 hyperplanes

have been used; this result agrees favorably with the results obtained from the new approach. However, the latter is produced with more than three times as many hyperplanes and is accomplished within a significantly shorter time. Hence, in a period of 13 s, the method of rotations has allowed for the introduction of 262 hyperplanes, producing a tightly bound region. For this example, the use of either method requires the introduction of a rather large number of hyperplanes for an adequate approximation of the AR to be made. Multiple curved and straight line segments are noticeable on the boundary, and it is apparent how curvature is approximated with many hyperplanes. Hence, it can be expected that the most gains in calculation time may be associated with boundaries having high curvature. We next examine a situation in which this is not the case. Hence, we now compare constructions for an example where multiple steady states allow for the occurrence of discontinuous regions to be observed from within the stoichiometric subspace. The reaction system is given by the following rates of formation:

[

rA(C) ) -k1 cA +

[

rB(C) ) k1 cA +

b(cAcB)2 1 + k2cAcB2 b(cAcB)2

1 + k2cAcB2

]

]

+ acA

with k1 ) 1 min -1, k2 ) 40 L3/mol3, a ) 100, and b ) 104 L3/mol3. We shall again assume that a single feed stream exists, with pure A, such that the feed stream concentration vector is Cf )

10554


Figure 4. AR construction with multiple steady states. Crosses (×) represent the CSTR locus for the feed concentration.

[1 0]T. A plot of the CSTR locus identifies two separate branches that exist for the specified feed concentration. It is possible, using standard techniques,2 to show that a CSTR followed by a PFR is required to obtain the AR associated with the kinetics provided. The results are given by Figure 4. The candidate AR found indicates a maximum achievable cB concentration of ∼0.61. For comparison, the CSTR locus of X values corresponding to the feed concentration has been plotted. The second branch is an isola, and, from the constructions, the significant contribution that it provides in the full range of total achievable concentrations is apparent. Candidate AR construction via inside-out methods requires knowledge of multiple steady states, as the discovery of only partial convex regions may result if knowledge of their existence is unknown. However, prior knowledge of such behavior is often difficult to predict. Indeed, behavior in isothermal reactions may vary considerably, even for noticeably similar kinetics,22 and occur frequently in adiabatic reactions.3 Both methods produce an accurate bound on the set of all attainable concentrations, despite the fact that there exist irregular areas formed by the isola. The AR exhibits minimal curvature, that is, the PFR trajectories vary somewhat linearly near the AR boundary, which allow for a rather simply shaped triangular region that may be approximated well with few hyperplanes. As a result, the construction times and accuracies produced by both methods are fairly competitive and little benefit is obtained with the use of plane rotations rather than translations. Consequently, the AR constructed by the original approach is approximated by 16 hyperplanes and was achieved within a time of 16 s while the revised method bounds the same region in a time of 4 s with 18 hyperplanes. Nevertheless, we notice that, even for kinetics in which the AR may be approximated well with few hyperplanes, opportunities for improvement still exist. Although such examples exhibit minimal overhead involved in vertex enumeration, we find that rotations allow for more-appropriate placement of hyperplanes in space and offer a more efficient means of candidate AR construction. 6.2. Temperature-Dependent Kinetics. Let us know consider an example in which the isothermal constraint is no longer required. The system of reactions we consider here again take the form of the familiar Van de Vusse kinetics; however, we shall now allow the rate constants to be of the Arrhenius form23

( )

ki(T) ) k0i exp -

Ei RT

(9)

The values for ki0 and Ei for each component are given in Table 2.

Table 2. Arrhenius Constants for Temperature-Dependent Van de Vusse Kinetics i

ki0

(Ei)/R (K)

1 2 3 4

4 1.5 6 0

500 800 0 0

The rate vector is no longer a sole function of C, but takes the form r(C, T) )

[

-k1(T)cA - 2k4(T)cA2 k1(T)cA - k3(T)cB

]

where it is now possible to vary the direction in which the rate vector points by a variation in temperature. Hence, for every point in concentration space, the rate vector is allowed to take on multiple directions, depending on the temperature at which it is evaluated. The handling of temperature-dependent kinetics is a fairly straightforward modification to the standard approach discussed in section 6.1. Instead of checking for a single tangent rate vector at a concentration residing on the plane, we now check for a range of rate vectors. This range corresponds to the range of allowable temperatures specified between given operating temperature limits. Much in the same way that a grid is generated in concentration space for calculation of discretized points on the plane, points in temperature space between the maximum and minimum permissible temperature limits are now also determined. Then, for each valid concentration generated on the plane, the range of rate vectors between the operating limits are all evaluated for tangency. Figure 5 shows the results of the construction for the nonisothermal case operating between a temperature interval of 300-1000 K. To maintain a path on the boundary of the AR, the temperature profile would need to follow the one suggested. It can be seen that the temperature profile declines rapidly over a small cA concentration range from ∼0.25 to 0.18. The optimum recommended operating temperature profile is in agreement with the recommendation suggested by Goddor et al.23 Significantly more construction time was required for this case, totaling 479 s, which is largely due to a fine grid of 500 points used between the temperature interval and small rotation angle of 1 × 10-4 radians. The boundary is composed of 150 hyperplanes. AR boundary outlines corresponding to isothermal operation at the minimum and maximum allowed operating temperature are also provided in Figure 5. Hence, the AR resulting from a design in which the operating temperature is maintained at 1000


Figure 5. Temperature-dependent kinetics. The ARs constructed for isothermal operation at the two operating limits are both smaller than that obtained by the optimal temperature profile.

K allows for a large portion of the total optimal set of concentrations to be achieved. This is expected because the optimal temperature profile suggests that, for the cA range of ∼1.0-0.25, the temperature that extends the region the most is one that is maintained at its maximum value. A design using only the lower temperature of 300 K, on the other hand, results in a region that is noticeably smaller between the entire cA range. What is apparent, nevertheless, is that isothermal designs maintained at either temperature do not result in a candidate region that is considered optimal, because a small range of concentrations near cA ) 0.1 to cA ) 0 exists that is only achieved by a varying temperature profile. Determination of the optimal profile is easily handled by this approach, because no other modifications to the underlying elimination method need to be changed. Constructions for nonisothermal kinetics allow for inclusion of many more reactor optimization problems, and ones which are much more likely to occur in reality when used with mass fractions described above. 6.3. Unbounded Construction. Finally, we consider finding the AR for the case where we wish to determine the smallest reactor volume; such examples are common in the design of batch processes.24 Therefore, we consider an AR construction in concentration-residence-time space, where the resulting feasible region is unbounded. The reaction is assumed to be an adiabatic decomposition under constant pressure, where the rate of formation is again assumed to follow the Arrhenius form

(

r(c) ) ck1 exp -

4000 8000 - (1 - c)k2 exp T(c) T(c)

)

(

)

where c is concentration, k1 ) 5 × 105, and k2 ) 5 × 108. For the sake of simplicity, we assume a constant specific heat and heat of reaction; the temperature for the process is then T(c) ) T0b + Tad(1 - c) where Tb0 ) 300 K and Tad ) 200 K. In this space, the rate vector is given as r(c) )

[ ] r(c) 1

We specify that the feed stream does not contain any products so that c ) 1 and construct the unbounded region in the same manner as the previous examples. The only additional alteration required is the direction of rotation. Given the fact that equilibrium is achieved as the residence time is increased, we now specify that hyperplanes must be rotated in a clockwise direction for the region to be bounded correctly. The results of

10555

Figure 6. Unbounded AR construction.

the construction are provided by Figure 6. For the specified rate, a minimum residence time of ∼0.09, at which point the minimum cA value was 0.16, has been reached. The region is defined by 78 hyperplanes, which took 43 s to construct. Construction of the AR face by face presents unique opportunities for the construction of a special class of convex polytope: unbounded polytopes. Elimination via a rotation implies no dependence or knowledge of all existing corners to the current polytope; hence, reliance on a closed polytope is not a requirement. Constructions of these types are unique to the revised method and are unachievable via the regular method. Both the original and revised bounding hyperplane algorithms were programmed in C++ using Microsoft Visual Studio 2008. The computations for all examples discussed were performed on an Intel Core i3-530 processor (2.93 GHz) and 4 GB RAM (DDR3 1600) under Windows 7 64-bit Edition software. 7. Remarks on Higher-Dimensional Construction Although the discussions in the previous sections have served well in discussing the workings of the modified approach, they have remained in R2. Constructions in Rn are currently not yet fully understood unfortunately; however, we provide the reader with several ideas of how they may be conducted and indicate, where necessary, the complications involved in generalized n-dimensional constructions. We begin with a basic understanding of the facets and edges associated with an n-dimensional convex polytope, hereby denoted by conv(S). By section 2, we have an n-face or n-facet of conv(S) that must be a collection of extreme points such that n linearly independent points in conv(S) lie on a hyperplane H residing in Rn. However, the choice of extreme points on H is not arbitrary, because we require that H be a support hyperplane to conv(S). Edges of the polytope are understood in much the same way as facets, although for edges, we require that n - 1 linearly independent points in an n-face lie on a support hyperplane. Rotations in Rn may be performed in an identical manner to that in R2, that is by use of an n × n rotation matrix Rn. However, the method of rotations is complicated by the existence of n!/[2!(n - 2)!] principal axes about which bounding planes may be rotated.25 Hence, arbitrary rotations in space may be achieved by linear combinations along the principal rotation axes. The rotation axis may not take on any arbitrary orientation though; instead, they must be aligned with edges of existing facets. In doing so, the method remains conceptually the same when extended to higher dimensions. Hence, we intend to enclose the AR within a rather large region P0, determined by stoichiometric constraints, and then seek to successively eliminate unattainable regions by the further addition of hyperplanes brought about by rotations. P0 is reduced to Pr by r trimming

10556


stages and, in this regard, Pr is a closer bound on all attainable concentrations. Once again, construction of Pr occurs facet-byfacet, utilizing existing n-faces and n-edges about which new hyperplanes may be introduced and rotated from. Generalization of the method to Rn may be achieved as follows. Construction begins by determination of the stoichiometric subspace S in the usual manner, producing the first and largest convex polytope belonging to the kinetics under consideration. The extreme points belonging to S are determined so that the facets and edges associated with it may be constructed. The orientation of a face in Rn provides an initial orientation of the new bounding plane, while the edges of the face provide an axis about which rotations may be achieved. Hence, additional hyperplanes are aligned such that all n points of a particular face lie in the plane and are rotated about an axis that belongs to an edge of that face. For a problem residing in Rd then, d points must be known before a facet may be determined, from which d possible edges belonging to the face may be chosen for rotations. The first elimination stage begins by identifying a face that passes through the feed point Cf and is used to orientate the first bounding hyperplane. The hyperplane is rotated until a tangent rate vector lying in H is found, at which point the position in space is recorded and a new face is constructed from the remaining d - 1 points (which comprise the rotation edge) and new tangent point. Construction of a new face establishes an orientation for the next hyperplane, while a new rotation edge is obtained from the establishment of further edges belonging to the new face. After each successive elimination step, a new face of a smaller convex polytope is determined, and construction continues until no further rotations are achievable to within a specified tolerance. Initially, points specified in the feed are the only ones that are known to be achievable and construction of the first face is determined by computing the initial extreme points of S. However, additional tangent points that are found from rotations generally may not be satisfied by eq 7 and, consequently, may not be achievable either. Nevertheless, we continue in the usual manner and append additional facets and edges to the current polytope from successive rotation steps. After each step, the introduction of new points will serve to break ties with the original extreme points. Eventually, a stage in construction is reached where the edges of the current face are no longer composed of the extreme points of the original face. A hyperplane may then be swept over the original extreme points where the test for tangency, with respect to the new plane, may be applied. Regardless of whether the point is achievable or not, a region may only be eliminated provided all rate vectors residing on the open negative half space satisfy the condition n · r(C) > 0 (specified in section 4), with respect to all hyperplane orientations. Hence, if the rate vector evaluated at the point fails to be tangent to the new plane, the region is not physically realizable with any reactor configuration and may be discarded from the set. Eliminations performed in this way may seem rather conservative when viewed in this regard and may not guarantee that tangent points will be achievable, although regions that are discarded are guaranteed to be unachievable. Uncertainty still remains as to whether these eliminations do indeed converge to the true AR (such questions presently are an active area of discussion), although the polytope resulting from hyperplane bounding will always enclose the AR and always be smaller than conv(S). Hence, a rather improved approximation on the limits of attainability may be determined from its construction.

The special setting of R2 allows construction to occur on an ordered (Hamiltonian) path, where each extreme point discovered is visited exactly once. However, higher-dimensional constructions are not generally ordered along such a path, and backtracking to previously found extreme points via the steady addition of new rotations is expected. The preceding statements may concern the reader regarding the benefit obtained with rotations as an efficient means of higher-dimensional constructions, given the fact that there are overheads involved in the bookkeeping of facets and, more importantly, the need for temporary extreme point elimination via backtracking and the uncertainty of achievability of further tangent points that may be found. However, the authors believe that a rather large savings associated with the revised method may be achieved because there is no need for a vertex enumeration step. It has been shown13,15 that the number of hyperplanes required for adequate AR approximation increases rapidly with increasing dimension; hence, the act of vertex enumeration may be a limiting step in the determination of higher-dimensional constructions. 8. Conclusion Originally posed for bounded isothermal constructions in concentration space, Abraham and Feinberg’s method of bounding hyperplanes has proven to be a robust method of AR construction, particularly for the determination of degenerate kinetics such as that in which multiple CSTR steady states occur. However, constructions performed in this way are hindered by computational complexities, which become ever more apparent with an increasing number of hyperplanes. It is found that the cost of hyperplane discretization, extreme point enumeration, and redundant hyperplane removal are computationally demanding stages in the original method that otherwise hinder a rather novel construction technique. A revised algorithm using much of the positive virtues linked with original has been presented and shown to construct regions that agree with the constructions of other authors.2,3,5,10,13,23 In much the same way as its predecessor, the method relies on the successive addition of bounding hyperplanes to separate a region (containing both achievable and unachievable parts) into two half spaces, such that one-half contains only unattainable concentrations. Orientation of the cutting plane has been revised to allow for a rotation about an axis, rather than a fixed orientation and translation through space. This allows for the same trimming mechanism found desirable in the original approach, while simultaneously eliminating the need for vertex enumeration and redundant hyperplane removal after the introduction of each elimination stage. The method demonstrates an improvement in both running time and construction accuracy, when compared to the original for the same problem; however, the most gains are observed for kinetics where the AR boundary is composed of many hyperplanes. These situations involve a large overhead in costly vertex enumeration stages. Construction via rotations of a plane has allowed for the approximation of unbounded ARs in concentration-time space as well as extended to allow for determination of nonisothermal ARs, and it has been shown to agree with the regions obtained via traditional methods. Extension of the method for the inclusion of such cases has allowed for the determination of a rather broader set of reactor network optimization problems, which occur frequently in many common reactor designs. Although the method has been implemented in two dimensions currently, a discussion on extension to n dimensions has been presented. The ideas that characterize construction in R2

Ind. Eng. Chem. Res., Vol. 49, No. 21, 2010 n

take on familiar roles in R . Extreme points of existing facets provide starting orientations (positions in space), while the edges of the same facet provide an edge about which the bounding plane may be rotated. Elimination performed in this way allows for the creation of new facets and edges about which new hyperplanes may be introduced and rotated to further advance the elimination process. Unfortunately, higher-dimensional constructions are not as simple to perform as those determined in two dimensions, and discussions about how these methods can be improved are still actively being investigated. This is due, in part, to the need for additional bookkeeping of facets and edges; however, it is chiefly due to tangent points, which are not generally guaranteed to be achievable. Acknowledgment The financial assistance of the National Research Foundation (NRF) toward this research is hereby acknowledged. Opinions expressed, and conclusions arrived at, are those of the author and should not necessarily be attributed to the NRF. Literature Cited (1) Horn, F. Attainable and Non-Attainable Regions in Chemical Reaction Technique. In Proceedings of the Third European Symposium, Chemical Reaction Engineering: Pergamon: London, 1964. (2) Glasser, D.; Hildebrandt, D.; Crowe, C. A geometric approach to steady flow reactors: The attainable region and optimization in concentration space. Ind. Eng. Chem. Res. 1987, 26, 1803–1810. (3) Hildebrandt, D. The Attainable Region Generated by Reaction and Mixing, Ph.D. Thesis, University of the Witwatersrand, Johannesburg, South Africa, October 1989. (4) Hildebrandt, D.; Glasser, D. The attainable region and optimal reactor structures. Chem. Eng. Sci. 1990, 45, 2161–2168. (5) Feinberg, M.; Hildebrandt, D. Optimal reactor design from a geometric viewpointsI. Universal properties of the attainable region. Chem. Eng. Sci. 1997, 52, 1637–1665. (6) Feinberg, M. Optimal reactor design from a geometric viewpoint. Part II. Critical sidestream reactors. Chem. Eng. Sci. 2000, 55, 2455–2479. (7) Feinberg, M. Optimal reactor design from a geometric viewpointsIII. Critical CFSTRs. Chem. Eng. Sci. 2000, 55, 3553–3565. (8) Rooney, W. C.; Hausberger, B. P.; Biegler, L. T.; Glasser, D. Convex attainable region projections for reactor network synthesis. Comput. Chem. Eng. 2000, 24, 225–229. (9) Kauchali, S.; Rooney, W. C.; Biegler, L. T.; Glasser, D.; Hildebrandt, D. Linear programming formulations for attainable region analysis. Chem. Eng. Sci. 2002, 57, 2015–2028.

10557

(10) Seodigeng, T.; Hausberger, B.; Hildebrandt, D.; Glasser, D. Recursive constant control policy algorithm for attainable regions analysis. Comput. Chem. Eng. 2009, 33, 309–320. (11) Burri, J. F.; Wilson, S. D.; Manousiouthakis, V. I. Infinite DimEnsionAl State-space approach to reactor network synthesis: Application to attainable region construction. Comput. Chem. Eng. 2002, 26, 849–862. (12) Manousiouthakis, V. I.; Justanieah, A. M.; Taylor, L. A. The ShrinkWrap algorithm for the construction of the attainable region: An application of the IDEAS framework. Comput. Chem. Eng. 2004, 28, 1563–1575. (13) Abraham, T. K.; Feinberg, M. Kinetic Bounds on Attainability in the Reactor Synthesis Problem. Ind. Eng. Chem. Res. 2004, 43, 449–457. (14) Chand, D. R.; Kapur, S. S. Algorithm for convex polytopes. J. Assoc. Comput. Mach. 1970, 17, 78–86. (15) Abraham, T. K. Kinetic bounds on attainability in the reactor synthesis problem, Ph.D. Dissertation, Ohio State University, Columbus, OH, 2005. (16) Balinski, M. An algorithm for finding all vertices of convex polyhedral sets. J. Soc. Ind. Appl. Math. 1961, 9, 72–88. (17) Matheiss, T. H.; Rubin, D. S. Survey and comparison of methods for finding all vertices of convex polyhedral sets. Math. Oper. Res. 1980, 5, 167–185. (18) Dyer, M. E.; Proll, L. G. An algorithm for determining all extreme points of a convex polytope. Math. Program. 1977, 12, 81–96. (19) Dyer, M. E. Complexity of vertex enumeration methods. Math. Oper. Res. 1983, 8, 381–402. (20) Hildebrandt, D.; Glasser, D.; Crowe, C. M. Geometry of the attainable region generated by reaction and mixing: With and without constraints. Ind. Eng. Chem. Res. 1990, 29, 49–58. (21) Zhou, W.; Manousiouthakis, V. I. Variable density fluid reactor network synthesissConstruction of the attainable region through the IDEAS approach. Chem. Eng. J. 2007, 129, 91–103. (22) Schlosser, P. M.; Feinberg, M. A theory of multiple steady states in isothermal homogeneous CFSTRs with many reactions. Chem. Eng. Sci. 1994, 49, 1749–1767. (23) Godorr, S.; Hildebrandt, D.; Glasser, D.; McGregor, C. Choosing Optimal Control Policies Using the Attainable Region Approach. Ind. Eng. Chem. Res. 1999, 38, 639–651. (24) Nicol, W.; Hernier, M.; Hildebrandt, D.; Glasser, D. The attainable region and process synthesis: Reaction systems with external cooling and heating. The effect of relative cost of reactor volume to heat exchange area on the optimal process layout. Chem. Eng. Sci. 2001, 56, 173–191. (25) Aguilera, A.; Pe´rez-Aguila, R. General n-Dimensional Rotations. Presented at the 12th International Conference in Central Europe on Computer Graphics, Visualization and Computer Vision, Plzen-Bory, Czech Republic, February 2-6, 2004.

ReceiVed for reView March 1, 2010 ReVised manuscript receiVed June 8, 2010 Accepted June 10, 2010 IE1004493

A Revised Method of Attainable Region Construction Utilizing Rotated

Recommend Documents