Nonlinear analysis in process design. Why overdesign to avoid

---

I n d . Eng. C h e m . Res. 1990, 29, 805-818 Films. Ind. Eng. Chem. Prod. Res. Deu. 1980, 19, 592. Samuelson, 0.;Soderholm, I. Determination of Carbonyl Groups in Cellulose by the Hydrazine Method. Suen. Papperstidn. 1963,66, 833.

805

Tokura, S.; Nishimura, S.; Nishi, N.; Nakamura, K.; Hasegawa, 0.; Sashiwa, H.; Seo, H. Preparation and Some Properties of Variously Deacetylated Chitin Fibers. J. SOC.Fiber Sci. Technol., Jpn. 1987, 43, 288.

Secrist, R. B.; Singh, R. P. Kraft Pulp Bleaching 11. Studies on the Ozonation of Chemical Pulps. Tappi 1971, 54, 581. Soteland, N. Bleaching of Chemical Pulps with Oxygen and Ozone. Pulp Pap. Mag. Can. 1974, 75, T153.

Received for review August 28, 1989 Revised manuscript received January 2, 1990 Accepted January 29, 1990

PROCESS ENGINEERING AND DESIGN Nonlinear Analysis in Process Design. Why Overdesign To Avoid Complex Nonlinearities? Warren D. Seider,* David D. Brengel, and Alden M. Provost Department of Chemical Engineering, University of Pennsylvania, Philadelphia, Pennsylvania 19104

Soemantri Widagdo Chemistry and Chemical Engineering Department, Stevens Institute of Technology, Hoboken, N e u Jersey 07030

Overdesigns are believed to be commonplace in the chemical industry. Usually they compensate for uncertainties in the model, allow for increases in capacity, and provide margins of safe operation. In some cases, they involve elements of conservatism to protect against regimes of complex nonlinear operation that can occur in exothermic and isothermal reactors, thermally coupled and azeotropic distillation towers, supercritical extractors, and other processes, as illustrated in this paper. T h e illustrated examples have their key parameters constrained and the degree of back-mixing reduced to avoid operation near or within regimes having multiple steady states and periodic or chaotic behavior. In designing these processes, it is important to locate the range of the design parameters over which complex operating regimes may occur. Bifurcation analysis and singularity theory are potentially helpful, but these tools have not been used often. Hence, this paper briefly reviews these areas and focuses on the roles they could assume in the design of nonlinear processes. A related focus is on algorithms for flexibility and resiliency analysis and model predictive control, which when applied to nonlinear processes may permit operation closer to the steady-state economic optimum, even when near or within complex regimes. This paper argues for the coordination of design, operations, and control optimizations to reduce the instances of overdesign. Many recent advances in the analysis of nonlinear systems are leading to better modeling of processes and, in turn, improved operations. Examples of the use of bifurcation analysis and singularity theory to characterize the steady state and dynamic performance of chemical reactors are now widespread, with applications to physical processes and recycle systems gaining in importance. It is noted, however, that bifurcation analysis and singularity theory are not often used in process design. Similarly, many new applications of mathematical programming for multiobjective and multilevel optimization are permitting much improved process identification and the location of local and global optima of systems subject to nonlinear algebraic and differential constraints. These applications often arise in process design. Although examples of nonlinear programs (NLPs) abound in the literature on design, key parameters are *Author t o whom correspondence should be addressed.

often constrained, and the degree of back-mixing is reduced to keep the processes from operating near or within regimes characterized by hysteresis and periodic or chaotic behavior. This may prevent operation near steady-state economic optima. However, with new nonlinear programming strategies being developed to permit more reliable m o d e l predictive control (MPC) near or within these operating regimes, an important opportunity and challenge faces process designers: how to utilize nonlinear analysis to obtain more economical designs that are flexible and controllable in regimes characterized by greater sensitivities to modeling errors (process/model mismatch) and changes in set points and in which good servo-control and the rejection of disturbances are more difficult to achieve. With innovative applications of nonlinear programming for multiobjective design, operations, and control optimizations, it should be possible to reduce sharply the instances of overdesign in the process industries. This paper begins with a brief review of nonlinear analysis, placing emphasis on the new developments and

C 1990 American Chemical Society 08S~-SSSS/9Q~2S29-0S05$02.50/0

806

Ind. Eng. Chem. Res., Vol. 29. No. 5. 1990

on the methods being applied for process design. (Note that Seider et al. (1990) present a more comprehensive review, and Seider and Ungar (1987) describe a graduate course on this subject.) Then, it turns to the causes of overdesign, with several examples that show how it may arise. The examples focus on chemical reactors and physical separators, although overdesign in one of the process units usually significantly affects the others in the flow sheet. Finally, this paper discusses how design, operations, and control optimizations can be coordinated during the design phase to permit operation closer to the steady-state economic optimum, even when near or within complex nonlinear regimes.

Nonlinear Analysis Much of‘the design literature is concerned with solving the material and energy balances, phase equilibria, transport equations, and chemical kinetic expressions that model a process flow sheet. Often these are the equality constraints t,o be satisfied in a design optimization. Often they include both temporal and spatial derivatives. The spatial derivatives, when they exist, are normally approximated with finite differences or polynomial trial functions. When the temporal derivatives are zero, for steady-state analysis, the resulting nonlinear equations (NLEs) are solved, for the most part, by using Newtonbased methods. For global convergence and parameterization, the homotopy-continuation methods are very effective, though more expensive. Packages such as HOMPACK (Watson et al., 19871, PITCON (Rheinboldt and Burkardt, 1983), and DERPAR (Kubicek and Marek, 1983) are widely available for the implementation of Newton and fixed-point homotopies. Recently, Seader et al. (1990) introduced a “boomerang” transformation which has not failed to locate all of the solutions for any of the sets of NLEs they tested. To solve the nonlinear equations,

flxl = 0

(1)

they utilize a fixed-point homotopy _H{x,tJ= tflx)

+ (1 - t)(x_- x”) = 0

(2)

which has the advantage that only one solution to H = 0 can exist at t = 0 but the disadvantage that 5 and t can extend to +a and return from -a along the paths, as illustrated for the solution of a simple quadratic equation in Figure l a The boomerang mapping ‘2Y1

Vi’ = -

1 + y,*

Vi

(3)

where y applies for and t , keeps the new variables bounded such that -1 5 y‘ I 1, as illustrated for the quadratic equation in Figure lb. With the variable space mapped to lie within finite bounds, all of the solutions can usually be found by tracing out a single path of finite length. As homotopy-continuation algorithms follow the steady-state solution branches, they bracket limit points but normally pass over higher order singular points, such as real bifurcation points, and Hopf bifurcation points. To locate real and Hopf bifurcation points, local stability analyses along the solution branches (evaluation of the eigenvalues of the Jacobian) can be performed, or the necessary and sufficient conditions for these points can be solved (Kubicek and Marek, 1983). At a limit point, an eigenvalue becomes zero, and a change in stability often occurs. Real bifurcation points have similar properties at

o“s”:,‘./ T(X*

ROOT 1

I

- 3x + 2) + ( l - T ) ( X - X O ) ¶

xo = 1 . 5

I -10

-2

-1

1

T

a. Untransformed homotopy path

PO IMT

X ’

- - - - - - - “5

I

xo

-

1.5

-1

-I

-.5 T ’

b. Transfomed path Figure 1. Boomerang transformation for solution of Ax) = x2 - 3x + 2 = 0. Reprinted with permission from Seader et al. (1990). Copyright (1990) Pergamon.

the intersection of two or more branches. At a Hopf bifurcation point, two complex eigenvalues become pure imaginary, the steady-state branch destabilizes or remains unstable, and a branch of periodic solutions may be born. Dynamic simulations, with the temporal derivatives included, are necessary to trace out the limit cycles. As a parameter of the system varies, a branch of periodic solutions can be traced, leading to secondary bifurcations and, in many cases, chaos. These parameterization studies are very important in characterizing the nature of the attractors: steady-state nodes or foci, periodic attractors, or strange (chaotic) attractors. To trace the periodic branches, a few packages are available, including DERPER (Holodniok and Kubicek, 1984),PEFLOQ (Aluko and Chang, 1984),and AUTO (Doedel, 1986). AUTO is the most complete package for bifurcation analysis, to our knowledge, in that it has programs to locate limit points and real and Hopf bifurcation points accurately and to trace out the steady and periodic branches. Singularity and catastrophe theories have been applied to chart the steady-state solution space without the extensive path following in homotopy continuation (Golubitsky and Schaefer, 1985). The singular points are located and their normal forms and universal unfoldings are expressed so as to determine the number of steady-state solutions in their vicinities. These theories have been

Ind. Eng. Chem. Res., Vol. 29, No. 5, 1990 807 applied for reactor models (Balakotaiah and Luss, 1982; Balakotaiah, 1982) that can be reduced to a single equation and unknown (e.g., conversion). They have been used sparingly for large systems because of the difficulties in performing the Liapunov-Schmidt reduction procedure. Recently, however, Kunkel(l988) has published a promising Gauss-Newton algorithm. When the temporal derivatives are not set to zero, dynamic simulation is often performed with initial-value integrators such as LSODE (Hindmarsh, 1980) and DASSL (Petzold, 1982a,b), both available in ODEPACK. Both implement the backward difference formulas (BDFs) of Gear. LSODE is intended for the integration of ordinary differential equations (ODEs) and DASSL for differential-algebraic equations (DAEs). Both estimate the local truncation errors and adjust the step size and order of the numerical approximation. Instead of backward-difference formulas, it is often preferable to introduce polynomial trial functions in time. This is especially effective for boundary-value problems (e.g., in tracking a limit cycle with an unknown period). The method of orthogonal collocation has also been applied (Doedel, 1986) as well as the Newton-Fox algorithm in the PEFLOQ program (Aluko and Chang, 1984). Nonlinear Programming. Practical design problems are, of course, formulated with one or more objective functions, for example:

subject to

gl_x&!,u,p,4 = 0 hlx,d,u,p,tl5 0

(Pl)

where 5, d, E, and p are vectors of state and design variables, manipulated inputs, and parameters, respectively. When derivatives are present in the constraints, their approximation is a key to the solution procedure, with polynomial trial functions often utilized in a collocation method. This is particularly crucial when the equality constraints in the process model are a system of differential and algebraic (with transcendental terms) equations (DAEs). When the index of the DAEs is greater than unity, problems can arise in obtaining consistent initializations and estimates for the truncation error. After the constraints are discretized, the Lagrangian is formed and the stationarity conditions are applied, resulting in NLEs known as the Kuhn-Tucker (K-T) necessary conditions: o_f + TTVg -+ &To> = 0 g=o _h,h, = 0

vi

(Pa

where .rr and are vectors of Lagrange and Kuhn-Tucker multipliers. During the past decade, two approaches have been applied, almost exclusively, to solve the K-T conditions arising from process design optimizations. One is successive quadratic programming (SQP), as implemented in the OPT program (Biegler and Cuthrell, 1985). The other is the reduced-gradient algorithm implemented in the MINOS program (Murtagh and Saunders, 1983). When the equality constraints, g{$ = 0, exhibit multiple solutions, however, the usual opfimality conditions may not apply (Doedel, p 22, 1986). Furthermore, our experience has shown that these optimization algorithms have difficulty converging highly nonlinear equality constraints. External algorithms can be used to solve these constraints in a modular fashion (satisfying them during each iteration

of the optimization algorithm). However, these combined algorithms can have difficulty rounding limit points and can jump between solution branches, as the decision variables are adjusted, in the vicinity of multiple steady states. While this is not surprising, it supports the need to take special care when optimizing these systems and can justify the implementation of algorithms that deal explicitly with these difficulties. Convergence to a nonglobal optimum is also likely. These NLPs are nonconvex, and consequently, convergence to the global optimum cannot be guaranteed. However, decomposition algorithms are being developed to locate the global optimum with a high degree of probability. These include Benders decomposition (Geoffrion, 1972) in which complicating variables are fixed in the NLP to locate upper bounds on the minimum, and a reduced master problem is solved to adjust the complicating variables and determine lower bounds on the minimum. Strategies designed to keep the subproblems convex are also being developed (Floudas et al., 1989). Furthermore, strategies that systematically alter the combinations of active inequality constraints, so-called active set strategies, show promise (Clark and Westerberg, 1990). There are several alternate approaches to solving the K-T conditions. These are based upon the homotopycontinuation equations for solving NLEs. In one strategy, implemented in the AUTO program (Doedel, 1986), stationary points are located by using successiue continuation. In brief, a parameter w is equated to the objective function, and the homotopy path with respect to w is traced. Local extrema are located at the limit points of the continuation path. For multivariable optimization, a set of adjoint equations is derived, which constitute the necessary conditions for optimality, and a continuation strategy is designed to satisfy these conditions. However, this method does .not explicitly deal with the inequality constraints. A more general method, suggested by Vasudevan et al. (19891, replaces the complementary slackness conditions (_hixi = 0, V i) with NLEs, using a theorem by Mangasarian (1976). The advantages of this formulation are that it guarantees feasibility and the proper sign for the K-T multipliers (h I0; L 0 for a minimum, I0 for a maximum). Note that solutions to the unmodified K-T conditions (P2) may be infeasible and, depending on the signs of the multipliers, may be maxima, minima, or saddle points. Although this formulation does not guarantee location of the global optimum, continuation with the boomerang transformation may be able to locate all of the feasible local optima (maxima or minima) and, thus, the global optimum. For the optimal control problem, where the constraints are DAEs, Vasantharajan et al. (1988) utilize global spline collocation (Villadsen and Michelsen, 1978) to approximate the ODEs. This procedure is also referred to as orthogonal collocation on finite elements (Carey and Finlayson, 1975). The finite elements are adaptively located during the optimization to minimize the largest residual at the midpoints of the elements. This equidistribution method is being extended to systems having index > 3 and to adjust adaptively the number and placement of the finite elements. Process Flow Sheet Analysis. Many of the steadystate process simulators (for example, FLOWTRAN, DESIGN 11, PROCESS, and ASPEN PLUS) now perform infeasible p a t h optimization using SQP; that is, they simultaneously adjust the recycle tear variables and the design variables to satisfy the stationarity conditions. Unfortunately, however, they are of limited utility for the design of processes that op-

808 Ind. Eng. Chem. Res., Vol. 29, No. 5 , 1990

erate near or within regimes characterized by hysteresis and periodic or chaotic behavior. They have no homotopy-continuation algorithms, no algorithms for stability analysis, and no algorithms designed to locate the global optimum. SPEEDUP (Pantelides, 1988) is a system for steady-state and dynamic simulation of chemical process flow sheets. It has many convenient facilities, including a library of steady-state and dynamic models of process units that, given the structure of the flow sheet, are assembled to solve the nonlinear DAEs simultaneously. It permits the user to provide sets of DAEs, design specifications, control variable profiles, and initial conditions. Symbolic manipulation is used to prepare analytical expressions for the partial derivatives. The user can also provide preprogrammed subroutines in FORTRAN. In addition, optimization of processes in the steady state is achieved by using a feasible path strategy with SQP or an infeasible path strategy with MINOS. Its limitations are that it does not utilize continuation algorithms and has no facilities for locating singularities or performing bifurcation analysis. It would be of limited use in optimizing a process design with hysteresis and periodic or chaotic regimes because of the difficulties in converging highly nonlinear constraints, as described above. Model Predictive Control of Nonlinear Processes. Several model predictive control (MPC) algorithms were introduced about a decade ago for the control of linear systems (Richalet et al., 1978; Cutler and Ramaker, 1979; Brosilow and Tong, 1978; Garcia and Morari, 1982). These minimize a combined objective that sums the squares of the predicted errors a t sampling instants over the predictive horizon, penalized for large changes in the manipulated variables. The ODEs that model the process are expressed as an impulse-response model, which are linear equality constraints in the NLP. The resulting NLP, with inequality constraints, when solved with an SQP algorithm, usually provides much improved servo-control and rejection of disturbances, as compared with conventional analog controllers. For nonlinear processes that operate near or in regimes characterized by hysteresis and periodic or chaotic behavior, better methods are needed to approximate the ODEs across the predictive horizon. Several recent contributions (Garcia and Prett, 1986; Economou et al., 1986; Renfro et al., 1987; Patwardhan et al., 1988; Li and Biegler, 1988; Morshedi, 1986; Brengel and Seider, 1989) show great promise for much improved control near or within these regimes. Flexibility and Resiliency Analyses. The ability of a process to operate a t various levels of capacity, with different feedstocks, variations in the ambient conditions, and changes in the process parameters (e.g., fouling of the heat exchangers), is critical to a successful design. For processes that operate at steady state, with an anticipated range of operating variations, Grossmann and Morari (1984) defined a flexibility index, which is a measure of the largest deviation from the nominal conditions of design that can be tolerated without violating the inequality constraints. These authors set up an NLP to locate the range of manipulated inputs over which all of the inequality constraints are satisfied, i.e., the feasible region. They embedded this NLP in a higher level NLP to maximize the flexibility index and showed how to create a multiobjective optimization that coordinates the economic and flexibility optimizations. Their framework for flexibility analysis was formalized, and computational algor-

ithms were presented by Swaney and Grossmann (1985a,b). Grossmann and Morari also considered the capability of a process to reject disturbances and respond to set point changes on a much smaller time scale. They introduced the term dynamic resiliency to characterize the ability of a process to respond well in the face of uncertainties in the process parameters, etc. Subsequently, Palazoglu and Arkun (1986) defined robustness indexes, which represent the sensitivity and invertibility of a linearized process model, thus indicating the limitations on dynamic operability inherent in the process, independent of any particular control law. They created a multiobjective optimization with a trade-off between the economics and robustness. This was extended by Palazoglu and Arkun (1987), who defined constraints for robustness and incorporated them in the coordinated economic and flexibility optimizations of Swaney and Grossmann. Much work remains to characterize properly the limitations on the dynamic operability inherent in any given highly nonlinear system.

Overdesign There are many instances of overdesign in the chemical industry. Often equipment is oversized to allow for uncertainties in the process model and eventual increases in capacity (or changing feedstocks), to provide margins for safe operation, etc. Somewhat less obvious, perhaps, are designs that circumvent operation near or within regimes characterized by hysteresis and periodic or chaotic operation. Such overdesigns could be considered prudent, if these regimes were known or expected to exist, as protection against unreliable and unsafe operation. In this section, several examples of overdesign are cited, which can be reduced significantly through coordinated design, operations, and control optimizations. First, consider a few generalizations worthy of recognition for chemical reactors. For a given conversion of reactant, the exothermic CSTR usually has the minimum volume but displays hysteresis and periodic or chaotic operation over broad ranges of residence time (or Damkohler number). Exothermic PFTRs, on the other hand, are normally devoid of these operating regimes but require substantially more volume. In fact, to avoid back-mixing (dispersion) in packed-bed tubular reactors, it is common to increase the Peclet number for mass transfer by increasing the flow rate, thereby increasing the volume to provide sufficient residence time. Batch and semibatch configurations are even less efficient than the PFTR when the throughput is sufficiently large to justify continuous processing. Yet, batch and semibatch reactors are often used for commodity chemicals, for example, in emulsion polymerization (to be discussed) and detergent manufacture. A second generalization concerns feed ratios in excess of the stoichiometric requirements for the principal reactions. It is best, of course, to convert stoichiometric amounts of the feed chemicals completely. When this is impractical, the most valuable reactants are often consumed in the presence of excess quantities of the other reactant(s). Overdesign is associated with large excesses because reactor sizes and the costs of separation and recirculation increase. However, for many reaction systems (e.g., the allyl chloride process below), smaller excesses lead to regimes characterized by hysteresis and periodic or chaotic behavior. Although operation near or within these regimes can be more attractive economically, larger excesses are often utilized to circumvent potential operating problems.

Ind. Eng. Chem. Res., Vol. 29, No. 5, 1990 809

c

L

W

I

c

LL

“ 2

50

AN Oa = Vk(TfVq Figure 2. Seven kinds of solution diagrams for A B. Reprinted with permission from Balakotaiah and Luss (1983). Copyright 1983 Pergamon.

-

0.15

0.25 DAMKOHLER NUMBER

0.35

--

Figure 4. Typical solution diagram for A B C. Solid curves are stable steady-state solutions; dashed curves are unstable. Dotted curves are branches of limit cycles (stable and unstable) between four Hopf bifurcation points. Reprinted with permission from Doedel and Heinemann (1983). Copyright 1983 Pergamon.

vicinity of the winged-cusp singular point (see Figure 3). When Hopf bifurcation points and periodic oscillations were included, the studies of Uppal et al. (1974) and others identified an additional 16 kinds of solution diagrams. Razon and Schmitz (1987) then examine the exothermic reactions A-B-C

“0

2 y

4 = E/ RT,

6

8

Figure 3. Solution diagrams in the vicinity of the winged-cusp singular point (see Figure 2). Reprinted with permission from Balakotaiah and Luss (1983). Copyright 1983 Pergamon.

Before proceeding, an underlying premise may be stated: less overdesign saves money but requires more complex operations. However, when good models, often developed in conjunction with pilot plant studies, are incorporated in model predictive controllers for nonlinear processes, much improved disturbance rejections and set point responses are obtained. These MPCs permit operation closer to the economic optimum and in regimes that hitherto were difficult (or impossible) to achieve. The next subsection considers overdesigns in exothermic reactors. First, the breadth of operating regimes associated with the general reactions A .- B, A B C, and A B, A C, is considered before specific examples of overdesign are presented. Exothermic Reactors. In an excellent review, Razon and Schmitz (1987) consider the “diabatic” (nonadiabatic) CSTR. For the classic reaction,

--

- -

A-B they review the work of Balakotaiah and Luss (1982), who used singularity theory to identify a third-order, wingedcusp singularity and displayed seven kinds of steady-state solution diagrams with the conversion, X, plotted as a function of the Damkohler number, Da,holding the other parameters fixed (see Figure 2 ) . They observed that a maximum of three steady states can occur a t a given Da and constructed a map of the solution diagrams in the

and review the work of Farr and Aris (1986),who identified a maximum of 7 steady-state solutions and 23 kinds of solution diagrams. In this regard, they note that “the complexity and number of bifurcation diagrams increases with the complexity of the system.” When extended to include periodic oscillations, Doedel and Heinemann (1983) obtained solution diagrams with two periodic branches (see Figure 4). Turning next to the exothermic reactions A-B -C Balakotaiah and Luss (1983) identified a maximum of 5 steady-state solutions and 48 kinds of solution diagrams, leading Razon and Schmitz to reflect “the complexities in dealing with higher dimensional state spaces present formidable obstacles to further progress.” A similar study was carried out by Jensen and Ray B in diabatic, (1982) for the exothermic reaction A packed-bed, tubular reactors, using a pseudohomogeneous model for dispersion. These authors constructed maps of the solution diagrams as the dimensionless heat of reaction and dimensionless heat-transfer coefficient vary, at several values of the Peclet numbers for heat and mass transfer, where Pe = Peh = Pe,. At high Pe (as Pe m) with no dispersion, plug flow occurs in the absence of hysteresis and periodic operation. At intermediate Pe = 5 , 14 kinds of solution diagrams were computed, with a maximum of 5 steady states and a maximum of 3 Hopf bifurcation 0, complete points. For the other extreme, as Pe back-mixing occurs (as in a CSTR). Here, only six kinds of solution diagrams were computed, with a maximum of three steady states and a maximum of two Hopf bifurcation points. Having reviewed, in general, the causes of hysteresis and periodic operations in CSTRs and tubular reactors, how overdesign arises in several specific processes will now be considered. ( 1 ) Allyl Chloride Process. Allyl chloride (AC) is

-

-

-

810 Ind, Eng. Chem. Res., Vol. 29, No. 5 , 1990 Table I. Heats of Reaction and Kinetics for Allyl Chloride Reactions (Biegler and Hughes, 1983) Wb ko, Ib-mol/ reaction Btu/lb-mol (h fts atm? EA/R, OR I -4 800 206 000 13 600 1 -79 200 11,; 3 430 21 300 3 -91 800 4.6 X lo8 1

-L ,::: Figure 5. Allyl chloride process. Results are for the optimal design of Vasantharajan and Biegler (1988)

produced commercially by the chlorination of propylene (Fairbairn et al., 1947): CHz=CHCH3

+ Cb

ki

CHz=CHCHzCI + Clo

Jk CH3CHCICHfiI (KP)

CHCl=CHCH&I (WP

-1

+ HCI

+ HCI

All three reactions are exothermic, with second-order kinetics and an Arrhenius dependence on temperature. Rate expressions, kinetic constants, and heats of reaction are provided in Table I. Vasantharajan and Biegler (1988) adjusted T, the reactor inlet temperature, 7, the propylene-to-chlorine ratio fed to the reactor, L , the reactor length, P,, the compressor outlet pressure, and four other design variables to maximize the net sales for the flow sheet in Figure 5 (illustrated for the optimal design). Propylene at 1.38 MPa (200 psia) and 338.7 K (150 O F ) is compressed to P, = 1.91 MPa (277 psia) and mixed with recycled fluid. Chlorine at 0.69 MPa (100 psia) and 299.8 K (80 O F ) is also compressed to P,. Heat is added to the propylene, which is mixed with the C1, stream to provide a reactor feed at T = 727.6 K (850 OF, upper bound). The reactor is modeled as a PFTR, based upon pilot plant data (Fairbairn et al., 1947). The propylene-to-chlorine ratio in the feed is high, 17 = 4.56, and they note that this “prevents the excess loss of allyl chloride to DCP=” and “helps to reduce coke buildup” in the homogeneous reactor. The reactor products are cooled to 394.3 K (250 OF) and compressed to permit condensation with cooling water in the depropenizer. The propylene stream is scrubbed to recover HCl, dried with caustic, and recirculated. Both DCP and DCP= are removed from the allyl chloride product in a finishing column. A similar study had previously been performed by Biegler and Hughes (1983) using a CSTR model to represent commercial reactors that promote back-mixing. Biegler and Hughes note that “they are really flame reactors, although the residence times are longer-about 0.1 sec.” In their study, the maximum net sales was obtained with the reactor temperature T = 699.8 K (800 O F , lower

Table 11. Heats of Reaction and Kinetics for Maleic Anhydride Reactions (Wohlfahrt and Emig, 1980) AHR(848K), ko, m3/ reaction kJ/mol (kg cat. s) EA,J l m ol 1 -1850 4 280 1 2 660 2 -3274 70 100 15 000 3 -1423 26 10 800 ~

I

_

_

bound) and 7 = 5.86. In an attempt to reduce the overdesign, we examined the performance of an adiabatic CSTR (minimum volume) as 7 is reduced. The results of parameterization using homotopy continuation at 7 = 1.5 and To= 450-600 K are particularly revealing. Figure 6a shows the variation of the reactor temperature, T , with the inverse residence time, LY = (RT/PJ(F/V), where Pf is the feed pressure (MPa), F the feed flow rate (kmolis), and V the reactor volume (m3). Note the Hopf bifurcation points (a). Figure 6b shows the unusual variation in the partial pressure of allyl chloride, Figure 6c the fraction of propylene converted, Figure 6d the fraction of chlorine converted, and Figure 6e the ratio of allyl chloride to DCP= partial pressures. With 17 = 1.5, high partial pressures of allyl chloride are accompanied by low conversions of Clz. This bifurcation analysis suggests that 7 = 1.5 is too low, with higher ‘7 necessary to maximize the venture profit. In a subsequent optimization (Brengel and Seider, 19901, the maximum venture profit was obtained with 7 = 3.152, TO = 483.3 K, T = 632.1 K, CY = 1.019, xc3= = 0.251, xci2 = 0.9986, PAC = 0.0831 MPa, p A C / p D C p = = 1.663, p A C / p D C p = 1.444, and Pf = 1 MPa. When compared with the Biegler and Hughes design, far less propylene is recirculated and the costs for separation are sharply reduced. Furthermore, a venture profit of $5.38 X 107/year is recorded, as compared with a venture loss for the previous design (Biegler and Hughes used an objective function that did not include capital and utility costs, which are appreciable for the larger recirculation rates). It should be noted that the costs of separation and recirculation were estimated to be proportional to the propylene flow rate raised to the 0.6 power (more detailed design calculations are currently underway). The new design has a unique steady state but operates closer to the regimes involving hysteresis and periodic behavior illustrated in Figure 6. A possible alternative for further reducing the overdesign is to recirculate propylene, chlorine, and HCl. Smaller 7 can be expected, with smaller depropenizer costs but higher costs for recirculation and reactor volume. A purge stream, followed by a scrubber and drier, would be necessary, and the purge-to-recycle ratio would be an important design variable. (2) Benzene to Maleic Anhydride. A similar process was reported by Westerink (1987) involving the partial oxidation of maleic anhydride over a Vz05 catalyst in a fluidized bed reactor. Like the allyl chloride process, two reactions in series and a single reaction in parallel are postulated: ki

C6H6

+ 91202

---.)

O=C

co-l C=O

+ 15 /202

I

(MA)

+ 302

k3

6C02

1 =cH

+

3H20

1. 4C02

+

H2O

+ 2C02 + 2H20

_

e.0

Ind. Eng. Chem. Res., Vol. 29, No. 5, 1990 811

='O

4.0

0.0

1-

3 :.

i, :.. .:

.. . .... 7.0

h

3

a.0

0.0

J

1.0

r

M P

2 I-

P

. ... ., ',

,

.... :. . ...

0.5

1.0

.

P

1 .O

5.0

4.0

'

I

I

1

I

I

1

0.0

0.5

1.0

9.5

1.0

1.5

0.0

0.0 0.0

b.

a (. 10-3)

a. 0.4

1.5

1.0

1.5

3.0

a (. 10')

1 .O

0.0

0.3

0..

L-

x

0.1

ii

..

x

..

0.4

0.1

0.1

0.0

0.0

I

I

I

I

0.9

1.0

1.5

1.0

a (,

C.

s'o

I

1.5

1

0.0

103)

1

4.0

a.0

k

p

1.0

aY

0.0

-1.0

e.

a (. 183)

Figure 6. Solution diagrams for the adiabatic CSTR to produce allyl chloride. bifurcation points ( 0 ) .

With a large excess of air (benzene is 1mol YO of the feed), the reactions are reported to be first-order and exhibit an Arrhenius dependence on temperature. Rate constants and heats of reaction are provided in Table 11. For the design of this reactor, Westerink applied an algorithm that seeks the maximum yield within the region of "safe unique operation". A series of constraints were derived that define this region, and the design parameters

7=

1.5. Stable

(-1 and unstable

(-e)

steady states and Hopf

were adjusted to avoid steady-state multiplicity. For a packed-bed tubular reactor, Westerink derived criteria to obtain an integral yield that lies within a desired fraction of the maximum yield in an isothermal reactor. Overdesign is inherent in the Westerink approach, which is most noteworthy because its explicit objective is to protect against operation within complex operating regimes. (3) Propylene Oxide to Propylene Glycol. Fogler

812 Ind. Eng. Chem. Res., Vol. 29, No. 5 , 1990

(1986) suggestas that the hydrolysis of propylene oxide + H20

CHZ-CHCH~

\'

-

CH&HCH,

(1

AH AH

be carried out in a CSTR in the presence of 0.1 wt % H2S04(catalyst) and methanol (presumably to prevent phase splitting). In excess HzO, Furusawa et al. (1969) found the rate of this exothermic reaction to be first-order with respect to propylene oxide (ko= 16.96 X 10l2h-l, EA = 32 400 Btu/lb-mol). Fogler considers a design with a single, high-conversion steady state. As the volume (residence time) is decreased, multiple steady states appear. Here again, design optimization can lead to more complex operating regimes. (4) Other Exothermic Reactions. For the most part, catalytic reactions such as

CO + 2H2 CO

+ HzO

Cu-Zn-Al

\

2 7

0 r

C02

-

co + y202 Pt

t

CHSOH

Cu-Zn-A1

(1

Figure 7. Schematic of fermentation process.

2 0

i

+ Hz a . o i

COZ

4.0

(4)

where p is the specific growth rate (maximum at S = K ) and S is the substrate concentration. In the absence of a model for the mass transfer of oxygen, a variable yield coefficient

y(SI = p{S]/a{S)= a

+ bS

,

,

4.0

5.2

:

,

,

,

,',

S.0

6.0

5.4

0.0

7.2

, 7.0

8.0

OoilO..

are carried out in fixed-bed reactors. When catalyst attrition is not an important problem, the alternative of fluidization can provide more efficient heat removal (eliminating hot spots), eliminate the periodic maintenance of fixed beds (e.g., catalyst replacement or decaking), etc. It can be expected to become an important alternative, as the new model predictive controllers permit more reliable operation in the vicinity of multiple steady states. Isothermal Reactors. In isothermal reactors, effects other than the rate of heat generation are responsible for hysteresis and periodic or chaotic regimes. Often the reactions involve adsorption at catalyst sites, with rate inhibition at higher concentrations of the reactants. In homogeneous systems, autocatalytic reactions often lead to periodic and chaotic operation (see, for example, Gurel and Gurel(l983) and Razon and Schmitz (1987)). Enzyme reactions are commonly observed to display these regimes as well as hysteresis. To circumvent these regimes, for these and other systems, designers may resort to overdesign. (1) Aerobic Fermentations. The aerobic fermentation of Saccharomyces cereuisiae on a sugar-cane molasses medium was carried out in a laboratory fermentor by Borzani et al. (1977), who demonstrated the existence of periodic regimes in which the cell mass builds until limited by the available oxygen and then is depleted until the oxygen becomes available in excess. Agrawal et al. (1982) modeled this cyclical behavior using a one-hump growth rate: p{S\ = kSe-S/K

, 4.4

(5)

was found to be necessary but not sufficient to reproduce

-0 5 0 0

0 3

o s

09

1 2

' 5

30

Figure 8. Steady state [(-) stable, unstable] and periodic (0) attractors of cell mass as a function of the Damkohler number for the fermentation reaction with y = 0.4 and = 0.1. (e..)

the periodic behavior. Note that (T is the substrate consumption rate. The chemostat in Figure 7 was modeled as a constant volume ( V )CSTR, with cell and substrate mass balances: dX - = -;x + pL(S)X dt (7)

These mass balances were dedimensionalized to yield dC1 dr

-=

-C,

+ Da(1

-

C2)e4/~C1

(8)

dC2 where C1 and C2 are the dimensionless cell mass [X/ (SFWSF))] and substrate conversion [(SF - s ) / s F ] , p = a / ( b S F ) , y = K/SF, 7 = t ( F / V ) , and Da is the ratio of the specific growth rate to the inverse of the residence time [p(SF]/(F/V)I. To identify the regimes of hysteresis and periodic operation, Agrawal et al. (1982) varied the parameters Da,

Ind. Eng. Chem. Res., Vol. 29, No. 5, 1990 813 Monomrr W,O rmdrif irr initiator

0.10 0.08

0.08

X 0.04 0.02

0.00

i

2

i

6

8

to

e ( min ) 1.0

.

X 0.4

0.2

0

Figure 9. Emulsion polymerization reactor. Reprinted with permission from Rawlings and Ray (1987). Copyright 1987 Pergamon.

P, and

y. Figure 8 is comparable to the results they obtained. Note that the stable steady state that yields the highest cell mass occurs at the Hopf bifurcation point, where Da = 0.6170. Beyond that, the steady-state attractors are unstable and the periodic attractors are stable with large periods, corresponding to the experimental observations. For the design optimization of a flow sheet containing the fermentation reactor, the maximum venture profit was located a t Da = 0.5597 (y = 0.4, = O.l), an unstable steady state (Brengel and Seider, 1990). However, this study confirmed that a model predictive controller for nonlinear processes (Brengel and Seider, 1989) can reject disturbances and provide satisfactory servo-responses at this steady state. A conservative overdesign would be at higher Da (higher volume and lower dilution rate), which would provide less cell mass and greater residence times. Note that, at the economic optimum, the venture profit exceeds that at the Hopf bifurcation point by 31.6%. It is also noteworthy that many industrial designs, especially in the pharmaceuticals industry, are in batch and fed-batch reactors, often to circumvent the perceived operating problems, even when the throughput is sufficient to justify continuous operation. (2) Emulsion Polymerization Reactions. The work of Penlidis et al. (1989) describes a classic case of overdesign in which batch or semibatch reactors have been used for the majority of industrial latex production. Although continuous reactors require smaller volumes, provide better heat removal, and offer potentially more consistent product quality, batch and semibatch reactors are often used to avoid the limit cycles that occur in continuous operation. Rawlings and Ray (1987) present a population balance model that quantitatively matches the experimental observations for the polymerization of methyl methacrylate. When oscillation occurs, it is explained by the following mechanism. Radicals from the monomer droplets (initiated with potassium persulfate; see Figure 9) enter into the micelles of emulsifier (e.g., sodium dodecyl

20

40

80

80

e ( min ) Figure 10. Steady-state conversion versus residence time. (0) Data of Schork, ( 0 )Hopf bifurcation, (1) dynamic simulation point, (- - -) unstable steady states, Sf= 0.03 mol/L, If = 0.03 mol/L, T = 40 "C, V , = 0.70. Reprinted with permission from Rawlings and Ray (1987). Copyright 1987 Pergamon.

sulfate) and initiate the growth of polymer particles. As the particles polymerize and grow, more emulsifier is required to stabilize the increased surface area until the emulsifier is consumed, after which no new particles are initiated. Gradually, the larger particles wash out of the reactor, and the micelles reappear as more emulsifier enters in the feed stream. The results of Rawlings and Ray, as compared with the data of Schork (1981) are illustrated in Figure 10 and show stable steady states at low residence times, d 5 5 min. A t d = 5 min, Hopf bifurcation occurs, and the steady states become unstable. Multiple unstable steady states are observed for 27 I d I 56 min, with mildly unstable oscillations along the high conversion branch, observed both experimentally and theoretically. Faced with the problem of oscillation, Penlidis et al. (1989) present several design alternatives: (1)overdesign with batch and semibatch reactors; (2) use of large excess quantities of soap (emulsifier) to eliminate (or dampen) the limit cycles, an unfavorable option due to cost, product contamination, and the generation of large numbers of smaller particles; (3) use of an advanced control scheme (linear quadratic stochastic optimal control) which, as can be expected, does not eliminate the oscillations but could stabilize the high-conversion steady states; and (4) seed a CSTR train using a very small PFTR, designed to initiate particle growth in small residence times, as illustrated in the split-feed system shown in Figure 11. This latter design eliminates the undesirable oscillations and is very flexible in controlling the particle size and conversion. Both theoretical and experimental studies (Penlidis et al., 1989) demonstrate that the split-feed system can be expected to provide an excellent design, sharply reducing the overdesign inherent in alternatives 1 and 2. A more careful comparison with alternative 3 is needed, in view of the

814

Ind. Eng. Chem. Res., Vol. 29, No. 5 , 1990 '4

w

Figure 11. Latex reactor train configuration. Reprinted with permission from Penlidis et al. (1989). Copyright 1989 Pergamon.

recent advances in model predictive controllers for nonlinear systems. (3) Ethylene Hydrogenation. The hydrogenation of ethylene occurs in an exothermic reaction:

C2H4 + H,

+

a. Tower configuration

C2H6

which was demonstrated experimentally by Furusawa and Kunii (1971) to proceed over a composite catalyst of CuO, Cr,O!, and MnO, with bimolecular kinetics (involving two reaction sites) and the intrinsic rate of reaction:

Their study verified that ignition and extinction accompany steady-state multiplicity in a single particle of catalyst. This occurs as pE varies about a critical partial pressure, above which the kinetics approach first-order and below which they approach negative first-order. Earlier, Matsuura and Kat0 (1967) examined theoretically the stability of well-stirred, isothermal reactors, which are comparable to well-cooled fluidized beds for solidcatalyzed reactions. They postulated bimolecular kinetics and studied the multiple steady states which may occur when a large excess of hydrogen relative to ethylene is fed to the reactor. They also studied the impact of dispersion in well-cooled, packed-bed tubular reactors and observed that the space-time conversion is greatly increased in a CSTR. Although the results of design optimization are not available, it seems clear that the fluidized bed is a promising alternative and that the hydrogen-to-ethylene ratio in the feed is a key design variable. Overdesign with packed-bed reactors using large excesses of hydrogen to ethylene should be far from optimal. Thermally Coupled Distillation Towers. In two important papers, Stupin and Lockhart (1972) and Tedder and Rudd (1978) successfully showed the various arrangements of interlinked distillation towers that can separate ternary mixtures into three products more economically than two towers in series. Recently, however, the potential for operating problems was demonstrated by Chavez C. et al. (19861, who obtained multiple steady-state solutions of the MESH equations for three configurations, one of which is illustrated in Figure 12a. For this Petlyuk system, they computed four sets of interlink flow rates over a sizable range of the reflux ratio. Calculations were performed for the separation of a benzene-toluene-oxylene system, with the following specifications: 95% molar purity of benzene in the distillate, 90% toluene in the middle product, 95% o-xylene in 380 kmol/h of the bottoms stream, and the reflux ratio [4.525,5.75]. As shown in Figure 12b, they computed a significant displacement

,,t

I

b. Interlink liquid flow rate as a function of the reflux ratio. Figure 12. Petlyuk system. Reprinted with permission from Chavez C. et al. (1986).

between t,he flow rates of the liquid interlink streams. However, the flow rates of the distillate and middle product, the boil-up, and the condenser load are negligibly different a t each reflux ratio. It is noteworthy that Lin et al. (1987) subsequently located all of the steady-state solutions from a single starting point using homotopy continuation. Unfortunately, the studies of Chavez C. et al. and Lin et al. did not include checks on the stability of the solutions and provided no experimental verification of the multiple solutions. For design calculations (to size the towers properly), it is important to recognize that multiple solutions are possible, to locate the solutions when they exist, to check for stability, and to perform dynamic simulations. Concerning overdesign, configurations of interlinked towers are not widely used in industry, except for crude petroleum tower systems. Perhaps this can be attributed to perceived operating problems in achieving product purities in the vicinity of the multiple steady states and the lack of confidence in the steady-state models. Instead, more expensive, conventional, two-tower sequences are often installed. Heterogeneous Azeotropic Distillation Towers. For a tower to dehydrate sec-butyl alcohol (SBA) with disec-butyl ether (DSBE; see Figure 13a), Kovach and Seider (1987a,b) showed the sensitive relationship between the aqueous reflux ratio and the interface between trays with

Ind. Eng. Chem. Res., Vol. 29, No. 5, 1990 815

n SBA

Organic Phase

."

DSBE

n

1

Figure 14. Aniline process flow sheet.

REBOILER SEA

a. Schematic of azeouopic distillation tower

didates for model predictive control algorithms. When these algorithms are perfected, it is expected that future towers could be built with far less overdesign. Phase Splitting for Product Recovery. Clark (1990) presents an excellent example in which the ability to split away a second liquid phase in a flash vessel can sharply reduce the downstream distillation costs. In this example, Clark analyzes the flow sheet in Figure 14 to hydrogenate nitrobenzene to aniline according to the reaction CeH5N02 + 3H2 CeH5NHZ 2Hz0

-

TOP

Ll

B

7

z h

ld

Ll

H

BOT Avg. Liq. Mole Fraction b. Liquid mole fractions for 33 and 12 theoretical trays. -single liquid phase, ..... two liquid phases.

Figure 13. sec-Butyl alcohol dehydration.

two liquid phases (LL) and trays with just one liquid phase (L). Their calculations were extended by Widagdo et al. (19891, who demonstrated that the tower can obtain comparable purity and recovery of SBA in the bottoms product with as few as one-third of its trays. This is illustrated in Figure 13b, where the liquid mole fractions are plotted for the industrial tower (33 theoretical trays) and for a hypothetical tower with 12 theoretical trays. Two points are worthy of emphasis. First, steady-state multiplicity is computed as a second liquid phase is added to each of the trays. To avoid problems of convergence near the limit points, a homotopy continuation algorithm was utilized for solution of the MESH equations. Second, the 33-tray tower permits the operator to respond to disturbances, as the L-LL phase interface moves through the stripping section. Much more rapid responses would be necessary for the 12-tray tower to prevent the interface from dropping through the bottom of the tower and contaminating the bottoms product with sizable amounts of water and DSBE. Such towers appear to be prime can-

+

The vapor from the flash vessel is recycled, and the liquid phase is sent to a distillation train. Note that, if the temperature and pressure of the flash vessel are adjusted properly, a water phase can be recovered, sharply reducing the cost of the first distillation operation. Such a design would require an advanced control strategy to maintain the second liquid phase in the flash vessel when disturbances are encountered in the feed composition, temperature, etc. In this process, overdesign can be expected in the dehydration tower to protect against control problems in the flash vessel. Supercritical Extraction. In supercritical extraction (SCE), a solvent such as carbon dioxide is compressed to just beyond its critical point and used to extract a valuable solute from a liquid stream, a slurry, or a solid. For the most efficient designs, the solute is recovered in a flash vessel at pressures just below the critical pressure of the solvent. Then, the solvent can be recompressed at low cost (low compression ratio) and recycled to the extractor. However, small changes in temperature and pressure in the critical region can lead to abrupt phase changes. For example, let the phase envelope in Figure 15a be at the composition of the feed to the flash separator, and note the abruptness of the phase changes in the retrograde region. Small changes in P and T can shift the flash vessel out of the two-phase region. Furthermore, in some cases, a branch of three-phase equilibria occurs, and a second liquid phase is encountered. Consider also the possibilities for steady-state multiplicity in the retrograde region. Figure 15a shows that two solutions exist for pressures specified between P, and P,. Similarly, there are two solutions for T ' I T I T-. Cygnarowicz and Seider (1989) showed that, at high solvent/feed ratios and at a given extractor pressure (>P, for COJ, the same amount of acetone can be extracted from water at different temperatures (see Figure 15b). These two solutions resulted in local and global minima of the utility cost for the SCE loop to dehydrate acetone. Carbon dioxide is an inert, nontoxic solvent and has great potential for the recovery of valuable solutes, such as @-caroteneand caffeine, in food processing (Cygnarowicz and Seider, 1990). Unfortunately, its use in SCE processes presents delicate phase-control problems. However, it may

816 Ind. Eng. Chem. Res., Vol. 29. No. 5 , 1990 L

I 1

Mcdcl Evoluuon

Mor Plant Data

W

a z

comdetc

VI VI W

Then, d y ” c

for sm-up. scrvo - and rcgulawr conrml

a n

1 TEMPERATURE

a. Phase envelope at constant composition

Figure 16. Role of model evolution and pilot plant data in design and control optimizations.

b. Moles of acetone exnacted at constant

solvenvfeed ratio (=8) Figure 15. Supercritical extraction. Reprinted from Cygnarowicz and Seider (1989).

be possible to overcome these with model predictive control algorithms. In this case, overdesigns involve extraction at higher pressures and solute recovery at lower pressures to avoid operation in the critical region.

New Design Tools Unfortunately, the generalized design tools available today do not permit the design engineer to study complex operating regimes easily. The examples above show that better tools are needed to explore operation near or within these regimes. More complete models are evolving, and design optimizations are gaining in complexity (more local solutions) as the design variables shift toward these regimes. Perhaps of even greater significance, extensions of these models for dynamic simulation are contributing to the implementation of improved MPCs which provide better control than previously possible with PID and other simple controllers. To develop more efficient processes, with less overdesign, better strategies are needed. For high-performance processes, it becomes important to screen and reject designs and operating regimes that are not sufficiently controllable. To accomplish this, the process and the MPC should be designed in parallel. This suggests that new design

strategies be developed in which the models developed for optimization, often in concert with pilot-plant studies, are extended and used, in the design stage, to predict their controllability. Figure 16 illustrates how evolution of the model, with the assistance of pilot-plant data, can permit design optimization to proceed in parallel with control optimization. The paragraphs that follow provide clarification. Process design usually begins with the analysis of simple structures using approximate models. Gradually, as the synthesis tree is pruned, more complex models are justified to represent the real phenomena more accurately. However, the more complex models usually have a richer solution space, and the number, type (steady, periodic, or chaotic), and stability of the solutions varies with the specifications and the parameters of the model. In design calculations, as in all mathematical modeling, the engineer must learn to beware of algorithms that converge to solutions that are not physically correct. Often, in the design stage, experimental data are not available and the design may continue based upon the wrong solution. For this reason, it is important to develop methods that can easily determine the number and kinds of solutions that exist (e.g., continuation methods). Then, when a rich solution space is encountered, steps can be taken to verify and refine the model experimentally through pilot-plant studies. Refinement of a model often involves the relaxation of assumptions, the addition of terms, and the reassessment of the parameters. In this way, the process/model mismatch is gradually reduced, although not eliminated. Design optimization can proceed in an attempt to maximize an objective function (e.g., profitability), penalized for poor controllability. To determine a measure of controllability, control optimization can proceed in parallel by using a dynamic model of the process. Eventually, the pilot plant can be outfitted with improved control systems, such as MPCs, to confirm that the most promising designs provide the projected quality of servo-response and disturbance rejection. Coordination of Economic, Flexibility, and Control Optimizations. To reduce overdesign. it seems clear that

Ind. Eng. Chem. Res., Vol. 29, No. 5 , 1990 817 the computing tools for homotopy continuation, bifurcation analysis, nonlinear programming, and dynamic simulation should be applied, preferably in a coordinated fashion. Beyond interfaces between the principal computer programs, to simplify the flow of data and computed results, it would be desirable t~ coordinate the economic and control optimizations, as well as the optimization that maximizes the flexibility of the process (in the face of uncertain specifications and parameters). Such a coordinated computing system, entitled PRODOC (Process Design, Operations, and Control), is being developed and will be described in a later paper (Brengel and Seider, 1990). Conclusions 1. Nonlinear processes may have steady-state economic optima near or within regimes characterized by hysteresis and periodic or chaotic behavior, as illustrated by the many examples in this paper. These examples support the contention that industrial practitioners often install conservative designs, in part, to avoid difficult operating regimes. Some are skeptical that processes can operate more economically near or within these regimes and would improperly classify such process designs as pathological. 2. To improve the characterization of these operating regimes, design engineers could apply the tools for nonlinear analysis such as algorithms for homotopy continuation, bifurcation analysis, and dynamic simulation, coupled with those for mathematical programming. 3. It is suggested that the design of model predictive controllers be coordinated with the generation of the process model and the optimization of the process economics. In this way, more economical processes can be generated that can be controlled effectively, even though they operate in more complex regimes of operation. Acknowledgment Partial funding was provided by the Design Theory and Methodology Program of the NSF under Grant DMC8613484 and is gratefully acknowledged. L i t e r a t u r e Cited Agrawal, P.; Lee, C.; Lim, H. C.; Ramkrishna, D. Theoretical Investigations of Dynamic Behavior of Isothermal Continuous Stirred Tank Biological Reactors. Chem. Eng. Sci. 1982,37,453. Aluko, M.; Chang, H.-C. PEFLOQ: An Algorithm for the Bifurcation Analysis of Periodic Solutions of Autonomous Systems. Comput. Chem. Eng. 1984, 8 (6), 355. Balakotaiah, V., Structure of the Steady-state Solutions of Lumped-parameter Chemically Reacting Systems. PhD Thesis, Univ. of Houston, 1982. Balakotaiah, V.; Luss, D. Structure of Steady-state Solutions of Lumped Parameter Chemically Reacting Systems. Chem. Eng. Sci. 1982, 37 ( l l ) ,1611. Balakotaiah. V.; Luss, D. Multiplicity Features of Reacting Systems: Dependence of the Steady States of a CSTR on the Residence Time. Chem. Eng. Sci. 1983, 38, 1709-1721. Biegler, L. T.; Cuthrell, J. E. Improved Infeasible Path Optimization for Sequential Modular Simulators-11: The Optimization Algorithm. Comput. Chem. Eng. 1985, 9 (3), 257. Biegler, L. T.; Hughes, R. R. Process Optimizatioh: A Comparative Case Study. Comput. Chem. Eng. 1983, 7 (5), 645-661. Borzani, W.; Gregori, R. E.; Vairo, M. L. R. Some Observations on Oscillatory Changes in the Growth Rate of Saccharomyces Cerevisiae in Aerobic Continuous Undisturbed Culture. Biotechnol. Bioeng. 1977, 19, 1363. Brengel, D. D.; Seider, W. D. Multistep Nonlinear Predictive Controller. Ind. Eng. Chem. Res. 1989, 28, 1812-1822. Brengel, D. D.; Seider, W. D. Coordinated Design and Control Optimization of Nonlinear Processes. In preparation, 1990. Brosilow, C. B.; Tong, M. Inferential Control of Processes: Part 11. The Structure and Dynamics of Inferential Control Systems. AIChE J . 1978, 24 ( 3 ) , 492.

Carey, G. F.; Finlayson, B. A. Orthogonal Collocation on Finite Elements. Chem. Eng. Sci. 1975, 30, 587. Chavez C., R.; Seader, J. D.; Wayburn, T. L. Multiple Steady-state Solutions for Interlinked Separation Systems. Ind. Eng. Chem. Fundam. 1986, 25, 566-576. Clark, P. A. Bilevel Programming for Steady-state Chemical Process Design: 11. Performance Study for Design Calculations. Comput. Chem. Eng. 1990, 14 (l),99-109. Clark, P. A.; Westerberg, A. W. Bilevel Programming for Steadystate Chemical Process Design: I. Fundamentals and Algorithms. Comput. Chem. Eng. 1990, 14 ( l ) ,87-97. Cutler, C. R.; Ramaker, B. L. Dynamic Matrix Control-A Computer Control Algorithm. AIChE 86th National Meeting, Houston, 1979. Cygnarowicz, M. L.; Seider, W. D. Effect of Retrograde Solubility on the Design Optimization of Supercritical Extraction Processes. Ind. Eng. Chem. Res. 1989, 28, 1497-1503. Cygnarowicz, M. L.; Seider, W. D. Design and Control of a Process to Extract Beta Carotene with Supercritical Carbon Dioxide. Biotechnol. Prog. 1990, in press. Doedel, E. J. AUTO: Software for Continuation and Bifurcation Problems in Ordinary Differential Equations; Concordia University: Montreal, 1986. Doedel, E. J.; Heinemann, R. F. Numerical Computation of Periodic Solution Branches and Oscillatory Dynamics of the Stirred Tank Reactor with A B C Reactions. Chem. Eng. Sci. 1983,38 (9), 1493. Economou, C. G.; Morari, M.; Palsson, B. 0. Internal Model Control. 5. Extension to Nonlinear Systems. Ind. Eng. Chem. Process Des. Deu. 1986, 25, 403. Fairbairn, A. W.; Cheney, H. A.; Cherniavsky, A. J. Commercial Scale Manufacture of Allyl Chloride and Allyl Alcohol from Propylene. Chem. Eng. Prog. 1947,43 (6), 280-290. Farr, W. W.; Aris, R. Reflections on the Multiplicity of Steady States of the Stirred Tank Reactor. Chem. Eng. Sci. 1986, 41 (6), 1385-1402. Floudas, C. A.; Aggarwal, A.; Ciric, A. Global Optimum Search in Nonconvex NLP and MINLP Problems. Comput. Chem. Eng. 1989, 13 (lo), 1117-1132. Fogler, H. S.Elements of Chemical Reaction Engineering; Prentice-Hall: Englewood Cliffs, NJ, 1986. Furusawa, T.; Kunii, D. Study of Multiple Steady States in Single Pellet of Catalyst. J . Chem. Eng. Jpn. 1971, 4 (3), 274-280. Furusawa, T.; Nishimura, H.; Miyauchi, T. Experimental Study of a Bistable Continuous Stirred-tank Reactor. J . Chem. Eng. Jpn. 1969, 2 (11, 95-100. Garcia, C. E.; Morari, M. Internal Model Control. 1. A Unifying Review and Some New Results. Ind. Eng. Chem. Process Des. Dev. 1982, 21, 308. Garcia, C. E.; Prett, D. M. Advances in Industrial Model Predictive Control. In Chem. Proc. Cont.-CPC III; Morari, M.; McAvoy, T . J., Eds.; CACHE, Elsevier: New York, 1986; p 245. Geoffrion, A. M. Generalized Benders Decomposition. J . Optim. Theory Appl. 1972, 10 (4), 237-260. Golubitsky, M.; Schaefer, D. G. Singularities and Groups in Bifurcation Theory; Springer-Verlag: New York, 1985; Vol. 1. Grossmann, I. E.; Morari, M. Operability, Resiliency, and Flexibility-Process Design Objectives for a Changing World. In Proc. 2nd Int. Conf. on Found. Camp.-Aided Proc. Des.; Westerberg, A. W., Chien, H. H. Y., Eds.; CACHE: New York, 1984; p 937. Gurel, 0.;Gurel, D. Oscillations in Chemical Reactions; SpringerVerlag: New York, 1983. Hindmarsh, A. C. LSODE and LSODI, Two New Initial Value Ordinary Differential Equation Solvers. ACM-Signum N e d . 1980, 15 (4), 10. Holodniok, M.; Kubicek, M. DERPER-An Algorithm for the Continuation of Periodic Solutions in Ordinary Differential Equations. J . Comp. Phys. 1984, 55, 254. Jensen, K. F.; Ray, W. H. The Bifurcation Behavior of Tubular Reactors. Chem. Eng. Sci. 1982, 37 (2), 199-222. Kovach, J. W., 111; Seider, W. D. Heterogeneous Azeotropic Distillation: Experimental and Simulation Results. AIChE J . 1987a, 33 (8), 1300-1314. Kovach, J. W., 111; Seider, W. D. Heterogeneous Azeotropic Distillation: Homotopy-continuation Methods. Comput. Chem. Eng. 1987b, 11 (6), 593-605. Kubicek, M.; Marek, M. Computational Methods in Bifurcation Theory and Dissipative Structures; Springer-Verlag: New York, 1983.

--

818

Ind. Eng. Chem. Res., Vol. 29, No. 5 , 1990

Kunkel, P. Quadratically Convergent Methods for the Computation of Unfolded Singularities. SIAM J . Numer. Anal. 1988, 25 (6), 1392. Li, W. C.; Biegler, L. T. Process Control Strategies for Constrained Nonlinear Systems. Ind. Eng. Chem. Res. 1988, 27, 1421. Lin, W.-J.; Seader, J . D.; Wayburn, T . L. Computing Multiple Solutions to Systems of Interlinked Separation Columns. AIChE J . 1987, 33, 886. Mangasarian, 0. L. Equivalence of the Complementarity Problem to a System of Nonlinear Equations. SIAM J . Appl. Math. 1976, 31 ( l ) ,89. Matsuura, T.; Kato, M. Concentration Stability of the Isothermal Reactor. Chem. Eng. Sci. 1967,22, 171-184. Morshedi, A. M. Universal Dynamic Matrix Control. In Chem. Proc. Cont.-CPC I I I , Morari, M., McAvoy. T. J., Eds.; CACHE, Elsevier: New York, 1986; p 547. Murtagh, B. A,; Saunders, M. A. MINOS 5.0 User’s Guide; Systems Opt. Lab, Stanford University: Stanford, CA, 1983. Palazoglu, A,; Arkun, Y. A Multiobjective Approach to Design Chemical Plants with Robust Dynamic Operability Characteristics. Comput. Chem. Eng. 1986, 10 (61, 567-575. Palazoglu, A,; Arkun, Y. Design of Chemical Plants with Multiregime Capabilities and Robust Dynamic Operability Characteristics. Compuc. Chem. Eng. 1987, 1 1 (3), 205-216. Pantelides, C. C. SPEEEDUP-Recent Advances in Process Simulation. Comput. Chem. Eng. 1988, 12 (7), 745. Patwardhan, A. A.; Rawlings, J. B.; Edgar, T. F. Model Predictive Control of Nonlinear Processes in the Presence of Constraints. AIChE Annual Meeting, Washington, D.C., 1988. Penlidis, A,; MacGregor, J. F.; Hamielec, A. E. Continuous Emulsion Polymerization: Design and Control of CSTR Trains. Chem. Eng. Sci. 1989, 44 (2), 273-281. Petzold, L. R. Differential/Algebraic Equations are not ODES. SIAM J . Sci. Stat. Comput. 1982a, 3 (3), 367. Petzold, L. R. A Description of DASSL: A Differential/Algebraic System Solver. Sandia Report SAND82-8637, 1982b. Rawlings, J. B.; Ray, W. H. Stability of Continuous Emulsion Polymerization Reactors: A Detailed Model Analysis. Chem. Eng. S C ~1987, . 42 i l l ) ,2767-2777. Razon, L. F.; Schmitz, R. A. Multiplicities and Instabilities in Chemically Reacting Systems-A Review. Chem. Eng. Sci. 1987. 42 (5), 1005-1047. Renfro, J. G.; Morshedi, A. M.; Asbjornsen, 0. A. Simultaneous Optimization and Solution of Systems Described by Differential/Algebraic Equations. Comput. Chem. Eng. 1987, 11, 503. Rheinboldt, W.C.; Burkardt, J. V. A Locally Parameterized Continuation Process. ACM Trans. Math. Software 1983, 9 (21, 215-235. Richalet, J. A.; Rault, A.; Testud, J. D.; Papon, J. Model Predictive Heuristic Control: Applications to Industrial Processes. Automntica 1978, 14, 413.

Schork, F. J. The Dynamics of Continuous Emulsion Polymerization Reactors. Ph.D. Thesis, University of Wisconsin, Madison, 1981. Seader, J. D.; Kuno, M.; Lin, W.-J.; Johnson, S. A,; Unsworth, K.: Wiskin, J. Mapped Continuation Methods for Computing All Solutions to General Systems of Nonlinear Equations. Comput. Chem. Eng. 1990, I 4 (l),71-85. Seider, W. D.; Brengel, D. D.; Widagdo, S. Nonlinear Analysis in Process Design-A Review. Submitted to AIChE J . 1990. Seider, W. D.; Ungar, L. H. A Course in Nonlinear Systems. Chem. Eng. Educ. 1987, Fall, 178-183. Stupin, W. J.; Lockhart, F. J. Thermally Coupled Distillation-A Case History. Chem. Eng. Prog. 1972, 68 (lo), 71. Swaney, R. E.; Grossmann, I. E. An Index for Operational Flexibility in Chemical Process Design. Part I: Formulaticn and Theory. AIChE J . 1985a, 31 (4), 621-630. Swaney, R. E.; Grossmann, I. E. An Index for Operational Flexibility in Chemical Process Design. Part 11: Computational Algorithms. AIChE J . 1985b, 31 (4), 631-641. Tedder, D. W.; Rudd, D. F. Parametric Studies in Industrial Distillation. AIChE J . 1978, 24, 303. Uppal, A.; Ray, W. H.; Poore, A. B. On the Dynamic Behavior of Continuous Stirred Tank Reactors. Chem. Eng. Sci. 1974, 29, 967-985. Vasantharajan, S.; Biegler, L. T. Simultaneous Solution of Reactor Models within Flowsheet Optimization. Chem. Eng. Res. Des. 1988, 66, 396-407. Vasantharajan, S.; Logsdon, J.; Biegler, L. T . Simultaneous Strategies for Optimization of Differential-Algebraic Systems. 1988 Annual AIChE Meeting, Washington, D.C., 1988. Vasudevan, G.; Watson, L. T.; Lutze, F. H. A Homotopy Approach for Solving Constrained Optimization Problems. In Proceedings of 1989 Amer. Cont. Conf., Pittsburgh; American Automatic Control Council: Green Valley, AZ, 1989; pp 780-785. Villadsen, J.; Michelsen, M. L. Solution of Differential Equation Models by Polynomial Approximation; Prentice-Hall: Englewood Cliffs, NJ, 1978. Watson, L. T.; Billups, S. C.; Morgan, A. P. Algorithm 652, HOMPACK: A Suite of Codes for Globally Convergent Homotopy Algorithms. ACM Trans. Math. Software 1987, 13, 281. Westerink, E. J . Heat Effects in Chemical Reactors: Design and Stability. Doctoral Thesis, Universiteit Twente, The Netherlands, 1987. Widagdo, S.; Seider, W. D.; Sebastian, D. H. Bifurcation Analysis in Heterogeneous Azeotropic Distillation. AIChE J . 1989, 35 (9), 1457-1464. Wohlfahrt, K.; Emig, G. Compare Maleic Anhydride Routes. Hvdrocarbon Process. 1980, 59, 83. Receiued for review October 24, 1989 Accepted December 8, 1989