Ind. Eng. Chem. Res. 1992,31, 300-312
300
Lopez, A. M.; Murrill, P. W.; Smith, C. L. Controller Tuning Relationships Based on Integral Performance Criteria. Instrum. Technol. 1967,14 (ll),57. Nijmeijer, H.; Schumacher, J. M. The Regular Non-Interacting Control Problem for Nonlinear Control Systems. SZAM J. Control Optim. 1986,24,1232. Pandit, H. G.; Rhinehart, R. R. Process Model-Based Control of a Fluidized Bed Reactor: A Comparison of Two Strategies. In Advances in Instrumentation; Instrument Society of America: Research Triangle Park, NC, 1989;Vol. 44,p 703. Pandit, H. G.; Rhinehart, R. R.; Riggs, J. B. Process Model-Based Control of a Nonideal Binary Distillation Column. Presented at the 1990 AIChE Annual Meeting, Chicago, IL, 1990. Peng, D.-Y.; Robinson, D. B. A New Two-Constant Equation of State. Znd. Eng. Chem. Fundam. 1976,15,59. Ramchandran, B.; Riggs, J. B.; Heichelheim, H. R.; Seibert, A. F.; Fair, J. R. Dynamic Simulation of a Supercritical Fluid Extraction Process. Presented at the 1990 AIChE Annual Meeting, Chicago, IL, 1990. Rhinehart, R. R.; Riggs, J. B. Process Control through Nonlinear Modeling. Control 1990,1 (7),86-90. Riggs, J. B. Nonlinear Process Model-Baaed Control of a Column with a Sidestream Draw-off. Presented at the 1989 AIChE Annual Meeting, San Francisco, CA, 1989. Riggs, J. B. Nonlinear Process Model-Baaed Control of a Propylene Sidestream Draw Column. Znd. Eng. Chem. Res. 1990,29,2221. Riggs, J. B.;Rhinehart, R. R. A Review of Nonlinear Approximate Models for Process Control and Online Optimization. Paper presented at the 1989 AIChE Annual Meeting at San Francisco, CA, 1989.
Riggs, J. B.; Rhinehart, R. R. Comparison between two Nonlinear Process Model-Based Controllers. Comput. Chem. Eng. 1990,14 (lo), 1075-1081. Riggs, J. B.; Watts, J.; Beauford, M. Advanced Model-Baaed Control for Distillation. Presented at the NPRA Computer Conference, Seattle, WA, Oct 1990. Smith, B. D.;Brinkley, W. K. General Short-cut Equation for Equilibrium Stage Processes. AZChE J. 1960a,6,446. Smith, B. D.; Brinkley, W. K. Supplement to the General Short-cut Equation for Equilibrium Stage Processes. Document 6352;American Document Institute, Photoduplication Service, Library of Congress: Washington 25,DC, 1960b. Smith, C. A.; Corripio, A. B. Principles and Practice of Automatic Process Control; John Wiley & Sons: New York, 1985; p 226. Tzouanas, V. K.;Luyben, W. L.; Georgakia, C.; Ungar, L. H. Expert Multivariable Control. 1. Structure and Design Methodology. Znd. Eng. Chem. Res. 1990,29,382. Van Winkle, M.Distillation; McGraw-Hill Co.: New York, 1967;pp 533-542. Williams, G. L.; Rhinehart, R. R.; Riggs, J. B. Inline Process Model-Based Control of Wastewater pH Using Dual Base Injection. Znd. Eng. Chem. Res. 1990,29,1264. Zheng, D.-L. Modeling and Optimization of the Batch Free-radical Polystyrene Polymerization Process. M.S. Thesis,Department of Chemical Engineering, Texas Tech University, Lubbock, TX, 1990. Received for reuiew December 11, 1990 Revised manuscript received August 5, 1991 Accepted August 15, 1991
Constructive Targeting Approaches for the Synthesis of Chemical Reactor Networks Subash Balakrishna and Lorenz T. Biegler* Department of Chemical Engineering, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213
The impact of the reaction system and reactor design on the character of the flowsheet has been appreciated by many researchers in process synthesis. In this paper, we develop bounds, or targets, on a given performance index of an isothermal reacting system. A targeting model based on optimizii flows between different reacting environments is described, which provides a sufficiently rich representation of the reactor network. The model is formulated as a dynamic optimization problem, and solution strategies are presented. A linear programming simplification, which corresponds to the segregated flow limit, is derived for the targeting model. Sufficient conditions for global optimality of this simplification are described. Moreover, even when these conditions do not apply, testa formulated as nonlinear programs can be applied to evaluate the suitability and enhancement of the target obtained by the segregated flow model. This approach leads to an iterative process for the targeting problem, based on two alternative methods. Also, the reactor network is straightforward to realize from the targeting solution. The approach is demonstrated on numerous isothermal reaction mechanisms and is compared to literature results. The extension of this approach to embedding the reactor target within a process flowsheet is also presented. 1. Introduction
Process synthesis has seen a great deal of progress since its inception, and the systematic synthesis of subsystems such as heat-exchanger networks and separation schemes has marked some of the major milestones. Here, physical insight and heuristics provide a firm basis for targeting, or bounding, good solutions to these problems. However, few results applied to reactor network synthesis are available, mainly because the reactions themselves are poorly understood and rate laws are hard to determine. The interplay of heat- and mass-transfer effects further
* Author to whom correspondence should be addressed.
complicates the situation. Even with homogeneous reacting systems, which we assume here, the problem is usually described by highly nonlinear differential algebraic equation models. Hence, heuristics and intuitive techniques can deal only with limited reaction systems. Despite these complexities, reactor network synthesis is a key task, because the reaction system usually determines the entire character of the process (Douglas, 1989). Among the recent attempts toward reactor network synthesis,superstructurebased approaches have played a very important role. Chitra and Govind (1985) optimized a structure consisting of recycle reactors in series. By placing heat exchangers at reactor inlets, they also allowed for interstage cooling or heating and treated yield-based ob-
0888-588519212631-0300$03.00/0 0 1992 American Chemical Society
Ind. Eng. Chem. Res., Vol. 31,No. 1, 1992 301 jective functions. Achenie and Biegler (1986)used a general superstructure of axial dispersion reactors in a series/parallel connection allowing for component splits and any objective function. Later, Achenie and Biegler (1990)also considered a general superstructure of recycle reactors with bypasses. The synthesis problem was converted to a nonlinear program and solved by deriving adjoints for a gradient-based optimization procedure. Kokossis and Floudas (1990)developed a mixed-integer nonlinear programming (MINLP) model to optimize a very general superstructure of isothermally operating plug flow reactors (PFR) and continuous stirred tank reactors (CSTR). By treating a PFR as the limit of a large number of CSTR units,they were able to keep the formulation free of differential equations. The superstructure-based approaches have the advantage of simultaneously generating the network along with the solution of the optimization problem. A more complete review of these is provided by Chitra and Govind (1985)and Kokossis and Floudas (1990). However, here the solution is only as rich as the superstructure chosen, and the nonuniqueness in reactor networks can play a very important role in superstructure-based approaches. For example, for isothermal systems, product concentrations attainable by two PFRs in parallel can also be attained by a single PFR with a sidestream. Targeting for reactor networks (i.e., deriving realistic bounds) was first suggested by Horn (1964)through the concept of attainable regions in concentration space. However, specific procedures for these are a relatively new development. Glasser and co-workers (Glasser et al., 1987; Hildebrandt et al., 1989;Hildebrandt and Glasser, 1990) attempt to find the entire region in concentration space that is attainable from some base concentration, by representing mixing and reaction using geometric concepts. They derived a set of necessary conditions for any region in space to qualify to be an attainable region. By drawing alternate PFR and CSTR trajectories, they were able to map the entire attainable region and develop a constructive technique for deriving the reactor network. The geometric concepts provide a great deal of insight into the reaction process and lead to some very powerful results. However, their method, though elegant, faces a dimensionality problem; geometric techniques are difficult to apply beyond three dimensions. Also, the necessary conditions have to be checked at infinitely many points, a procesa that may be very cumbersome for multicomponent problems. Achenie and Biegler (1988),on the other hand, posed the reactor synthesis problem as the following two-part problem: (a)Find a bound on the performance index (targeting) irrespective of the reactor type and configuration. (b)Find the network of reactor types that will achieve this target index. For part a, they modified and extended a two-compartment mixing model by Ng and Rippin (1965). By postulating a variable residence time distribution and mixing profiles between two compartments of segregation and maximum mixedness, they were able to simulate a variety of mixing states and were successful in generating targets for many isothermal/nonisothermal reacting systems. Using the superstructure-based approach discussed earlier, they were able to realize the network. In this study, we develop a simpler and more efficient formulation for reactor targeting by considering special cases of the Achenie and Biegler model (1988)and by applying some of the attainable region concepts of Glasser et al. (1987).The resulting procedure consists of solving
W l
AAA
Xexit
I 1 1 1 1
Figure 1. Three-compartment mixing model.
small linear and nonlinear programs and leads to a straightforward and constructive technique for generating optimal reactor networks. This technique is illustrated on a battery of example problems on reactor optimization from the literature, included among which is an example of integrated optimization of the reactor within a flowsheet. Practical results show significant improvements in performance indexes and thus indicate a lot of promise for this approach toward reactor network synthesis. In the next section, we extend the isothermal targeting model to multiple regions of maximum mixedness. This is a general formulation, described as a MINLP model. Following this, we consider only the segregated flow component of this model, which can often be solved as a linear program (LP). In addition, sufficient conditions are derived where the segregated flow case yields the optimal reactor network. These conditions are based on concepts of attainable regions. Moreover, when the sufficient conditions do not hold, these concepts lead to a constructive network synthesis algorithm, which is described in section 4. This procedure is illustrated in section 5 with numerous examples. Finally, section 6 offers directions for future work. 2. Development of a Targeting Model In a broad sense, reactor synthesis can be accomplished through superstructure- and target-based synthesis routes. As described earlier, in a target-based approach, we ask the following question: What is the maximum value of the performance index irrespective of reactor type and configuration? Achenie and Biegler (1988)have shown that a two-compartment mixing model of segregation and maximum mixedness provides a sufficiently rich representation for predicting targets for isothermal systems; it has also been extended to nonisothermal systems. They modified a two-compartment mixing model of segregation and maximum mixedness (Ng and Rippin, 1965) to formulate an optimal control problem for deriving the target. Here the original Ng and Rippin model was extended so that the transfer coefficient was an arbitrary function of age. However, it was noted that though they were able to provide targets for many complex reacting systems, the formulation has some restrictions regarding the variety of mixing states. (However, it should also be noted that other mixing models can also be treated with this approach.) In this study, we propose two, more extensive, targeting models, the multicompartment approach and the constructive approach, which encompass a wider variety of mixing states and reactor network alternatives. We defer discussion of the constructive approach until the next section, as it stems from the former approach as well as the attainable region concepts of Glasser et al. (1987). In the multicompartment model approach the model essentially consists of one segregated zone and multiple zones
302 Ind. Eng. Chem. Res., Vol. 31, No. 1,1992
of maximum mixedness. A schematic of the model is shown in Figure 1. Since segregation (seg) and maximum mixedneas (max) have been recognized as the two extremes of reactor behavior, a mixing model between these two environments can simulate a wide variety of mixing states. Within the segregated zone (seg), only molecules of the same age are allowed to mix with each other, whereas in the maximum mixed environment (rnax), molecules of different ages but the same residual lives are well mixed. The residence time distribution function (RTD) denoted by f(t) is defined as the fraction of molecules that have a residence time t within the reactor network. Thus for isothermal systems within the segregated environment (and single-valued rate vectors in concentration space), the residence time distribution determines the product composition, as the concentration is just a function of the age of the molecule (Zwietering, 1959). In the maximum mixed environments, the concentration is a function of the residual life. The question we ask then is “What should the residence time distribution function and the transfer functions between these environments be so that some performance index for the entire system is optimized?” The key point is that the model does not assume any specific reactor types. In Figure 1, Xo is the concentration of the feed that is entering the reactor network. Some molecules leave the seg environment and proceed directly to the reactor exit, while others may choose to travel through the plethora of max environments and mix with molecules of different ages before they exit. Xeritis the concentration of molecules leaving the reactor exit. Transfer can occur only from the seg environment to the max environment (shown by broken arrows), because molecules that are mixed cannot spontaneously unmix. We also assume that the rate of transfer at any age of molecules from one environment to another is directly proportional to the amount of molecules in the donating environment. With this assumption, the targeting model can be formulated as an optimal control problem, where the control profiles include the residence time distribution and the mixing or transfer functions. As shown in Figure 1,the transfer coefficient from the segregated zone to the ith max zone is given by hi(a),where a is the age of the molecule, and the transfer coefficient from the ith max zone to the jth max zone is given by hi,(X),where X is the residual life. Since these are general functions of age and residual life, various degrees of mixing are simulated. Thus, optimization with respect to hi, hij, and the RTD profile leads to a general representation of the performance of the reacting system. Let X,, denote the dimensionless concentration of molecules in the seg environment, and Ximarthe concentration of molecules in the ith max environment. Also, let m, = mass contained in the seg environment as a function of age, mi = mass contained in the ith max environment, R(X) = rate vector, g(a)= Ci.fihi(s)ds, 7 = mean residence time, and J = objective function. Following the derivation in Appendix A, the isothermal formulation for the targeting problem is given by
P1:
m m J(Xexitt7) dm,/da = -Chi(a) m,(a) L
dmi/da = hi(a)m,(a)
+ EXmhji(X) mj(a)dX J
Xseg(a=O)= xo Ximax(~=-) = finite
Jmtf(t) dt = 7 logical constraints: (i) y, 5 yi-l for i = 1, ...,j =$ unit j exists only if unit j (ii)
- 1 exists
hi I U i ~ i hij 5 Uoij Yij 5 Yi, Y j
(iii) where yi and yij are binary variables: yi = 1 if maxi exists, 0 otherwise yij = 1 if there is transfer between maxi and maxj The model is a mixed-integer differential algebraic optimization problem and is more general than the twocompartment model of Achenie and Biegler. However, the difficulties associated with the discrete decisions are actually circumvented because of the iterative approach we propose in the next section. Then, for all practical purposes, one may solve this model as a differential algebraic optimization problem (DAOP). The advantage of the multicompartment model, P1, lies in its ability to cover a very wide variety of mixing states. On the other hand, P1 can be difficult to solve and the generation of the network may not be always possible from the mixing profiles. To address these questions, we first look at some simplifications to model P1. These serve as a starting point for the development of the constructive approach in section 4. 3. The Segregated-Flow Approximation A complex model, such as the multiple-compartment model, is not needed in many cases. As a matter of fact, we see that models even simpler than the Achenie and Biegler model (1988) may provide globally optimal solutions in some cases. It seems, therefore, that a stagewise approach toward the reactor synthesis problem should be a viable alternative, where one could solve simpler models, derive sufficiency conditions for optimality of these models, and then go on to increasing orders of complexity whenever these conditions are not satisfied. This is precisely the
Ind. Eng. Chem. Res., Vol. 31, No. 1,1992 303
Figure 2. Segregated flow model.
theme of the approach that we discuss in this section. The simplest model for targeting of reador networks is the segregated flow model. Here, we assume that only molecules of the same age are perfectly mixed and that molecules of different ages mix only a t the reador exit. The performance of such a model is completely determined by the residence time distribution. A schematic of this model is shown in Figure 2. This model is essentially a limiting case of the targeting model P1 as all the mixing functions hi, hi. 0. By enforcing this in P1, the isothermal formuhion for this targeting problem can be given by
-
P2
m= J(XeXit,7) dX,/dt
= R(X,,)
Here, the objective function, J , is specified by the designer and can be any function of Ltand T. One can see that the differential equation system can be uncoupled from the rest of the model and can be solved offline if Xo is a constant vector. Once the X.,'s are determined, if Gaussian quadrature is applied to the above model on finite elements on the domain [O,tmJ, we get P3:
max J(Xexild
Figure 3. Example for PFR projections.
sufficiency conditions for this model. SufficiencyConditions. We first introduce the concept of attainable regions before developing the suffciency conditions. Horn (1964) was the f i s t to suggest the idea of attainable regions for reador optimization. Glaaser et al. (1987) defined an attainable region as a region in concentration space that is attainable from some base concentration by using the processes of mixing and reaction and cannot be extended any further by these processes. The sag flow model deeeribed by problem P3 itself forms a basis to generate an attainable region. We now develop conditions for the closure of this space with respect to the operations of mixing and reaction by means of a PFR, a CSTR, or a recycle reador (RR). Let A represent the region depicted by the constraints of P3. Our aim is to develop conditions that can be checked very easily for the reaction system in question so that, if these conditions are satisfied, one needs only to solve P3 for the reactor targeting problem. We will analyze these conditions based on projected PFR trajectories in two-dimensional space. A PFR is an n-dimensional trajectory in concentration space parametzic in time which is generated by the solution of the initial value differential equation system in P2. By pairwise grouping of the eoncentrations and analyzing the projections of the actual PFR trajectory on this R2space, we develop some interesting properties. Figure 3 illustrate a PFR trajectory in three-dimensional space and ita projections. Note that the solid line represents the actual PFR trajectory and the dotted lines represent the projected trajectories. A reaction system consists of a set of components i = 1,2, n. Associated with each is a reaction rate ri. Let I = (i]denote the complete set of reacting species. This set can be further classified into Il = (a) = {a E I: r. 5 0)
...,
Here, fCij) and X (ij)correspond to the values in the ith finite element a J j t h quadrature point; i is the index of fmite elements;j is the index of the quadrature points, w, and a, are constants that denote the weights of Gauss-Legendre quadrature and the fixed length of the ith finite element, respectively. It is clear that the above formulation is linearly constrained and the solution of P3 gives us a global optimum in segregated flow for any concave objective function. However, we can reduce the above problem to a linear program for both yield and selectivity objective functions by applying suitable transformation techniques (Appendix B). The above problem can then be solved very easily given that efficient algorithms already exist for linear programs. The solution to this problem provides a good lower bound to the targeting problem. Also, we see that, in some cases, the seg flow model is sufficient for the reactor synthesis problem. The solution to this simple formulation thus could be chosen to represent the first stage in the iterative approach suggested earlier. It therefore seems appropriate at this juncture to discuss the
1 2 = f i ] = ( i E I :r j Z O ) I3 = (k]= {k E I: rh Z 0 or r, 5 0) In the most common case, our objective function is an explicit function of elements of 13,which can represent products of the reacting system. Property 1. If the projected (onto R2)PFR trajecbriea in concentration space are all such that {Xi:i E 11,131 defines a concave trajectory with respect to {X.)and no two PFR projections meet within the bounds of possible concentrations, then each of these projected spaces is closed under the operations of mixing, additional PFRs/CSTRs/RRs, or any combination of these. This property enables us to analyze the sufficiency of projeded concentration spaces. The proof is baaed on the bounding properties in Glasser et al. (1987) and is outlined in Appendix C.
304 Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992
Corollary 1. For any system that satisfies the conditions of property 1,any reactor starting from the projected space will have an outlet concentration such that, if AXi is the change in concentration of component i accomplished by the reactor and AX, is the corresponding change in the concentration of A, then
where rto and raOare reaction rates at the reactor inlet. On the other hand, if the X , vs X , projections are convex but still satisfy the condition that no two PFR projections meet, then the inequality above is reversed. Details are provided in Appendix C. At this point, the utility of this property with respect to problem P3 deserves attention. A careful look at P3 reveals that the shaded region in the projected space (for example, the XI-& space) is exactly the projection on the XI-& space of the feasible region of P3. The concave PFR projection defines the concentrations in segregated flow, and the interior is a convex combination of all boundary points created by the residence time distribution function. This gives a new interpretation to the residence time distribution as a convex combinator. For any convex objective function to be maximized, the solution to the seg flow model will always lead to a boundary point of A. The demonstration of these sufficiency results will become evident in section 6 when we discuss the example problems. Fortunately, since we have closed-form expressions for the slopes through the rate equations (though we do not have closed-form expressions for the trajectories themselves), checking for concavity is easy. However, it is clear that we require certain concavity properties among the projections, i.e., we need to know a priori whether certain projections are concave. This task can be simplified by developing certain properties for mutual concavity in the reacting system. Property 2. If the projection of any element of {XI) vs (X,) defines a conuex curve and the projection of any element of (&) v8 (X,)defines a concaue trajectory, then the are concave up to the projected curves of (xk)vs (XI) stationarypoint in this projection, and this stationary point is unique if strict concavity holds in the xk-x, projection. The proof is outlined in Appendix D. Corollary 2. If the projections of { X ) vs {X,)define concave curves and the projections of vs {X,) define concaue trajectories, then the projected curves of (X,) vs (Xi}are concaue beyond the stationary point in this projection and this stationary point is unique if strict concavity holds in the xk-x,, projection. The proof is based on the same arguments as for property 2.
(dk}
4. A Constructive, Stage-Wise Approach for Reactor Synthesis The propertiea described above enable us to identify the nature of different projections and can be very useful, as shown in the example problems. It is essential to note, however, that property 1is only a suffcient condition, but not always necessary. In other words, the seg flow model can be optimal even if these conditions are not satisfied. If the LP region (for P3) is not sufficient, we would like to generate a comtructiue procedure for extending this region, by superimposing various reactors on the region provided by the linearly constrained formulation. The main idea for the constructive approach is Given a candidate for the attainable region, can extensions to this region be generated? If yes, then create
the extension and, on its conuex hull, check for further extensions that improve the objectiue function. This procedure is continued until there are no extensions that improue the objective function. The development of an efficient algorithm for this constructive approach (based on formulation and solution of small nonlinear programs) is the main focus of this section. As described before, G h r et al. (1987) have developed a set of neceeaLLly conditions baaed on geometric arguments for an attainable region. However, their approach, though insightful and elegant, is limited to three dimensions, and the implementation is based on physically drawing alternate PFR and CSTR trajectories, which can be very cumbersome. We seek procedures that can be implemented and manipulated more easily and take advantage of some of the properties of the seg flow model described by problem P3. The main insight in the proposed algorithm is that the RTDs are convex combinators and the region enclosed by the segregated flow model is always convex. The aim now is to develop an algorithm by which, given a candidate for an attainable region, we should be able to check whether it can be extended to our advantage. Here, we restrict ourselves to PFR, recycle reactor, and CSTR extensions only. Consider the feasible region for P3, 34, as the first candidate for the attainable region, which is a convex combination of the concentrations in segregated flow. Each combination of the RTDs and the concentrations gives a unique point in the feasible region. The following subproblem can also be embedded within P3, in order to check whether, starting from any point in A, a reactor can provide an extension that gives a point outside of A. Here, we illustrate a recycle reactor extension, since it includes the PFR and CSTR extensions as special cases. Simpler formulations for purely CSTR or PFR extensions can be developed along the same lines. Recycle Reactor (RR) Extension. If J , > Jp3,then the recycle reactor provides an advantageous extension outside of A. P4 m a Jm(xexit)
X,,(t=O) =
RSexit + X P , Re + 1
(3)
In addition, if any of the projections Xi vs Xi Q’E 1 2 ) are concave and if we wish to eliminate searching in the interior of these projections, we can always include the following inequality, which arises out of geometric arguments (corollary 1):
Here, J , is the value of the objective function at the exit of the recycle reactor. Equation 1describes the concen-
Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992 305
t
,+
xmodei
LP region
-
Figure 5. Illustration for extansion of convex hull.
Form nsw c m x hul of mnqnmtim
CCf(ii) + Ffmodel(k) i l
1.0
In this formulation, Xm&l(k) is a constant vector and refleda the concentration at the exit in the models previously
chosen. A convex combination of this with the seg flow region described by P3 gives the fresh feed point for the recycle reactor we are looking for, Xupdate.Xeit then represents the concentration at the exit of the recycle reactor; and if AXexi,) > J(Xmdel),the earlier model chosen is sufficient. The control variables essentially are the f s and fm&l(k), which are the linear combinators used to provide a convex candidate. A careful look at the formulation rev& that we are checking for completeness of the convex hull of the region found by the model. Figure 5 gives a geometrical interpretation to the solution of P5. If the solution of P5 indicates that the objective function can be improved by extending the attainable region, we consider a more complex model. Here, we therefore consider our two alternative schemes: the constructive upproach (using P3 and P5) and the multicompartment approach (using incrementally more complex versions of Pl). An obvious, but slight, increase in complexity comes from the recursive solution of P5. In the constructive approach, the expression for Xu,, automatically includes all the points in the convex hulls attained so far by previous recycle reactor extensions. Thus,this region is made up of the segregated flow model as well as favorable recycle reactor extensions from previous solutions of P5. We continue to check for extensions and terminate when there are no further extensions that improve the objective function. Note that with this constructive approach, the algorithm allows the reactor network to be synthesized readily. It is also clear that this approach shares many characteristics with the geometric approach of Glasser et al. A n important difference is that their approach searches for all possible extensions of candidate attainable regions and requires checking of an infinite number of points on the convex hull. Because this can only be done efficiently on a geometric basis, their approach has been applied only to two or three dimensions. Our constructive approach, on the other hand, automatically finds only those extensions that improve the objective function. Since it is an optimization-based approach, it is not limited by problem dimensionality or the addition of constraints, as will be demonstrated by the examples. We also note that because problems P3 and P5 are relatively small and simple optimization problems, constructive solutions for the attainable region can be found very quickly. One potential disadvantage to the constructive approach, however, is that a particular extension which does not improve the objective function may still enlarge the attainable space enough so that a point from within this extension can improve the objective function beyond what we started with. This nonmonotonic increase in the objective is a limitation of the proposed approach. Also, we note that even though the attainable space of concentrations is always convex, P5 is not necessarily a convex program and therefore will not necessarily find the global optimum that we are seeking. (However, termination at local optima is not serious because the solution of P3
306 Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992
usually yields an excellent starting point and multiple starting points can be tried, if necessary, to ensure a global optimum has been found for P5.) As our second or multicompartment approach, we consider the multiple-compartment model given in section 2. Because of the stagewise nature of the algorithm in Figure 4, it is clear that the level of complexity of the model can be increased simply by adding another mixing compartment, if it is needed. Rather than solving an MINLP for the complete version of P1, and risking problems with nonunique, or local optima, we see that the stagewise solution to P1 guarantees steady improvement in the objective and uses only the simplest model that it needs. For the multiple-compartment model, we briefly consider solution techniques for problem P1 in the next section.
5. Solution Techniques Many techniques can be used to solve problem P1. Common adjoint techniques that have been used efficiently (Achenie and Biegler, 1988) are not very attractive in this case, because the derivation of the adjoint system can be cumbersome. Furthermore, constraints can be difficult to handle, and the method requires repeated solution of differential algebraic equations consisting of the model and the adjoint equations. We wish to take advantage of the nature of the problem and therefore choose a simultaneous solution and optimization strategy. By applying orthogonal collocation on finite elements, we can reduce the problem to a nonlinear program (Cuthrell and Biegler, 1987). The integrals within the optimization problem make this approach very attractive, because they are now automatically evaluated at Gaussian quadrature points. To solve the above problem, the control profiles (f and h ) were assumed to be piecewise constant on each element, and the state profiles (i.e.,),X were approximated by Lagrange basis polynomials in each element. Substituting these approximations into the general DAOP (Pl), we can reduce P1 to a nonlinear program (NLP) in terms of the values of the state and control variables at the collocation points in each element (the shifted roots of an orthogonal Legendre polynomial). As with the targeting problems in the constructive a p proach, the resulting nonlinear program was solved by using MINOS 5.2 (Murtagh and Saunders, 1978), which is interfaced to the GAMS modeling system (Brooke et al., 1988). Thus, in this formulation the modeling equations need not be solved completely except at the optimum. Also, as all variables and equations are part of the optimization problem, state variable constraints can be imposed directly, and there is no need to derive the adjoint equations. However, the price one pays for the advantages noted above is that we are now solving a larger problem and need to develop a good initialization method. Fortunately, the solution to the linearly constrained segregated flow approximation (P3) gives an initial guess for the residence time distribution function and is a very good initializer to the reformulated nonlinear program. The solution to this reformulated NLP of P1 then gives us the RTD and mixing profiles corresponding to the optimum performance index. 6. Example Problems In this section we present seven example problems to illustrate our approach, and give a brief discussion and comparison with literature results. The first three examples consist of nonlinear reaction systems for which the sufficient conditions for optimality of the seg flow model (linear programs in these cases) hold. Ths remaining examples involve nonlinear reaction systems where these
1 .o
0.288
0
t (secs)
Figure 6. RTD for example 1.
conditions are not satisfied; the performance of the stagewise approach is illustrated by these examples. Example 1. The isothermal Van de Vusse reaction involves four species for which the objective is the maximization of the yield of intermediate species B, given a feed of pure A. The reaction network is given by ki
A-0-C
kz
lk3
n
u
The reaction rate vector for components A, B, C, and D, respectively, is given by R(C) = [-k,cA - k & ~ ~klCA , - k& k & ~ , k&~'] where kl = 10 s-l, k, = 1 s-l, k, = 0.5 L mol-' s-l, CAO = 0.58 mol/L Dedimensionalizing with respect to CAo,we get
R(X)= [-1oxA - 0.29XA2,
lox, - XB,
XB, o.29xA2]
where XA = CA/CAO, XB = CB/CAO The objective function is the yield of component B. Following the algorithm in Figure 4, we find that the solution of P3, which is an LP, gives a globally optimal reactor network. This is because the reaction kinetics satisfy the conditions in property 1(i.e., the XB vs X A curves are all concave). The verification of concavity is simple, because we just have to verify the monotonicity of the slopes, given in closed form for the reactor synthesis problem. Here, the slope of the PFR trajectories in XB w XAspace is just given by rate(B)/rate(A), and a simple program to find the maximum of the second derivative can be solved. The optimal value of the objective function is given by XBexit 0.7636. This is the globally optimal solution and can be realized by a PFR with a residence time of 0.288 s. The residence time distribution predicted by the LP formulation is shown in Figure 6. A brief comparison with results from the literature is presented below. author Chitra and Govind (1981) Achenie and Biegler (1986) Achenie and Biegler (1988) Kokoseis and Floudaa (1990) P3 formulation (LP)
yield 0.752 0.7531 0.757 0.752 0.7636
residence time 0.2551 0.2965 0.2370 0.2539 0.2880
The results from the different approaches are nearly the same. However, the slightly better objective function predicted by the LP formulation may be attributed to the better approximation procedure when the differential equations are solved offline. Example 2. The a-pinene problem is a reaction network that consists of five species and has the following reaction
Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992 307
hex
Figure 7. RTD for example 2.
network. The objective function here is the maximization of the selectivity of C over D given a feed of pure A. A-E
C
s
1"
D
k7
ki
(7) ks
s
B
k3
The reaction vector for the components A, B, C, D, and E, respectively, is given by R(X)= [-(k1 + k2)X~- 2 k & ~ ~ , -~ k,Xc, -k6X~ k & ~ , k&A2 + k 4 X ~ k2XA k6Xb - k3XD - 2k4XD2+ 2k7&, klXA] where Xj = CJCAO and CAO = 1 mol/L
+ +
kl = 0.333 84 s-l, k2 = 0.266 87 s-l, k3 =0.14940 s-l, k4 = 0.18957 s-l, k5 = 0.009598 s-l, k6 = 0.294 25 S-l, k7 = 0.011 932 5-l This problem was solved by using P3 and the RTD predicted is shown in Figure 7. Interestingly, the RTD settles at the upper bound in the time horizon. Thus,we find that the selectivity in this problem increases monotonically with age within the segregated flow model. For instance, the selectivity obtained with a t,, of 60 s was 1.48. This example corroborates some of the strengths of a targetbased perspective that we mentioned earlier. Achenie and Biegler (1988) arrived at an optimal selectivity of 0.2336 by placing a bound on the residence time of 6 s. Kokossis and Floudas (1990), using a superstructure-based approach, arrived at a complex network of a CSTR and a PFR (represented by sub-CSTRs) with recycles and intermediate feeds to get a maximum selectivity of 1.4020. The superstructure-based approaches, therefore, may suffer from the problem of nonuniqueness, in that there is no difference between many configurations. Moreover, P3 (LP in this case) is based on a finite approximation to an infinite residence time domain. Thus, from the solution it is easy to detect monotonicities, etc. (monotonicities are clearly a possibility when the solution to the seg flow model occurs at the maximum age). Example 3. The isothermal Denbigh reaction involves five species and has the following reaction diagram: kz
ki
A-B-C /ks
/h
D
E
We now examine two cases. (a) Consider the yield of intermediate B as the objective function. In this case, again the conditions in property 1 are satisfied; i.e., the XB-XA curves are all concave and thus and the LP formulation gives the global optimum. The optimal X Bis 0.2578 ~ ~ in~ a PFR ~ with a mean residence time of 0.39 s. (b) Consider XB/XDas the objective function. Here, there are intersections in the XB-XD space. However, since the XBv8 XAtrajectories satisfy the conditions in property 1, we have for any recycle reactor line (from corollary 1)
The dedimensionalized reaction rate vector .for components A, B, C, D, and E, respectively, is given by -x~(k4X + ~k2), R(X) = [-XA(k3 + klXA), 0.5klX~~ k2XB9 k 3 x A ~ k4XB21 where kl = 6.0 s-', kz = k3 = 0.6 s-', k4 = 0.6 s-l, CAO = 6 mol/L, CDO = 0.6 mol/L
AXA I O Again,the XD-XA curves are convex, and we get (corollary 1)
It can be shown that it is sufficient to consider only those points such that rB0 > 0. Therefore, from (7) by multiplying and dividing by AXD
>O
Using (81, we have for rgo 1 0, then:
(2)(
o):
However, from property 2 we must have that the XBvs XDc w e s are concave as long as rB 1 0. From concavity, we infer that inequality (9) can lead to an extension that is only inside the LP region (converse of corollary 1). Hence, the selectivity obtained with the LP formulation is the globally optimal one. This optimal selectivity is 1.322 in a PFR at a mean residence time of 0.209 s. Both Achenie and Biegler (1988) and Kokossis and Floudae (1990) obtained similar results, though they did not address the question of global optimality. Example 4. Here, we consider the Denbigh problem again, with an altered objective function and rate vector as follows: R(X) [-x~(k3+ kiXA), kiXA2- X B ( ~ ~ +Xkz), B k2XB9 ~ X A Ik4XB21 where kl = 6.0 s - ~ , k2 = k3 = 0.6 s-~,k4 = 0.6 s-', CAO = 6 mol/L The objective here is to maximize the production of C subject to 95% conversion of A. The conditions for P3 optimality are not satisfied here. Following the stagewise approach, we observe a CSTR extension from the seg flow model, which gives a C yield of 3.726 mol/L with a CSTR residence time of 3505 s. The residence time in the segregated environment is 0.766 s with a Dirac S function for the RTD, which corresponds to a PFR. No further extensions are observed when solving P4 again. By mixing with the feed 95% conversion of A can be achieved, which corresponds to a C yield of 3.54 mol/L. Glasser et al.
308 Ind. Eng. Chem. Res., Vol. 31, No. 1,1992 . rl.
Fa
-tirod
I
0
Figure 8. RTD and mixing functions for example 6.
(1990) observe similar results with a CSTR with infinite space time. Example 5. The Trambouze reaction involves four components and has the following reaction scheme: k,
A-B
with k l = 0.025 mol/L,
kz
A-C
k3
A-D
k2 = 0.2 min-', k3 = 0.4 L/(mol min)
We wish to maximize the selectivity of C to A defined by Xc/(l - XA). Again, the conditions of P3 optimality are not satisfied, and following the stagewise approach, we arrive at a CSTR extension from the feed point of the segregated flow model, indicating a single CSTR with a selectivity of 0.5 and residence time of 7.5 s. No further extensions are observed by solving P4. Achenie and Biegler (1990) observe a selectivity of 0.4999 in a two-CSTR combination. Kokossis and Floudas (1990) report many solutions to this problem; Glasser et d. (1987) observe that this problem has an infinite number of optimal solutions with a selectivity of 0.5. Example 6. Here, we revisit the Van de Vusse reaction with altered rate constants. The objective function again is the yield of intermediate species B. The rate vector is given by R(X)= [-XA - 20xA2, XA- 2xB, 2xB, 2oxA2] This was solved by using both the multiple-compartment approach and the constructive approach. In this case, the seg flow model gives a yield of 0.061. However, the sufficiency conditions for the LP formulation are not satisfied. Following the algorithm (Figure 41, we find that a twocompartment model is sufficient. The optimal value of the objective is given by Beit = 0.0705 with a residence time of 0.25 s. The control profiles, namely, h (transfer function from seg to max environment) and f, are as shown in Figure 8. Glasser et al. report a yield of 0.071 with a geometric approach. Using the constructive approach, by using recycle reactor extensions, we observe a B yield of 0.069. The network corresponds to a recycle reactor from the feed point (recycle fraction = 0.772; residence time in the plug section of RR = 0.1005 s) in series with a PFR with a residence time of 0.09 s. The lower yield obtained by the constructive approach in this example can be attributed to a nonmonotonic increase in the objective, as we discussed earlier in the paper. Example 7 (Integration of the Reactor with the Flowsheet). It is well-known that the reaction chemistry determines the character of the entire process (Douglas, 1989). The recycle structure and downstream processing steps for the flowsheet are directly influenced by the conversion of raw materials to selective product. Most of the studies on reactor network synthesis, however, have considered the reactor in isolation. The targeting approach, coupled with the simultaneous solution strategy presented before, allows for the integration of the reactor with the flowsheet. Though this integrated approach is independent of a particular reactor network, it is effective
I
Figure 9. Williams and Otto flowsheet.
because the capital cost of the reactor is very low compared to raw material and downstream processing costs. Here, we replace the reactor within the flowsheet by a suitable targeting model. The initial and final conditions for this model may be related with the variables of the flowsheet such as the feed rate and recycle ratio. Using the stagewiae approach, we can find the best reactor network for all the possible initial concentrations dictated by the constrainta posed by the flowsheet variables. Our results indicate significant savings when one considers the integration of the reactor with the flowsheet. To illustrate this concept, we consider the Williams and Otto flowsheet problem, which has been often described as a typical flowsheet optimization problem (see Vasantharajan and Biegler (1988) for a review of this problem). The schematic of the flowsheet is shown in Figure 9. The plant consists of a reactor, a heat exchanger to cool the reactor effluent, a decanter to separate a waste product G, and a distillation column to separate product P. A portion of the bottom product is recycled to the reactor, and the rest is used as fuel. The following reactions are involved in the manufacture of compound P: A+B-C (reaction 1) C+B-P+E (reaction 2) P+C-G (reaction 3) The rate vector for components A, B, C, P, E, and G, respectively, is given by R(X)= [-kiXAXB, (hxA+ k2Xc)X~, 2kiX~x 2kzXBXc ~ - k3XpX~, k&BXc - 0.5kJpXc, Bkpx~Xc, 1.5k,XpXc] where
kl = 6.1074 h-l wt fraction, k2 = 15.0034 h-' wt fraction, k3 = 9.9851 h-' wt fraction The X s here denote the weight fractions of the components. In Figure 9 Fa and Fbare the flow rates of fresh A and B, respectively; Fgis the flow rate of waste G; and Fpis the fixed exit flow rate of pure P out of the plant. Researchers have solved this problem by assuming the reactor to be a CSTR and optimizing the volume of this CSTR to maximize the rate of return on investment. We replaced the CSTR by a one-compartment segregated flow targeting model embedded within the flowsheet. The objective function, namely, the return on investment (ROI) of 130% has been typically obtained for this problem with the fixed CSTR model. With a one-compartment model integrated within the flowsheet, an ROI of 278% was obtained. We now look for CSTR extensions from the onecompartment model by solving P4 for a CSTR extension by including all the constrainta imposed by the flowsheet. No CSTR extensions that improve the ROI are observed. It must be noted here that in this case, since the space of initial conditions for the seg flow model is infinite but bounded by flowsheet constraints, the differential equations in P2 cannot be uncoupled and have to be solved
Ind. Eng. Chem. Res., Vol. 31,No. 1, 1992 309 Science Division of Argonne National Laboratory for facilities and support.
x. Figure 10. PFR projection.
simultaneously. The optimal network is just a PFR with a residence time of 0.0111 h. These results indicate that significant savings can be obtained by integrating the reactor with the flowsbeet, even with very simple targeting models. 7. Conclusions
A constructive targeting procedure for reador network synthesis has been developed. Here, we highlight a few basic differences between superstructure- and target-based approaches to reactor synthesis. Superstructure approaches have the advantage of generating the network along with the solution of the optimization problem; capital costa can also be incorporated easily. However, as indicated in example C, the nonuniqueness in reador synthesis can sometimes give target-based approaches an upper band. The main idea in our targeting approach is based on successively generating points of the feasible region for the reaction system with optimization over the convex hull of these points. The example problems indicate that, in some cases, the segregated flow model is sufficient to describe the network. When the seg flow model is not sufficient, simple nonlinear programs can be solved to enhance the target. Thus a stagewise approach, where one could solve simpler models and verify their sufficiency, becomes more attractive than the complete solution of P1. We describe two approaches for increasing the complexity of the targeting model. With the constructiue approach, convex regions are extended recursively and synthesis of the network is straightforward. With the multicompartment approach, additional compartments are added at each stage. These targeting methods are illustrated on several literature examples with very good results. With the exception of example 6, the results, in general, favor the constructive approach over the multiple-compartment approach, the advantage mainly b e i i the ease of solution and realization of the network with the former approach. The challenge still remains in extending the above techniques to nonisothermal reactors and reaction-separation systems. We are currently considering reactionseparation systems where, again, one can define a generalized residence time distribution function and mixing functions. However, unlike the reaction kinetics that constrain the attainable space, there are fewer constraints with separation. Here, often we will be able to attain the entire stoichiometric space in concentration. The choice of a suitable economic objective function thus plays a much more imporbnt role in the success of this approach. These issues will be addressed in a future publication. Acknowledgment This work was supported by the National Science Foundation and the Computer Aided Process Design Laboratory, an industrial consortium of the Engineering Design Research Center at Carnegie Mellon University. Thanks are also due to the Mathematics and Computer
Nomenclature a’ = recycle fraction C = molar concentration vector C,, = reference concentration (usually inlet concentration of reactant) g = integral of cumulative transfer function from seg zone = Z;lih;(s)ds f = residence time distribution function f. = linear combinator of all concentrations from PFR section of the recycle reactor fmh,ca, = weight for concentration obtained from solution of model at kth iteration h, = transfer function from the seg zone to the ith max mixed zone h;j = transfer function from the ith max mixed zone to the jth max mixed zone J = objective function for reactor synthesis J, = objective function attained hy a recycle reactor extension Jp3= objective value after solution of P3 1 = lower bound on XeXit n.= mass contained in the seg environment (function of age) mi = mass contained in the ith mas environment rj, rj = reaction rates for components i, j, etc. ra, rN = reaction rates for components i, j, etc., at reactor inlet conditions R(X)= rate vector (dedimensionalized to 8-9 Re = recycle ratio t = residence time u = upper bound on X . U,= upper bound h; U,= upper bound on hjj V = reactor volume w = weights for Gaussian quadrature X = dimensionless concentration vector (DCV) = (C/Cmf) X j = dimensionless concentration for component j X , = Concentrations available fmm solution of P3 (parame!& in f , AXj = change in X j due to reaction Xo = dimensionless concentration vector at reactor entry Xeeg= DCV within the segregated environment Xi, = DCV within the max mixed zone i XeXit= DCV at reactor exit Xu DCV for fresh feed to recycle reactor extension in $5
-
Xf= fresh feed concentration (dimensionless) for recycle reactor
X, = concentration profiie in plug section of recycle reactor y; = integer variable denoting existence/nonexistenceof max
unit i yij = integer variable denoting existencelnonexistenee of
transfer between rnax unit i and max unit j
Greek Letters a = age of a molecule X = residual life for molecule T = mean residence time Appendix A. Derivation of Targeting Model Pl Based on the nomenclature described in the paper, we can derive the targeting model P1 as follows: Within the seg environment, since hi(.) is the rate of transfer of material at age 01 to the different zones, by a mass balance, we have
dm./da = -Ehi(a) m.(a) where mJ0) = 1.
(-4-1)
310 Ind. Eng. Chem. Res., Vol. 31, No. 1,1992
Similarly, a mass balance on the ith rnax zone gives dmi/da =
The exit concentration is then given by the following weighted average:
hi(a)m,(a) + CJmhjimj(a)dX - ZJmhijmi(a)dX i
I
o
(A-2) These equations give us expression for the bulk flow of material in each region of the model. The concentrations can now be calculated a function off and mixing profiles by writing mass balances over the mixing environments. Within the seg environment, since the mole fractions are described by the age of the molecule, we have the following species balance: dXseg/da = R(Xseg)
It can be easily shown that, from (A-1)
Also, for the residence time distribution we must have Jmf(t) dt = 1
Xseg(0) = XO Within each rnax environment, the species balance will be performed over a differential residual life. The volume of material in the segregated zone with residual life between X and X + dX is given by
The logical constrainta follow in a straightforward manner: (i) y j 5 yi for i = 1, ...,j-1 * unit j exists only if unit j - 1 exists
-
(ii) hi I U,yi In the rnax zone it is given by V(dh/r)Smmi(a) 0 f(a+X) da
transfer to ith max occurs only if it exists hij I Ugij
yij Iyi, yj; VI, U, are upper bounds on hi and hij
The moles of the components within the seg and rnax environments within X and X + dA are then given by V(dh/T)~mX,,g(a)m,(a) f(a+X)da
(iii)
C y j 1 1 * at least one rnax unit exists
where yi and yij are binary variables: yi = 1 if maxi exists, 0 otherwise yij = 1 if there is transfer between maxi and maxj
and
The mass balance in the rnax zone i can now be written as (accumulation of species in MMJ = (rate of supply from seg) (rate of supply from all MMjj+i) - (rate of loss of MMjjzi) (rate of consumption by reaction)
+
CCwjf(ij)a(i) = 1
+
Note that mixing occurs only at the same value of A because of the assumption of maximum mixedness. Thus, we get -~(X1m,(~)Jmmi(a) d . f(a+h) da) = JmXseg(a) hi(a) m,(a) f(a+A) da CXim,(X) hij(h)Smmi(a) 0 f(a+X) d a I
Appendix B. Transformation Techniques for Linearization of the Seg Flow Model P3 The formulation for seg flow as given by P3 is max J(Xexit,T)
+
This equation can be solved by using the Zwietering boundary condition (Zwietering, 1959),i.e. dXi,, (A=..) =0 dX However, Glasser et al. (1986) have shown that, for many reaction systems, different values of Xim&==) yield the same values of Xim,(0) needed in the expression for the average exit concentration. Thus, any finite value for Xim,(X=m) will solve the above system.
i l
7
CCWjf(ij)t(ij) a(i) i
j
Xexit = CCwjf(ij)&eg(iJ) d i ) i l
(B-1)
(B-2) 03-31
We wish to transform the above model to a linear program when the objective function involves expressions such as selectivity or a general form such as J(Xexit,r) = aT[XexitT,TI T/bT[XexitT,~] where a and b are constant vectors and the superscript T is the transpose. For convenience in bookkeeping, we can replace w j f ( i j ) a(i) by just a new f ( i , j ) . Also, the double subscripts can be removed without loas of generality to get the equivalent problem: max J(Xexit,T) Cf(i)= 1 (B-1’) i
Cf(i)t(i)
(B-2’)
= Cf(i)Xseg(i)
(B-3’)
r =
i
&it.
I
Let us define the following scalars with the assumption that bT[XeetT,7lT# 0:
Ind. Eng. Chem. Res., Vol. 31, No. 1,1992 311 z = l/bTIXexjF,~]T,yT= 72
and the following vectors: = Xesitz, Yf = fz where f is the row vedor of the f(i)’s,and where yfi denotes the elements of yr. Multiplying (B-1’), (B-2’), and (B-3’) throughout by z gives ~x
Zf(i)z= z * Cyti
=z
E
1
72
because any extension beyond A from the interior means an extension from the boundary, since the mixed feed point, if one exists, must always lie outside A. From (C-21, it is clear that points A, B,and C are collinear. It follows from the definition of a concave function $(x) (which in this case is the function for V,, $(X,)) that, if X,, > X,C, then for any X , > X,, we must have
= C f ( i ) z t ( i )* Y T = CYfit(i) 1
1
Also z = l/bTIX,sitT,~]T * bT[XeXitTz,7zlT = 1 bT[yXT,yTlT =1
-
Finally J = aT[yxT,yT1
Thus the reformulated model is max J = aT[yST,yTIT CYfi = z
(B-1”)
1
Y, = CYfit(i)
(B-2”)
i
This means that, for all points X , > X&, $ lies below the line connecting (X,c,$(X,c)) and (X,B,$(X,B)). Since B and C are points on a concave trajectory, we must have from (C-3) that the PFR trajectory is bounded above by line AB. Now, if a reverse PFR (a PFR with the rate vedor negated) is employed from point B,then since X , has to increase monotonically, V2 must meet V1 because V2 always lies below AB and A E V1. This again contradicts our condition that the PFR trajectories cannot intersect. Since A is any arbitrary point on the PFR trajectory and only those extensions are possible such that X , decreases, this contradiction must hold for every point on A. Hence, our assumption that a recycle reactor could extend the attainable space is false, and the projected space is closed under recycle reactor operation. Property l a (CSTR Extensions). It is clear that in the case of the CSTR, both points B and C coincide, and the slope of line BC becomes the slope of the PFR trajectory at point B. Here, ((3-3) becomes
(B-3”) bT[yXT,yTlT =1
(B-4”)
This is the linear programming representation of the original problem. Appendix C. Proof of Property 1 Consider any such projected space of the trajectories (Figure 10). Let X,,u E I,, be the abcissa and Xkbe the ordinate. The proof is essentially geometrical and is based on the bounding properties in Glasser et al. (1987). The attainable space in this projection is shown by the shaded region, and let us call it A. Given that all the PFR trajectories are concave and do not meet each other, we have to prove that no extension of this space is possible. If we reatrict ourselves to PFR, CSTR, and recycle reactors, then it is enough to show that there is no recycle reactor extension, because the other two reactors can be shown to be special cases of the recycle reactor. The recycle reactor is defined by dX/dt = R(X) (C-1) X(t=O) = aXf + (1- u?X,
(C-2)
where Xf = fresh feed concentration, A; X, = exit concentration, C; and a’ = fraction of recycle E [0,1]. Consider any point A of the projected space (i.e., choose any point for the fresh feed). Assume that there is a recycle reactor extension from this point such that x, E complement (A) and is at the end of the PFR curve, V2. Pointa B and C denote the mixed feed point and the exit point of the recycle reactor, respectively. Since C E complement (A), the mixed feed point should also belong to complement (A) because otherwise the PFR curves V1 and V2 would meet, a situation that contradicts the condition that the projected PFR curves do not intersect. Furthermore, it is enough to analyze the boundary of A,
and
and again, V2 is bounded above by AB, because of which there cannot be any CSTR extension that goes out of A based on the arguments for proof of property 1. Property l b (Implications in Two Dimensions). In two dimensions, the projections are the true PFR trajectories themselves. The PFR trajectories can never meet because of the uniqueness of rates at a point; and thus if all the PFR trajectories are concave, we have found the complete attainable region. Proof of Corollary 1. Consider Figure 10 and let us assume that we start a reactor from any point on V1. From property 1, since V1 is a concave trajectory, the reactor cannot lead us to any point outside the shaded space; that is, the endpoint of the reactor must lie within the projected seg flow feasible region (shaded space). Also, any reactor can only decrease A. Therefore, only the left half of the shaded space is accessible from any point on VI. In order that this condition be satisfied, the line joining the initial point to the exit point must have a slope that is at least greater than the slope of the tangent to V1at that point. For example, if we start from A, the slope of the reactor inlet-outlet line w i l l be greater than that of line AB. However, the slope of the tangent is the ratio of the rate vectors. Therefore, whenever property 1holds, any reactor must satisfy the following condition: a
i
-1AXa
rio F ~ O
On similar grounds, it can be shown that for a convex
312 Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992
projection without intersections, the inequality above is reversed.
tionary point in the unique.
Appendix D. Proof of Property 2 If we choose projections Xj vs X, and Xk vs X, to be the given convex and concave projections, respectively, we need to prove the result for the Xk vs Xj projections. Let (Xao,Xio,XkO}be any reference point on the actual single PFR trajectory. Then, from the conditions imposed in property 2, we have that if {X,, xj,xk}is another point on this PFR trajectory X, 1 Xa0 implies Xj IXjo (D-1)
Literature Cited
From the convexity of this projection function ( X j vs X,), we must have increasing slopes with increasing X,. It must be noted that we are analyzing only a single PFR trajectory which is completely defined by the space time. Therefore, for X, > X,,
(g)o
dXj 0 I- I dXa
Again, by the concavity of the x k vs X, curve for X, XaO
>
(D-3) Also, dX,/dX, = rm/rnfor any components m and n. No partial derivatives are involved here, simply because all the X s are implicit functions of one variable, i.e., space time for the fixed PFR. Therefore, using implicit differentiation on (D-3), we have
0
from (D-2)
(;)( );
0
For (rk/rj)o I0, substituting the above result in (D-4) and for all Xi IXj0 we have
that is, the Xk vs Xj projections are concave up to the stationary point in this curve, or, in other words, concavity is guaranteed until the point (rk/r')o= 0 is reached. h, since dxk/dXj = ( d x k / d x a ) ( d x , , ] d x , ) , the left-hand side can vanish at most at only one point because dX,/dX, vanishes at most at only one point to satisfy the strict concavity condition, and dXa/dX.does not vanish from physical arguments in the generd case. Hence, the sta-
xk-x, projection of property 2 is
Achenie, L. E. K.; Biegler, L. T. Algorithmic Synthesis of Chemical Reactor Networks Using Mathematical Programming. Znd. Eng. Chem. Fundam. 1986,25,621. Achenie, L. E. K.; Biegler, L. T. Developing Targets for the Performance Index of a Chemical Reactor Network. Znd. Eng. Chem. Res. 1988,27,1811. Achenie, L. E.K.; Biegler, L. T. A Superstructure Based Approach to Chemical Reactor Network Synthesis. Comput. Chem. Eng. 1990,14, 1, 23. Brooke, A.; Kendrick, D.; Meeraus, A. GAMS: A User's Guide; Scientific Press: Redwood City, CA, 1988. Chitra, S. P.; Govind, R. Yield Optimization for Complex Reactor Systems. Chem. Eng. Sci. 1981,36,1219. Chitra, S. P.; Govind, R. Synthesis of Optimal Reactor Structures for Homogenous Reactions. AZChE J. 1986,31(2),177. Cuthrell, J. E.;Biegler, L. T. On the Optimization of Differential Algebraic Process Systems. AIChE J. 1987,33,282. Douglas,J. M. A Hierarchical Decision Procedure for the Synthesis of Complex Planta; Foundations of Computer Aided Process Design; CACHE Corp.: Snowmass, CO, 1989. Glaseer, D.; Crowe, C.; Jackson, R. Zwietering's Maximum Mixed Reactor Model and the Existence of Multiple Steady States. Chem. Eng. Commun. 1986,40,41. Glasser, D.; Crowe, C.; Hildebrandt, D. A Geometric Approach to Steady Flow Reactors: The Attainable Region and Optimization in Concentration Space. Znd. Eng. Chem. Res. 1987,s (9),1803. Hildebrandt, D.; Glasser, D. The Attainable Region and Optimal Reactor Structures. Proceedings of the ZSCRE Meeting, Toronto, 1990. Hildebrandt, D.; Glasser, D.; Crowe, C. The Geometry of the Attainable Region Generated by Reaction and Mixing: with and without constraints. Znd. Eng. Chem. Res. 1990,29(l), 49. Horn,F. Attainable Regions in Chemical Reaction Technique. The Third European Symposium on Chemical Reaction Engineering; Pergamon: London, 1964. Horn, F. J. M.; Tsai, M. J. The Use of Adjoint Variables in the Development of Improvement Criteria for Chemical Reactors. J. Opt. Theory Appl. 1967,I (2), 131. Jackson, R. Optimization of Chemical h c t o r s with Respect to Flow Configuration. J. Opt. Theory Appl. 1968,2 (4),240. Jackson, R.; Glaaser, D. A General Mixing Model for Steady Flow Chemical Reactors. Chem. Eng. Commun. 1986,42,17-35. Kokossis, A. C.; Floudas, C. A. Synthesis of Isothermal ReactorSeparator-Recycle System. A n n u l AZChE Meeting, San Francisco, CA; AIChE: New York, 1989. Kokossis, A. C.; Floudas, C. A. Optimization of Complex Reactor Networks-I. Isothermal Operation. Chem. Eng. Sci. 1990,45 (31,595. Levenspiel, 0. Chemical Reaction Engineering; Wiley: New York, 1972. Mu.rtagh, B.;Saunders, M. MINOSIAugmented Supplementary User's Manual; Report OR/78/6; University of New South Wales: Sydney, Australia, 1978. Ng, D. Y. C.; Rippin, D. W. T. The Effect of Incomplete Mixing on Conversion in HomogeneousReactions. Proceedings of the Third Symposium in Chemical Reaction Engineering; 1965. Vasantharajan, S.; Biegler, L. T. Large Scale Decomposition Strategies for Successive Quadratic Programming. Comput. Chem. Eng. 1988,12(10, 1087. Zwietering, N. The degree of mixing in continuous flow systems. Chem. Eng. Sci. 1969,11,1. Receiued for review November 26, 1990 Reuised manuscript received June 21, 1991 Accepted July 5, 1991
