Verification of Controllers in the Presence of Uncertainty: Application

Jul 3, 1996 - A process design framework for considering the stability of steady state operating points and Hopf singularity points in chemical proces...
0 downloads 9 Views 287KB Size
Download PDF
Ind. Eng. Chem. Res. 1996, 35, 2277-2287

2277

Verification of Controllers in the Presence of Uncertainty: Application to Styrene Polymerization Evangelia Gazi,† Warren D. Seider,* and Lyle H. Ungar‡ Department of Chemical Engineering, University of Pennsylvania, Philadelphia, Pennsylvania 19104-6393

This paper describes a methodology for verifying the stability and performance of closed-loop nonlinear systems, as arise in chemical plants with process controllers, when the process model is incompletely known, with uncertainty both in the parameter values and in the form of the equations. A systematic way of generating and summarizing all possible behavior patterns of the system is described. Temporal logic is used to formalize the posing of and automatic answering of qualitative questions about the response of the system. The methodology is demonstrated for controller verification in the design of a CSTR for styrene polymerization. Introduction In recent years, increased competition has motivated the design of high-performance processes that operate in regions of complex dynamics. To achieve more profitable designs that satisfy strict safety regulations, controllers are increasingly being designed during process synthesis. However, in the design stage, experimental data are often scarce and the models exhibit substantial uncertainties. For this reason, it is of practical importance to develop tools that assist the engineer in designing economical and safe plants and verifying that the control system will respond adequately to anticipated disturbances in the presence of uncertainty. Because the state-variable trajectories produced by a numerical simulation may not be indicative of either typical or worst-case behavior, a methodology for the automatic generation of all of the possible closed-loop behavior patterns, providing quality assessment of the performance and stability of a potential closed-loop system when the process model is incompletely known, is presented herein. With this off-line control-quality assessment, the engineer can locate possible design flaws, guarantee that the process will not move into hazardous operating regimes when disturbances or faults are encountered, and avoid the large overdesign factors that are often used to compensate for model uncertainties. In this paper, the methodology for verifying the performance and stability of the closed-loop system when the process model is incompletely known is demonstrated for the polymerization of styrene in a CSTR. Since the material properties of polymers are largely determined by the kinetic mechanism and the operation of the polymerization reactor, a good understanding and close control of these processes is essential to the production of high-quality polymers. Thus far, advances in control theory have provided new tools to attack the polymerization control problem using differential equation models of the process with uncertain parameters (having parametric uncertainty). In polymerization models, however, the kinetic mechanisms are not well understood, and, hence, the form of the equations is uncertain as well (having nonparametric uncertainty). * Author to whom correspondence should be addressed. Fax: (215) 573-2093. E-mail: [email protected]. † Present address: Chemical Industry Institute of Toxicology, 6 Davis Dr., Research Triangle Park, NC 27709. E-mail: [email protected]. ‡ E-mail: [email protected].

S0888-5885(95)00436-2 CCC: $12.00

Although techniques exist to deal with parametric uncertainty, the representation and simulation of nonlinear models with inexact functional relationships have not been addressed as effectively in the literature. Furthermore, the existing verification techniques (Moon et al., 1992; Alur et al., 1993; Henzinger and Ho, 1995) apply only to discrete-event or discrete-continuous systems but cannot be used for the verification of continuous nonlinear systems. The methodology proposed herein is applied for the control-quality assessment of a styrene polymerization reactor. Although the control of this reactor has been studied extensively (Congalidis et al., 1989; Hidalgo and Brosilow, 1990; Soroush and Kravaris, 1993), results for a limited number of simulations for specific parameter values have been reported; that is, the controllers have not been verified in the presence of uncertainty. The remainder of the paper is organized as follows. In the next section, two techniques are reviewed for the solution of semiquantitative models to predict all of the possible behavior patterns. One involves qualitative analysis, implemented in the NSIM program, that uses symbolic manipulation to build bounding equations, and the other is a nonparametric Monte Carlo technique. Then, temporal-logic operators are introduced to formalize the posing and automatic answering of questions of the type, “Can the system operate adequately if a disturbance/fault occurs?” Finally, the algorithm is applied for control-quality assessment in the design of a CSTR for styrene polymerization. The process model is presented, and a control structure is suggested and verified to perform well in response to anticipated disturbances in spite of the model uncertainties. Analysis of Inaccurate Models Several qualitatively distinct behavior patterns may be consistent with a semiquantitative model, one in which the parameters and terms are uncertain. For example, an increase in the reactor temperature may result either in a desirable increase in the reaction rates or in undesirable instability, depending on the values of the parameters. Hence, it is important to know the full range of behavior patterns consistent with an uncertain model. Thus far, there have been two approaches to the development of tools for solving semiquantitative models to predict all of the possible behavior patterns. The first involves a “worst-case”, deterministic analysis, in which ordinary differential equations (ODEs) are derived whose solution is guaranteed to bound all of the solutions of the inexact model, © 1996 American Chemical Society

2278

Ind. Eng. Chem. Res., Vol. 35, No. 7, 1996

Figure 1. Bounding envelopes for f{x}.

implemented in the NSIM algorithm. The second is a statistical analysis that extends the Monte Carlo technique to systems with nonparametric uncertainty. The two techniques are described briefly below. It is assumed that the plant model contains unknown functions, y ) f{x}, which are bounded between envelopes, such that fl{x} e y e fu{x}, and are monotonic (Figure 1). The monotonicity assumption is very important for the analysis but is not overly restrictive, since any nonmonotonic function can be constructed by combining monotonic ones. In addition, it is often known that physical quantities are related to each other monotonically, even though the relations may not be known accurately; e.g., reaction rates increase monotonically with temperature, and friction factors increase monotonically with velocity. Kay and Ungar (1993) present a systematic method for determining the bounding envelopes for multivariable monotonic functions, given a data set. NSIM Algorithm. The QSIM algorithm (Qualitative SIMulator; Kuipers, 1986) provides a framework for developing models in the form of qualitative differential equations (QDEs), which are abstractions of ODEs, where the values of the variables are described qualitatively and the functional relationships between the variables may be known incompletely. The constraints are qualitative expressions involving the common mathematical operations, such as addition, multiplication, and differentiation, and monotonic functions. Given a QDE and a qualitative description of an initial state, QSIM derives a tree of qualitative state descriptions, where the paths from the root to the leaves of the tree represent the possible behavior patterns of the system. In NSIM (Numerical SIMulator; Kuipers, 1989; Kay and Kuipers, 1993), an extension of QSIM, parameters may be assigned lower and upper bounds and monotonic functions may be bounded by lower and upper envelopes, thus forming semiquantitative differential equations (SQDEs). These permit NSIM to eliminate behavior patterns that are inconsistent with the numerical information, while describing the remaining behavior patterns more quantitatively. NSIM produces bounds on the trajectories of the state variables that are guaranteed to include all of the behavior patterns consistent with the approximate model. Given a SQDE and an initial condition, NSIM derives and numerically integrates an “extremal” system of ODEs, whose solution is guaranteed to bound all of the solutions of the SQDE. This is referred to as NSIM simulation. The extremal system is generated automatically, using symbolic manipulation, based on a set of rules (Kay and Kuipers, 1993). Minimal and maximal equations, which bound the first derivatives of the state variables, are derived, defining a dynamic envelope for each variable, and, hence, an nth-order system is transformed to a 2nth-order system. As a result, the extremal ODEs are not generally members of the class

of ODEs represented by the SQDE. Therefore, the dynamic envelopes do not necessarily have the same shape as the behavior patterns of the SQDE. The extremal system is integrated using a standard RungeKutta technique. It is proven theoretically that, when the extremal system bounds the solution of the SQDE at t ) 0, it bounds the solution at all times (Kay, 1991). A SQDE is an abstraction of an ODE with incompletely specified functions and parameter values. Stated differently, one NSIM simulation bounds the solutions of an infinite number of ODE models. Since the NSIM bounds are guaranteed to bound all of the solutions to the SQDE, NSIM can be used for automatic proof construction. However, the main disadvantage of this technique is that the bounds are often too conservative, as illustrated by Gazi et al. (1994). Nonparametric Monte Carlo Technique. To extend the Monte Carlo technique to handle nonparametric uncertainty, a large number of monotonic functions are randomly generated within the envelopes that bound the incompletely known monotonic functions. Each generated function is incorporated into the ODEs that model the process, which, when integrated, produce trajectories of the state variables that are systematically checked for interesting behavior patterns, as described in the next section. This is accomplished using the following algorithm: (1) M intervals in x are formed. (2) At x ) x1, y1 is selected randomly such that fl{x1} e y1 e fu{x1}. At x ) xk, yk is selected randomly such that fl{xk} e yk e fu{xk} and yk g yk-1, k ) 2, ..., M. (3) The points (xk, yk), k ) 1, ..., M, are connected with straight lines. (4) The system of ODEs is integrated. (5) Steps 2-4 are repeated. In this analysis, it is important to ensure a satisfactory coverage of the space of possible functions. Gazi et al. (1996) show that the monotonicity constraint clusters the generated functions close to the upper envelope, as x increases, when the y’s are selected from uniform distributions. To avoid this, y1 and yM are selected randomly first, and then the x coordinate is subdivided recursively. A large number of simulations (one with each generated function) is required, and the results are guaranteed to include all possible behavior patterns only in the limit as the number of samples approaches infinity. Controller Verification. During process and controller design, because the specific trajectories produced by numerical simulations may be unreliable, qualitative descriptions of the possible closed-loop behavior patterns enable the engineer to establish bounds on a desired nominal region of operation and automatically verify offline whether, for example, the controller can return the process to a steady state within the nominal region in response to anticipated disturbances and faults. Given the nominal, semiquantitative, model for the process and the controller, as well as the potential faults, the proposed verification scheme is as follows. The engineer asks a question of the type, “What may happen if this fault occurs?” In response, the fault model is automatically built and integrated (using either of the techniques described above). Then, temporal-logic operators (Emerson, 1990) are used to express formally the usersupplied questions about the system behavior and to provide qualitative answers. A similar scheme has been applied for the verification of discrete-event control systems by several authors. For

Ind. Eng. Chem. Res., Vol. 35, No. 7, 1996 2279 Table 1. Temporal-Logic Operators modal

temporal

logical

necessarily possibly

always eventually until

and or not

example, in Moon et al. (1992), the model that describes the process and its control software is in the form of a state transition graph, where the state variables may only have true or false values. Computational tree logic (CTL; Emerson, 1990) formulas are used to answer questions about the operability, reliability, and safety of the system. Given a CTL formula, the state transition graph is unwound into a tree with the initial state at the root, to be checked systematically to determine the truth of the question. Unfortunately, this methodology cannot be extended to dynamic systems, for which a transition-graph model cannot be generated. Kuipers and Shults (1994) apply CTL*, an extension of CTL, to QSIM trees. As mentioned earlier, the QSIM algorithm simulates QDEs to generate a tree of qualitative states representing a branching-time description of the possible behavior patterns of the system. This permits temporal logic to be applied in much the same way as for discrete-event systems, but with qualitative or quantitative variables and derivatives, representing a much broader class of systems. The application of this methodology is also limited, however, due to the intractable branching of the QSIM trees for complicated dynamic systems. While CTL* is applied by Kuipers and A° stro¨m (1994) to verify the heterogeneous control of a CSTR with a reversible exothermic reaction, the assumption that the net rate of reaction is zero (quasisteady-state operation) was necessary to obtain a tractable tree. In a parallel contribution, Alur et al. (1993) introduce a model-checking procedure for the automatic verification of hybrid systems (having discrete and continuous variables) that are linear or linearly bounded. Henzinger and Ho (1995) extend this approach to nonlinear systems by discretizing the continuous variables and bounding the nonlinear functions in each region with piecewise linear envelopes. However, this technique is illustrated for the control of a nonlinear railroad gate with only one continuous state variable. Its extension, although theoretically straightforward, is likely to yield conservative bounds similar to NSIM (Gazi et al., 1994). Herein, temporal-logic operators are used to pose questions about the Monte Carlo simulations. The trajectory of the state variables produced from each Monte Carlo run may be viewed as one behavior, with each time point representing a state. Since there is no branching, a simplified CTL is utilized, with the temporal-logic operators shown in Table 1. Note that Gazi et al. (1996) provide a methodology for representing the Monte Carlo simulations in “QSIM-like” trees and generating qualitative descriptions of the possible behavior patterns. Through the modal operators, the user can specify that the answer to the question must be true for all Monte Carlo simulations (“necessarily”) or at least one (“possibly”). Note that the verification algorithm may be applied in parallel with the Monte Carlo simulations; i.e., checked as each trajectory is generated and classified as true or false. When the modal operator is “possibly”, the simulations may be terminated after one behavior is found, with an affirmative response returned. Alternatively, when the modal operator is “necessarily” and a false behavior is found, providing a

Figure 2. Schematic of a jacketed CSTR for styrene polymerization.

counterexample, a negative response is returned. Thus, unnecessary simulations are avoided and the results of all Monte Carlo runs need not be stored. When using NSIM, the computed bounds are checked using the temporal and logic operators, but the modal operator is set to “necessarily”. Questions often involve one of three temporal operators, referring to the full time range (“always”), or a range of time (“eventually”, “until”), as well as the logical operators, which are the conventional Boolean connectives. For an example question, expressed using the temporal-logic operators, see the next section. Polystyrene Reactor Process ModelsDesign and Operating Conditions. The model for the free-radical polymerization of styrene in a jacketed CSTR is presented below. As shown in Figure 2, the CSTR has three feed streams for pure styrene monomer, azobis(isobutyronitrile) (AIBN) initiator dissolved in benzene, and benzene solvent. Heat generated by the exothermic reactions is removed by cooling water in the jacket. The exit stream contains polymer, unreacted monomer, initiator, and solvent. In the free-radical polymerization, the following elementary reactions are postulated (Jaisinghani and Ray, 1977): kd

Initiation

I 98 2R ki

M + R 98 P1 Propagation

kp

Px + M 98 Px+1

Termination ktd

Px + Py 98 Tx + Ty (disproportionation) ktc

Px + Py 98 Tx+y (combination) The two initiation reactions involve the decomposition of initiator to produce radicals, R, which react with the monomer molecules to initiate new live (radical) polymer chains, P1. During the propagation step, monomer molecules are added, one at a time, to the live-polymer chains. The growth of the chains terminates when the propagating radicals lose their activity through either of the termination reactions, resulting in dead-polymer chains, T. The active-polymer chains may collide with, and transfer radicals to, monomer or solvent molecules.

2280

Ind. Eng. Chem. Res., Vol. 35, No. 7, 1996

However, for this reaction system, these collisions may be neglected. To obtain the full chain-length distribution of the polymer, a very large number of mass balances (ODEs) must be solved. In practice, however, it is sufficient to calculate only the leading moments of the chain-length distribution, where the kth moment of the live-radical and dead-polymer chain-length distributions are defined, respectively, as: ∞

λk )

∞

∑xkPx x)1

µk )

∑xkTx x)1

(1)

Since the radical concentrations typically range between 10-4 and 10-8 mol/L and reach steady values in less than 1 s (Penlidis et al., 1992), the quasi-steadystate approximation can be applied for λk, k ) 0, 1, 2, .... Then, the following mass balances for the initiator (2) and monomer (3), the overall mass balance (4), the energy balances for the reactor holdup (5) and jacket (6), and the balances that define the moments of the dead-polymer distribution (7-9), as described by Penlides et al. (1992), coupled with the lower and upper bounds on the expressions for kp and kt (see the section on Sources of Uncertainty), provide a model to design the reactor and to project its controllability in the design stages:

dI FiIf - FtI ) - k dI dt V

(2) (3)

dV ) Fin - Ft ) Fi + Fm + Fs - Ft dt

(4) (5)

dTc Fc(Tcf - Tc) hA ) + (T - Tc) dt Vc FccpcVc

(6)

Ftµ0 ktcλ02 dµ0 )+ + ktdλ02 dt V 2

(7)

dµ1 Ftµ1 )+ ktλ0λ1 dt V

(8)

dµ2 Ftµ2 )+ ktcλ12 + ktλ0λ2 dt V

(9)

where

P ≡ λ0 )

[ ] 2fkdI kt

1/2

(10)

is the concentration of the live-polymer chains of any length, and

λ1 )

2fkdI + kpMλ0 ktλ0

λ2 ) λ1 +

2kpMλ1 ktλ0

The initiator efficiency, f, is the fraction of initiator radicals that undergoes the addition of any number of monomer units. For AIBN, it is estimated to lie within the interval 0.6-0.7. The most important quantity that characterizes the molecular-weight distribution (MWD) is the weightaverage molecular weight, Mw, defined as: ∞

Mw )

dM FmMf - FtM ) - kpMP dt V

hA dT FinTf - FtT -∆Hr ) + (T - Tc) k MP dt V Fcp p FcpV

Figure 3. Heat generation and removal in the CSTR.

(11)

(12)

∑ i)1

∞

wiMi )

∑ i)1

∞

NiMi2/

NiMi ∑ i)1

(13)

where Ni is the number, and wi is the weight fraction of polymer chains with molecular weight Mi, respectively. Using the definition of λi and µi, it can be shown that:

µ2 + λ2 Mw ) Wm µ1 + λ1

(14)

where Wm is the molecular weight of the monomer. To assure product quality and control of these highly exothermic reactors, two important design and operating constraints are maintained: (1) the volume fraction of solvent is kept at 60%, to avoid a significant increase in the gel effect (Choi, 1986), and (2) the reactor temperature is bounded by 150 °C (423 K) since the rate of thermal initiation becomes significant at temperatures higher than 100 °C (373 K) and probably dominates catalytic initiation at temperatures higher than 150 °C (Biesenberger and Sebastian, 1983). Figure 3 shows a plot of the heat generated and the heat removed, under steady-state conditions, as a function of temperature and illustrates the three steady states, for the parameters, thermophysical properties, and the design and operating variables in Tables 2 and 3. Note that similar behavior is exhibited for a wide range of conditions. Because the restrictions above exclude operation at the high-conversion, stable steady state and the conversion at the low-conversion, stable steady state is too low, operation at the middle-conversion, unstable steady state is required. To achieve reliable operation at an unstable steady state, in the face of typical disturbances, a reliable control system must be designed, as discussed in the next section. In this regard, it is noted that the µk, k ) 0, 1, 2, do not appear in eqs 2-6, and, hence, the reactor model is comprised of two subsystems. When a measurement of the molecular-weight distribution is not available, eqs 2-6 relate the measured temperatures

Ind. Eng. Chem. Res., Vol. 35, No. 7, 1996 2281 Table 2. Nominal Values of Parameters, Thermophysical Properties, and Design and Operating Variables (Hidalgo and Brosilow, 1990)a f Ad Ed At Et Ap Ep -∆Hr hA Fcp Fccpc a

0.6 5.95 × 1013 s-1 14 897 K 1.25 × 109 L mol-1 s-1 843 K 1.06 × 107 L mol-1 s-1 3557 K 16 700 kcal mol-1 70 kcal K-1 s-1 360 kcal L-1 s-1 966.3 kcal L-1 s-1

kd ) Ade-Ed/T, kt ) Ate-Et/T, and kp ) Ape-Ep/T.

Table 3. Design and Operating Variables at Steady State flowrates: initiator ) 0.03 L s-1 solvent ) 0.1275 L s-1 monomer ) 0.105 L s-1 cooling water ) 0.131 L s-1 volumes: reactor ) 3000 L jacket ) 3312.4 L feed concentrations: initiator ) 0.5888 mol L-1 monomer ) 8.6981 mol L-1 feed temperatures: reactor ) 330 K cooling water ) 295 K steady-state values at unstable steady state: initiator concentration ) 0.0486 mol L-1 monomer concentration ) 2.337 mol L-1 reactor temperature ) 343.1 K jacket temperature ) 316.2 K

to the flow rates of the feed streams and cooling water, and the molecular weight is controlled implicitly. With measurements to characterize the molecular-weight distribution, the second subsystem (eqs 7-14) also applies as it relates the weight-average molecular weight, Mw, to the input flow rates, and a control loop that directly adjusts Mw may be implemented. Control of the Polymerization Reactor. The control of polymerization reactors is challenging because of: (a) the inability of on-line sensors to provide fast sampling rates with small time delays, (b) the complex nonlinear behavior of polymerization reactors, and (c) the extreme sensitivity of the steady states to small changes in the operating conditions (Ogunnaike and Ray, 1994). Due to increased competition and emphasis on product quality, as well as the need to operate at unstable steady states, the methods of recipe control, in which the operating conditions are updated periodically, given product analyses from the laboratory, have become inadequate. The lack of reliable and easily-implemented on-line measurements has motivated several studies of stateestimation techniques to estimate unmeasured polymer properties. For example, Jo and Bankoff (1976) introduced an extended Kalman filter to estimate the conversion and weight-average molecular weight in the solution polymerization of vinyl acetate using temperature measurements in a glass CSTR. Recently, studies were carried out to determine relationships, both qualitative and quantitative, among measurements of the density, viscosity, refractive index, and similar properties and the desired polymer properties (McGreavy, 1994). To improve the control algorithms, three recent studies are particularly noteworthy. The first, by Hidalgo and Brosilow (1990), combines model-predictive and

Figure 4. Proportional gain as a function of reactor temperature.

coordinated control strategies for the solution polymerization of styrene in a jacketed CSTR at an unstable steady state but does not attempt to control the molecular weight of the polymer. In a second study, Congalidis et al. (1989) implement a feedforward-feedback controller, using PI controllers, to operate a solution copolymerization reactor. However, frequent on-line measurements of Mw are assumed to be available and the reactor is operated at a low-conversion, stable steady state. Finally, in a comprehensive study, Soroush and Kravaris (1993) introduce empirical equations to infer conversion from on-line measurements of the density and temperature and implement so-called globallylinearizing control (GLC) for the solution polymerization of methyl methacrylate. The first objective of the controller implemented herein is to stabilize the reactor using four possible control efforts: Fi, Fm, Fs, and Fc. Since it is desired to maintain a nearly-constant volume fraction of solvent in the reactor, the ratio of the solvent/monomer flow rates, Fs/Fm, is controlled. Also, since the initiator concentration primarily influences the polymer chain length (Congalidis et al., 1989), the flow rate of the initiator, Fi, controls the molecular-weight distribution. The remaining two control actions, the cooling-water and monomer flow rates, are used to control the reactor temperature, as suggested by Hidalgo and Brosilow (1990). The cooling-water flow rate is the primary manipulated variable, with the monomer flow rate used as a secondary control action when the cooling-water valve approaches saturation. Due to the nonlinearities of the reactor, especially in operation at an unstable steady state, some form of nonlinear control is desirable. However, most nonlinear controllers require a precise process model, in some cases with provision for parameter updates to reduce process/model mismatch (e.g., Brengel and Seider, 1989, 1992). To simplify the temperature controller, especially for the projections of controllability needed during process design, sliding-gain scheduling of a PID controller (Shamma and Athans, 1990) is implemented herein. Like a nonlinear controller, the controller gain varies with temperature, as shown in Figure 4, but the controller is much simpler in that it does not utilize the process model at each sampling instant. The controller is tuned by trial and error, since the saturation bounds on the manipulated variables do not permit the use of IMC tuning rules (Zafiriou and Morari, 1989). Furthermore, the following rules are adhered to:

2282

Ind. Eng. Chem. Res., Vol. 35, No. 7, 1996

Table 4. Tuning Parameters controlled manipulated variable variable

Kc

conditions

τI

τD

106

103

1.4 × 104

102 Fc saturated, or |T - Tsp| g 2 K, or T > 380 K 4 not in the vicinity of Mw,sp

T

Fc

T

Fm

given by Figure 4 0.3 105

Mw

Fi

10-6

(1) 0 e Fc e 1.0 L/s, and 0 e Fm e 1.0 L/s. (2) When the temperature is close to its setpoint, little control action is applied, to avoid driving the reactor away from the unstable steady state. (3) As the temperature deviates further from its setpoint, the control action is strengthened, to avoid switching to the undesirable, stable steady states, which are competing attractors. (4) Fm is not manipulated unless necessary, since it also affects the molecular-weight distribution as well as the production rate (Congalidis et al., 1989). (5) When |T - Tsp| g 2 K, it is anticipated that Fc will become saturated before the disturbance is rejected, and, hence, Fm is also manipulated. (6) When the temperature exceeds 380 K, both Fc and Fm are manipulated, to keep the reactor below the region of thermal initiation. (7) Large changes in Fc, that create sharp temperature changes, are avoided by requiring that |∆Fc| e 0.2 L/s. Table 4 and Figure 4 provide the parameters for the temperature controller. As mentioned earlier, on-line measurements of Mw are rarely available. Instead, the intrinsic viscosity (Braun et al., 1971), η, is measured and Mw is inferred using a power-law expression, η ) aMwb, where a and b are determined by regression of experimental data. A formal definition of the intrinsic viscosity is given in Appendix A. The flow rate of the initiator, Fi, is used to control the molecular weight. Note that, as the concentration of the initiator increases, the rate of formation of active-polymer chains increases, and, therefore, for a constant propagation rate, the average molecular weight decreases. Since the relationship between Fi and Mw does not exhibit strong nonlinearities, PID control of Mw is adequate, with the controller parameters shown in Table 4. Note that, because the manipulation of Fi implicitly couples the two subsystems, only small changes in Fi are allowed (|∆Fi| e 0.002 L/s), to avoid large departures in the temperature and monomer conversion. In addition, a dead-band, involving no control action, is introduced near the Mw setpoint, due to the uncertainties in Mw, as discussed in the Results section. Sources of Uncertainty. Although styrene is one of the most widely-studied monomers, there remains considerable uncertainty in the model, especially in the rate constants of propagation, kp, and termination, kt. This is illustrated by the large range of values reported in the Polymer Handbook (Brandrup and Immergut, 1989). The overall chain-termination rate constant, kt, is the sum of the combination and disproportionation contributions, kt ) ktc + ktd, with each assumed to have an Arrhenius dependence on temperature (Schmidt and Ray, 1981). Furthermore, it is assumed often that ktc = ktd, which may be a poor approximation for some systems. Since their relative contributions to kt are not estimated easily and may vary with temperature and

Figure 5. Gel effect correlations for the bulk polymerization of methyl methacrylate at 70 °C [(1) Cardenas and O’Driscoll, 1976; (2) Friis and Hamielec, 1976; (3) Ross and Laurence, 1976; (4) Marten and Hamielec, 1979]. Reprinted with permission from Schmidt and Ray (1981).

composition, there is nonparametric uncertainty in the overall termination rate constant, kt. This is accounted for herein by defining the following envelopes:

ktl ) 1.0 × 109e-1000/T

ktu ) 2.6 × 109e-820/T (15)

Note that although the bounding envelopes have identical functional forms, with different At and Et, the methods for nonparametric uncertainty are utilized to compute the results presented below. Since ktc and ktd are not measurable, they are estimated as a fraction of kt, with ktc ) νkt and ν in the range 0.4-0.6. In addition, as the monomer conversion increases, the termination rate constant has the special form kt ) kt0gt, where kt0 is the rate constant at very low monomer conversions, while gt, a function of temperature and composition, falls due to strong diffusion limitations at higher conversions. This phenomenon is known as the gel or Trommsdorff effect. Although many models successfully reproduce operating data, there is no universally-accepted relationship between gt and the molecular weight, polymer concentration, and temperature, all of which determine the gel effect (Biesenberger and Sebastian, 1983). Several common correlations for gt, as a function of the monomer conversion, x, are shown in Figure 5 for the bulk polymerization of methyl methacrylate (MMA). Nearly identical correlations exist in the literature for styrene (see Schmidt and Ray, 1981, for a thorough review), all of which are bounded by the following envelopes:

gtl ) e-8x,

gtu )

{

1 e-2.85(x-0.3)

0 e x e 0.3 (16) 0.3 e x e 1

Note that the kinetics are very sensitive to these parameters. In this regard, Knorr and O’Driscoll (1970) were the first to show that the gel effect can explain steady-state multiplicity under isothermal conditions for styrene polymerization. Uncertainties in the propagation rate constant are important, not only because it is proportional to the rate of polymer conversion but because it strongly affects the MWD. Despite its significance, kp has proven to be a difficult parameter to determine quantitatively. Herein, nonparametric uncertainty is accounted for using the envelopes:

kpl ) 9.0 × 106e-3700/T

kpu ) 2.2 × 107e-3500/T (17)

Ind. Eng. Chem. Res., Vol. 35, No. 7, 1996 2283

Figure 6. P&ID for styrene polymerization in a CSTR.

Figure 7. NSIM bounds on the reactor temperature for a setpoint increase by 5-8 K.

Since kp is not strongly affected by the gel effect, it is assumed that gp = 1, except for extremely viscous solutions close to the glass transition temperature. In this regard, the results presented below are at normal temperatures and viscosities, and, hence, gp ) 1. A major source of uncertainty is associated with the relationship between the intrinsic viscosity and the weight-average molecular weight, Mw. Since the parameters of the power-law expression, η ) aMwb, are adjusted to represent scattered experimental data, there is considerable uncertainty in the form of the equation. Furthermore, Mw is a large number, increasing the sensitivity of the correlation to uncertainty. Even worse, there can be significant changes in the MWD, as well as η, while Mw remains constant; that is, the relationships among the MWD, η, and Mw are not well understood. Because of this, the following envelopes are used:

ηl ) 0.0012Mw0.71

ηu ) 0.0013Mw0.73

(18)

to project the performance of a control loop in which the intrinsic viscosity is measured and the flow of initiator is adjusted. Results. As discussed above, the verification technique has been applied to the polymerization reactor, with the controllers shown in Figure 6. First, the control system is tested for the servo problem. The temperature setpoint is increased by 5-8 K, with the setpoint on the weight-average molecular weight, Mw,sp, decreased by 7000-10 000, corresponding to the new steady states. Note that, for each Monte Carlo simulation, the setpoint is set to Tsp ) Tsp, nominal + x, where x

Figure 8. Monte Carlo bounds on the response for a set point increase by 5-8 K.

is a random number chosen from a uniform distribution on the interval 5-8. Figures 7 and 8 show the response of the closed-loop system, as predicted by NSIM and Monte Carlo simulations, respectively, in answer to the question (presented in the QSIM/CTL* syntax) “If the reactor temperature setpoint increases by 5-8 K, is it (necessarily ((always (in nominal region)) and eventually (steady-state)))?” where the “nominal region” of operation is defined such that 340 e T e 360 K, and 20 000

2284

Ind. Eng. Chem. Res., Vol. 35, No. 7, 1996

Figure 9. Monte Carlo bounds on the response for an increase in the feed temperature by 3-5 K.

e Mw e 35 000. Stated differently, if the reactor temperature setpoint increases by 5-8 K, it must be determined whether it is true for all simulations (necessarily) that the system remains in the nominal region for all time intervals during the simulation (always) and that during the last time intervals (eventually) it returns to steady state? The NSIM bounds diverge, showing that unstable trajectories are possible, and, hence, the answer “false” is returned. However, this is attributed to the conservatism of the NSIM bounds. Even when the uncertainties are small, NSIM incorrectly predicts instability, due to the high nonlinearity of the system. To the contrary, the Monte Carlo simulations predict tight bounds for all of the state variables and, hence, provide an affirmative response. For this reason, Monte Carlo simulations are utilized to obtain the remaining results in this section. Next, the controller is verified for disturbance rejection. In response to a 3-5 K increase in Tf, Figure 9 shows that tight bounds are predicted. Similar behavior is predicted for a decrease in the heat-transfer coefficient, h, by 5-10 kcal m-2 K-1 s-1, to account for fouling of the heat-transfer surfaces, as shown in Figure 10. It should be noted that multiple disturbances, that occur either simultaneously or sequentially, may also be simulated and verified. Finally, Figure 11 shows the importance of the MwFi control loop. To examine the effect of this loop, all

Figure 10. Monte Carlo bounds on the response for a decrease in the heat-transfer coefficient by 5-10 kcal m-2 K-1 s-1.

uncertainties are eliminated, except those associated with the η-Mw relationship (18). First, as a reference, Figure 11a shows the response to a 3-5 K increase in Tf when perfect Mw measurements are available. Also, for comparison, Figure 11b shows the response without Mw control. Note that Mw drops by approximately 2000. As expected, the presence of the Mw-η loop does not affect the performance of the temperature controller. With viscosity measurements, Figure 11c shows the bounds on the Mw trajectory. For each Monte Carlo simulation, a random function Mw ) Mw{η} is generated within the envelopes, and, given a measurement, η*, an estimate for Mw, Mw*, is calculated, together with the upper and lower values corresponding to η*, Mw,u* and Mw,l*, respectively. When Mw,l* e Mw,sp e Mw,u*, no control action is applied, because an action based on an imprecise Mw* in this region may drive the reactor in the wrong direction. Outside this region, the control action is based on the error |Mw* - Mw,sp|. Figure 11c shows a significant improvement in the Mw trajectory when viscosity measurements are available, as compared to the response without control in Figure 11b. The technique is well-suited for these comparisons, as well as for the design and quality assessment of controllers when direct measurements are not available.

Ind. Eng. Chem. Res., Vol. 35, No. 7, 1996 2285

Conclusions A methodology has been introduced for the off-line verification of the stability and performance of controllers, when the process model has substantial parametric and nonparametric uncertainty. Two techniques were reviewed for the simulation of such semiquantitative models. Although more expensive computationally than NSIM, with guarantees only in the limit of infinite samples, the Monte Carlo method was shown to produce bounds that are not conservative. The verification algorithm was developed to help the engineer during process design by examining the numerous and detailed Monte Carlo simulations that describe the possible behavior patterns of the system. It checks these behavior patterns against temporal-logic queries and provides qualitative answers to questions about the responses of the system to anticipated disturbances and set point changes. In this way, it has been shown to reject potential unsafe and unstable behavior patterns in the design stages. This algorithm applies to any linear or nonlinear dynamic system, thus greatly extending the applicability of existing verification techniques. It has been demonstrated herein to be effective for control-quality assessment in the design of a CSTR for styrene polymerization. Based on the design specifications, a slidinggain scheduling controller is designed to stabilize the reactor at an unstable steady state and is verified to perform well in response to anticipated disturbances and changes in operating conditions. Furthermore, it is verified that, in the absence of direct measurement of the molecular weight, the control of the molecular weight is improved with viscosity measurements, even though the relationship between molecular weight and viscosity exhibits substantial nonparametric uncertainty. Acknowledgment The authors are grateful to B. J. Kuipers and H. Kay for their assistance with NSIM and to T. Ogunnaike for advice in characterizing the uncertainties in the polymerization of styrene. Helpful discussions with R. IrizarryRivera and M. Soroush are also acknowledged. Partial funding provided by NSF Grant Nos. IRI-9216714, DDM-9114080, and IRI-9216714 is gratefully acknowledged. Nomenclature

Figure 11. Weight-average molecular weight for an increase in the feed temperature by 3-5 K.

To generate the envelopes above, 10 000 Monte Carlo simulations, requiring approximately 180 cpu minutes on a DECstation 5000, were performed. Note that, in all cases, the envelopes remained unchanged after about 1000 simulations.

a ) coefficient in the power-law expression for the intrinsic viscosity, η A ) heat transfer area of CSTR, m2 Ad ) preexponential factor for the initiation rate constant, s-1 Ap ) preexponential factor for the propagation rate constant, L mol-1 s-1 At ) preexponential factor for the termination rate constant, L mol-1 s-1 b ) coefficient in the power-law expression for the intrinsic viscosity, η cp ) mean heat capacity of reactor fluid, kcal g-1 K-1 cpc ) heat capacity of cooling jacket fluid, kcal g-1 K-1 Ed ) activation energy for the initiation rate constant, K Ep ) activation energy for the propagation rate constant, K Et ) activation energy for the termination rate constant, K f ) initiator efficiency; monotonically increasing function Fc ) flow rate of cooling jacket fluid, L s-1

2286

Ind. Eng. Chem. Res., Vol. 35, No. 7, 1996

Fi ) flow rate of initiator stream, L s-1 Fin ) total flow rate into the reactor, L s-1 Fm ) flow rate of monomer stream, L s-1 Fs ) flow rate of solvent stream, L s-1 Ft ) flow rate of exit stream, L s-1 g ) gel effect factor h ) overall heat-transfer coefficient, kcal m-2 K-1 s-1 I ) concentration of initiator, mol L-1 If ) concentration of initiator in feed, mol L-1 kd ) rate constant for initiation of radicals, s-1 ki ) rate constant for initiation of live-polymer chains, L mol-1 s-1 kp ) propagation rate constant, L mol-1 s-1 kt ) overall termination rate constant, L mol-1 s-1 kt0 ) overall termination rate constant at low monomer conversions, L mol-1 s-1 ktc ) rate constant for termination by combination, L mol-1 s-1 ktd ) rate constant for termination by disproportionation, L mol-1 s-1 Kc ) PID controller gain M ) concentration of monomer, mol L-1; number of intervals in x Mf ) concentration of monomer in feed, mol L-1 Mw ) weight-average molecular weight, g mol-1 P ) total concentration of live macroradicals of any length, mol L-1 Px ) concentration of live macroradicals of length x, mol L-1 t ) time, s T ) temperature of reactor, K Tc ) temperature of cooling jacket fluid, K Tcf ) temperature of cooling jacket fluid in feed, K Tf ) temperature of reactor feed, K Tx ) concentration of dead-polymer chains of length x, mol L-1 V ) reactor volume, L Vc ) volume of cooling jacket, L wi ) weight fraction of polymer chains with molecular weight Mi Wm ) monomer molecular weight, g mol-1 x ) conversion of monomer; independent variable y ) dependent variable Greek Letters -∆Hr ) heat of polymerization reaction, kcal mol-1 η ) intrinsic viscosity, L g-1 λi ) ith moment of the live-radical distribution µi ) ith moment of the dead-polymer distribution ν ) ktc/kt F ) mean density of reactor fluid, g L-1 Fc ) density of cooling jacket fluid, g L-1 τD ) derivative time constant in PID controller, s τI ) integral time constant in PID controller, s Subscripts l ) lower bound sp ) set point u ) upper bound

Appendix A For the determination of the molecular weight using the viscosimetric method, only the relative increase in viscosity is important. The viscosity of the solution of the polymer, η, and that of the solvent, η0, are measured, and the specific viscosity is determined:

ηsp )

η - η0 η0

(A.1)

When the measurements are made in capillary viscometers and at low concentrations, the flow times, t of the solution and t0 of the solvent, provide a good approximation:

ηsp )

t - t0 t0

(A.2)

By dividing ηsp by the concentration of the polymer in the solution, the reduced specific viscosity, ηsp/c, is obtained. However, since this quantity depends upon the concentration, the intrinsic viscosity, [η]:

[η] ) lim cf0

ηsp c

(A.3)

is used to characterize the polymer. Infinite dilution measurements are impossible, and, hence, viscometric measurements are carried out at low concentrations and the results are extrapolated to zero concentration. Note that, since ηsp is dimensionless, [η] has the dimension of reciprocal concentration. For simplicity, in eq 18 and elsewhere, the symbol η, instead of [η], is used to represent the intrinsic viscosity. Literature Cited Alur, R.; Henzinger, T. A.; Ho, P.-H. Automatic Symbolic Verification of Embedded Systems. Proc. 14th Annu. IEEE Real-Time Syst. Symp. (RTSS-93) 1993, 2-11. Biesenberger, J. A.; Sebastian, D. H. Principles of Polymerization Engineering; Wiley: New York, 1983. Brandrup, J., Immergut, E. H., Eds. Polymer Handbook; WileyInterscience: New York, 1989. Braun, D.; Cherdron, H.; Kern, W. Techniques of Polymer Syntheses and Characterization; Wiley-Interscience: New York, 1971. 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. Comput. Chem. Eng. 1992, 16 (9), 861-886. Cardenas, J. N.; O’Driscoll, K. F. High-conversion Polymerization. I. Theory and Application to Methyl Methacrylate. J. Polym. Sci. Chem. 1976, 14, 883-897. Choi, K. Y. Analysis of Steady State of Free Radical Solution Polymerization in a Continuous Stirred Tank Reactor. Polym. Eng. Sci. 1986, 26, 975-981. Congalidis, J. P.; Richards, J. R.; Ray, W. H. Feedforward and Feedback Control of a Solution Copolymerization Reactor. AIChE J. 1989, 35 (6), 891-907. Emerson, E. A. Temporal and Modal Logic. In Handbook of Theoretical Computer Science; van Leeuwen, J., Ed.; Elsevier Science Pub. B.V./MIT Press: Cambridge, MA, 1990; pp 9951072. Friis, N.; Hamielec, A. E. Gel-effect in Emulsion Polymerization of Vinyl Monomers. ACS Symp. Ser. 1976, 24, 82-87. Gazi, E.; Seider, W. D.; Ungar, L. H. Controller Verification using Qualitative Reasoning. In Proceedings of 2nd IFAC Workshop on Computer Software Structure Integ. AI/KBS Sys. in Proc. Cont., Lund, Sweden, 1994. Gazi, E.; Seider, W. D.; Ungar, L. H. A Non-parametric MonteCarlo Technique for Controller Verification. Automatica 1996, submitted for publication. Gazi, E.; Seider, W. D.; Ungar, L. H. Automatic Analysis of MonteCarlo Simulations of Dynamic Chemical Plants. Proceedings of ESCAPE-6, May 25-27, Rhodes, Greece; May 1996, in press. Henzinger, T. A.; Ho, P.-H. Algorithmic Analysis of Nonlinear Hybrid Systems. Proc. Conf. Comput.-Aided Verif. 1995, in press. Hidalgo, P. M.; Brosilow, C. B. Nonlinear Model Predictive Control of Styrene Polymerization at Unstable Operating Points. Comput. Chem. Eng. 1990, 14 (4/5), 481-494.

Ind. Eng. Chem. Res., Vol. 35, No. 7, 1996 2287 Jaisinghani, R.; Ray, W. H. On the Dynamic Behavior of a Class of Homogeneous Continuous Stirred Tank Polymerization Reactors. Chem. Eng. Sci. 1977, 32, 811-825. Jo, J. H.; Bankoff, S. G. Digital Monitoring and Estimation of Polymerization Reactors. AIChE J. 1976, 22, 361. Kay, H. Monitoring and Diagnosis of Multi-tank Flows using Qualitative Reasoning. MS Thesis, The University of Texas at Austin, Austin, TX, 1991. Kay, H.; Kuipers, B. J. Numerical Behavior Envelopes for Qualitative Simulation. Proc. Nat. Conf. Artif. Intell. (AAAI-93) 1993, 606-613. Kay, H.; Ungar, L. H. Deriving Monotonic Function Envelopes from Observations. Working Papers from the Seventh International Workshop on Qualitative Reasoning about Physical Systems, Rosario Resort, Washington, 1993; pp 117-123. Knorr, R. S.; O’Driscoll, K. F. Multiple Steady States, Viscosity, and High Conversion in Continuous Free-Radical Polymerization. J. Appl. Polym. Sci. 1970, 14, 2683-2696. Kuipers, B. J. Qualitative Simulation. Artif. Intell. 1986, 29, 289338. Kuipers, B. J. Qualitative Reasoning: Modeling and Simulation with Incomplete Knowledge. Automatica 1989, 25, 571-585. Kuipers, B. J.; A° stro¨m, K. Composition and Validation of Heterogeneous Control Laws. Automatica 1994, 30 (2), 233-249. Kuipers, B. J.; Shults, B. Reasoning in Logic about Continuous Systems. In Principles of Knowledge Representation and Reasoning: Proceedings of the Fourth International Conference (KR-94), Morgan Kaufmann, San Mateo, CA, 1994. Marten, F. L.; Hamielec, A. E. High Conversion Diffusioncontrolled Polymerization. ACS Symp. Ser. 1979, 104, 43-70. McGreavy, C., Ed. Polymer Reactor Engineering; Blackie Academic & Professional, VCH Publishers: New York, 1994.

Moon, I.; Powers, G. J.; Burch, J. R.; Clarke, E. M. Automatic Verification of Sequential Control Systems using Temporal Logic. AIChE J. 1992, 38 (1), 67-75. Ogunnaike, B. A.; Ray, W. H. Process Dynamics, Modeling, and Control; Oxford University Press: Oxford, U.K., 1994. Penlidis, A.; Ponnuswamy, S. R.; Kiparissides, C.; O’Driscoll, K. F. Polymer Reaction Engineering: Modelling Considerations for Control Studies. Chem. Eng. J. 1992, 50, 95-107. Ross, R. T.; Laurence, R. L. Gel Effect and Free Volume in the Bulk Polymerization of Methyl Methacrylate. AIChE Symp. Ser. 1976, 72, 74-79. Schmidt, A. D.; Ray, W. H. The Dynamic Behavior of Continuous Polymerization Reactors. I. Isothermal Solution Polymerization in a CSTR. Chem. Eng. Sci. 1981, 36, 1401-1410. Shamma, J.; Athans, M. Analysis of Gain-Scheduled Control for Nonlinear Plants. IEEE Trans. Autom. Control 1990, 3518, 898. Soroush, M.; Kravaris, C. Multivariable Nonlinear Control of a Continuous Polymerization Reactor: An Experimental Study. AIChE J. 1993, 39 (12), 1920-1937. Zafiriou, E.; Morari, M. Robust Process Control; Prentice-Hall: Englewood Cliffs, NJ, 1989.

Received for review July 14, 1995 Accepted February 7, 1996X IE9504361

X Abstract published in Advance ACS Abstracts, May 15, 1996.